WO2025160583A1 - Tunable shear wave front generation for biomechanical measurement of biological tissue - Google Patents

Abstract

Some aspects of the described technology may include a method, including inducing a wave field in a biological sample, the wave field comprising waves having a wave frequency, scanning the sample at a scanning velocity to perform a plurality of spatial measurements of the sample, determining spatial data from the plurality of spatial measurements, amplitude demodulating the spatial data based on the wave frequency and the scanning velocity, and generating displacement data for the biological sample based on the amplitude demodulated spatial data.

Description

TUNABLE SHEAR WAVE FRONT GENERATION FOR BIOMECHANICAL MEASUREMENT OF BIOLOGICAL TISSUE

CROSS-REFERENCE TO RELATED APPLICATIONS

[0001] This application claims the benefit of US Provisional Application No. 63/625,677, filed January 26, 2024, and US Provisional Application No. 63/710,475, filed October 22, 2024, the contents of each of which are incorporated in their entirety.

STATEMENT REGARDING FEDERALLY SPONSORED RESEARCH OR DEVELOPMENT

[0002] This invention was made with government support under award numbers 5R01EB033306-02 and 5P41EB015903-12 (agreements 2021A011939 and 2021A012865) from the National Institutes of Health and the National Institute of Biomedical Imaging and Bioengineering. The government has certain rights in the invention.

BACKGROUND

[0003] Optical coherence tomography (OCT) has enabled high resolution imaging of the subsurface micro-structure of tissues. By leveraging interferometry to measure light reflected from a sample, OCT can acquire depth-resolved cross-sectional images with a lateral resolution of a few microns and penetration depth of a few millimeters. Due to the transparency of ocular tissue, OCT has seen the most translation in ophthalmology; however, its compatibility with fiber optics makes it also well suited for imaging internal organs through catheters or endoscopes.

[0004] Beyond ophthalmology, there are many other disease pathways and clinical targets that also concur with biomechanical stiffness, such as tumors due to fibrosis, pulmonary fibrosis due to thickening of alveolar walls, and atherosclerosis due to heterogeneous plaque composition and calcification.

[0005] To image stiffness, optical coherence elastography (OCE) is a functional extension of OCT which tracks internal displacements due to external forces. These forces range from quasi-static compression, to passive excitation from natural bodily motions such as the heartbeat, all the way to dynamically actuated shear waves. The measured displacements that occur as a result of external forces are then correlated to a mechanical model, from which parameters such as elastic modulus or shear modulus can be derived.

[0006] Wave-based OCE methods leverage the properties of shear waves traveling through tissue primarily to extrapolate shear modulus. It may be advantageous for imaging in vivo because, unlike quasi-static compression methods, it does not require sensitive reference measurements or computationally-intense finite-element modeling. On the other hand, passive shear wave excitation from natural bodily motions such as the heartbeat is may also be too slow and uncontrollable depending on the tissue target of interest. Comparatively, reverberant elastography is a wave-based method which uses multiple contact points of excitation. Instead of trying to avoid shear waves traveling in all directions, internal reflections from tissue boundaries within a sample are beneficial towards generating a diffuse wave field. A diffuse wave field is, in essence, an objective shear wave speckle field produced by the coherent superposition of fully randomized shear waves. Then, correlations of the resulting velocity field— which represent an underlying relationship between shear wave speckle size and shear wavelength— enable measurements of localized shear wave number. Given the frequency of excitation used to generate shear waves, one can directly calculate the shear wave speed and shear modulus from the shear wave number.

SUMMARY

[0007] Some aspects of the described technology may include a method, including: inducing a wave field in a biological sample, the wave field comprising waves having a wave frequency; scanning the sample at a scanning velocity to perform a plurality of spatial measurements of the sample; determining spatial data from the plurality of spatial measurements; amplitude demodulating the spatial data based on the wave frequency and the scanning velocity; and generating displacement data for the biological sample based on the amplitude demodulated spatial data.

[0008] Further aspects may include a system, including: a driver to induce a wave field in a biological sample, the wave field comprising waves having a wave frequency; an imaging system to scan the sample at a scanning velocity to perform a plurality of spatial measurements of the sample; a processor; and a non-transitory computer readable medium storing instructions executable by the processor to: determine spatial data from the plurality of spatial measurements; amplitude demodulate the spatial data based on the wave frequency and the scanning velocity; and generate displacement data for the biological sample based on the amplitude demodulated spatial data

[0009] Further aspects may include a non-transitory computer readable medium storing instructions executable by a processor to control a driver to induce a wave field in a biological sample, the wave field comprising waves having a wave frequency; receive sample scan data comprising a plurality of spatial measurements of the sample from an imaging system having a scanning velocity; determine spatial data from the plurality of spatial measurements; amplitude demodulate the spatial data based on the wave frequency and the scanning velocity; and generate displacement data for the biological sample based on the amplitude demodulated spatial data.

BRIEF DESCRIPTION OF THE DRAWINGS

[0010] FIGS. 1A-1G illustrate example elastography systems.

[0011] FIGS. 2A-2C illustrate examples of spatial data amplitude demodulation.

[0012] FIGS. 3A-3B illustrate example elastography systems with indirect contact-induced wave fields.

[0013] FIGS. 4A-F illustrate an example of spatial data amplitude demodulation via Fourier analysis.

[0014] FIG. 5 illustrates example Bessel function sums from different wave field configurations.

[0015] FIG. 6 illustrates an example system that may be used to implement an elastography system.

[0016] FIG. 7 illustrates further aspects of an example system that may be used to implement an elastography system.

[0017] FIGS. 8A-8E illustrate example aspects of an experiment demonstrating the described technology.

[0018] FIGS. 9A-9I illustrate example aspects of an experiment demonstrating the described technology.

[0019] FIGS. 10A-10C illustrate example aspects of an experiment demonstrating the described technology.

[0020] FIGS. 11A-11F illustrate example aspects of an experiment demonstrating the described technology.

[0021] FIGS. 12A-12F illustrate example aspects of an experiment demonstrating the described technology.

[0022] FIGS. 13A-13F illustrate example aspects of an experiment demonstrating the described technology. [0023] FIGS. 14A-14E illustrate example aspects of an experiment demonstrating the described technology.

[0024] FIGS. 15A-15D illustrate example aspects of an experiment demonstrating the described technology.

[0025] FIGS. 16A-16G illustrate example aspects of an experiment demonstrating the described technology.

[0026] FIGS. 17A-17C illustrate example aspects of an experiment demonstrating the described technology.

DETAILED DESCRIPTION

[0027] Despite its utility, existing OCT systems lack the ability to robustly and conveniently measure biomechanical contrast— a feature highly correlated with ocular pathologies such as keratoconus, intraocular pressure dysregulation, and glaucoma.

[0028] Compared to other wave-based methods, which use a single point of excitation and a single wavefront to measure shear wave speed along one dimension, reverberant elastography relies on a diffusely generated shear wave field and is therefore theoretically more robust to motion. It is also conceptually a 3D volumetric technique, although it requires large correlation windows on the order of the shear wavelength— typically, >1 mm— and thereforeO exhibits subpar resolution along the correlation dimensions— generally corresponding to the two lateral dimensions.

[0029] Despite its overall benefits, reverberant elastography still has certain limitations that hamper its widespread clinical translation. In particular, since shear waves travel much faster than conventional OCT raster-scanning speeds, state-of-the-art wave-based elastography methods require reproducible synchronization (or pseudo-synchronization) of shear wave excitation with imaging to acquire MB-mode scans (a B-scan composed of M-mode scans at each lateral position) at all spatial locations. In the cornea, reverberant elastography has been demonstrated with a 100x100 lateral sampling volume necessitating an acquisition time of 1 minute. The nature of the MB-mode scan implies that motion during the entire volumetric acquisition will destroy the coherency of the diffuse shear wave field and thus hamper the analysis. This severe limitation may explain the lack of a demonstration of reverberant elastography in vivo. It is also crucial to understand that merely having a faster A-line rate does not necessarily enable shorter acquisition times: the technique relies on observing a number of periods (at least one but more commonly up to 8 periods) of the shear wave field, which is a function of the shear wave frequency, not the A-line rate. Thus, current limitations in reverberant elastography makes this method impractical for clinical translation.

[0030] The stiffness and compliance of biological tissues are vital characteristics that can be altered in the course of pathologic change. There is increasing evidence that elasticity plays a key role in corneal pathologies such as keratoconus, myopia, and complications from refractive error corrective surgery. Optical coherence tomography is a widely established ophthalmic imaging technology. Along with virtually every other conventional medical imaging modality, it has been adapted to tissue elasticity measurements. One promising approach measures the velocity of shear waves propagating through tissue to quantify traditional elasticity metrics, although its widespread clinical adoption has been complicated by slow image acquisition, sensitivity to motion, expensive custom hardware, and requirement for contact with the cornea. Conventional scanning protocols using OCT for displacement measurements, (such as continuous raster-scanning) typically use arduous phase-locked MB-mode scanning protocols. However, these methods are slow and more sensitive to motion.

[0031] Disclosed herein is an apparatus for tunable elastic wave front generation in biological tissue to enable rapid, real-time, in-vivo quantification of bio-mechanical properties. First, elastic waves are generated by an input energy source which causes displacement on or within the biological tissue. The geometry, frequency or multiple frequencies, amplitude, and duration of the applied displacement may all be tunable to generate a desirable wavefront.

[0032] Aspects of the described technology may address these challenges by enabling rapid and robust measurements of fully coherent 3D shear wave displacement fields in human subjects in vivo. By continuously raster-scanning, sensitivity to motion may be reduced while improving lateral sampling by over an order of magnitude compared to other OCE methods. An autocorrelation-based phase gradient method for shear wave number retrieval may be insensitive to the makeup of the shear wave field. Indirect piezoelectric contact via the eyelid may be used to generate random shear wave fields in the cornea, without the need for careful alignment or focusing. The described technology supports elasticity measurements of the cornea in subjects with a conventional frame-rate OCT system, all without the need for custom lasers or additional complex hardware. These results open the door to many appealing clinical applications in need of higher fidelity imaging tools ranging from detection of keratoconus, to screening for refractive-error correction surgeries, to understanding and preventing myopia in pediatric and adolescent populations.

[0033] Certain examples may be described herein with respect to a right-handed coordinate system defined such that the optical axis or depth direction is z (increasing in depth), the inplane lateral dimension (e.g., direction of scans) is x, and out-of-plane lateral dimension (e.g., direction of translation between scans) is y. The described technology is not limited to or by any particular coordinate systems, any suitable coordinate system may be employed or the described operations or data may be performed or provided in a coordinate-free manner.

[0034] FIGS. 1A-G illustrate an example system 100 for tunable elastic wave front generation for biomechanical measurement of tissue. System 100 may include an energy source to provide motive power from which displacements in a sample will originate. The energy source is coupled into a sample 105, generating internal displacements in the sample 105 that can be measured in order to characterize its biomechanical properties. System 100 may further include an imaging system to measure displacements induced in the sample. Embodiments of an imaging system that could be used to for displacement measurements may include but are not limited to optical coherence tomography (OCT), full-field OCT, linefield OCT, holography, high speed cameras, ultrasound, MRI, etc.

[0035] In the illustrated example, a displacement field is induced in a sample 105. Sample 105 may comprise an in-vivo location on/in a subject's body (e.g., a human subject or a nonhuman subject) or an in-vitro sample. Examples of sample 105 include, but are not limited to, lung tissue, an eye, skin, blood vessels, gastrointestinal tracts, internal organs, bone structures, etc. While sample 105 is illustrated, for ease of illustration, samples 105 in applications may have any geometry.

[0036] System 100 may include a motive energy source 101, 110 from which displacements may originate. In some cases, the motive source may include a single actuator 101 to generate displacements originated from a single location. In other cases, such as illustrated, the motive source may include multiple actuators devices. In some cases, various aspects of the motive energy source, such as the number of may be configurable. For instance, system 100 might be used with two actuators 101, 110 to measure a first sample 105 and might be used with three actuators to measure another sample 105, etc. [0037] In the illustrated example, source 101, 110 may comprise a plurality of piezoelectric actuators 101, 110. In such examples, a piezoelectric actuator 101, 110 may take the form of a bender, stack, ring, arc, or any othertunable geometry depending on the desired conditions of the elastic wave field and the corresponding sample of interest. For example, a ring may generate focused elastic waves, whereas a piezoelectric stack generates a relatively lower amplitude and normal elastic force while also maintaining a relatively smallerfootprint. Other actuators to generate displacements may include motors (such as those with an eccentric mass), microelectronic mechanical system (MEMs) actuators, servos, solenoids, pneumatic cylinders, air puff generators, sound/ultrasound transducers (e.g., to generate sonic or ultrasonic pressure waves, acoustic radiation force, cavitation, etc.), lasers (e.g., via a photothermal effect), cavitation devices, etc. System 100 may include one or more of any of such energy sources in any combination, each of which may also be driven independently.

[0038] System 100 may further include contacts 103, 108, 111, 118 to physically contact sample 105 and induce a wave field in sample 105 based on motion of actuators 101, 110. The illustrated example, piezoelectric actuators 101, 110 may comprise bending actuators 101, 110. The bending of the actuators results in motion in the z-axis (e.g., depth axis) generating elastic waves originated at contact locations 102, 109, 115, 116. In further examples, contacts 103, 108, 111, 119 may be coupled to other actuators as described above, such as eccentric motors, MEMs actuators, servos, solenoids, gas cylinders, etc.. Contacts 103, 108, 111, 119 might be mechanically coupled to an actuator, electronically coupled to an actuator, pneumatically coupled, etc. For instance, an air-puff generator may be positioned to actuate contacts 103, 108 by blowing on them from above (e.g., the illustrated upper planar surface), etc.

[0039] In some examples, contacts 103, 108, 111, 118 may be interchangeable to provide a chosen wavefront geometry. For instance, contacts 103, 108, 111, 118 may be dimensioned to provide a contact location geometry, such as a line, square, rectangle, circle, blunt point, arc, grid, speckle, or any othergeometry. In the illustrated example, contacts 103, 108 provide straight line contact locations 102, 109 and contacts 111, 118 provide arcuate contact locations 115, 116. System 100 may support diverse contact geometries to generate a wavefield which is more chaotic than an individual wavefront, yet still well-defined in view of known excitation geometry. For example, in an embodiment where contact locations 102, 109 have a surface area with a linear edge (such as line, rectangle, square, or otherwise), plane waves may be generated from these geometries, which may support 2D processing of the displacement field. In an embodiment where a contact location is a blunt single point, spherical waves may be generated for analysis, which may support ID analysis of the displacement field among other strategies.

[0040] System 100 may include a single contact 101 to generate a displacement field originated from a single location. System 100 may further include plural contacts 103, 108, 111, 118 having multiple geometries (such as multiple points, two lines parallel to one another) to intentionally generate a superposition of their respective displacement fields. In a specific embodiment where two or more lines generate elastic waves in tissue, a semi- reverberant field is generated, which enables both heterogenous tissue characterization and high-speed imaging and processing.

[0041] As another example, FIG. ID illustrates an example system 100 comprising an internal examination device 119, such as an endoscope, bronchoscope, catheter, colonoscope, etc. For example, device 119 may be used to examine samples such as an inner or outer ear, lung or respiratory tract (e.g., nose, mouth, trachea, bronchi, etc.), urinary tract, gastrointestinal tract, etc.

[0042] In this example, actuator 118 may comprise an actuator to induce the displacement field via radiative coupling. For example, actuator 118 may comprise an acoustic transducer (e.g., ultrasonic or sonic) and lens/lens array. In this example, acoustic waves may also be absorbed into tissue to generate shear waves. One or more lenses in an array may be used to focus acoustic energy in a controlled geometry onto or into the biological sample.

[0043] As another example, actuator 118 may comprise a light/heat source for photothermal excitation. Actuator 118 may further include use a spatial light modulator or diffractive optical elements to produce tunable and controllable shear wave fronts in biological tissue. Absorption of photothermal energy in tissue causes displacement and shear wave excitation. Actuator 118 may further include one or more waveguides to provide a controlled delivery of photothermal radiation to induce displacements with known or controllable properties.

[0044] As a further example, actuator 118 may comprise a vibration source, such as a piezoelectric actuator, motor, servo, MEMS vibrator, etc. As anther example, actuator 118 may be located externally to device 119 (e.g., in unit 117, etc.) and device 119 may induce vibrations in the sample via mechanical force transfer from the external actuator. As another example, actuator 118 may be placed in contact with an externally accessible region of a subject while device 119 images the internal sample. For instance, an actuator such as ring actuator 124 may may generate displacements in distal lung tissue via coupling of an energy source with other more proximal tissues, such as bronchi, which are easier to access from outside of the body. In a similar embodiment, displacements may also be generated indirectly in distal lung tissue through excitation of rib bones which are also easier to access from outside of the body.

[0045] FIG. IE illustrates an example system 100 for elastography of a sample 105 comprising a blood vessel (e.g., artery, vein, etc.). In this example, actuator 121 may comprise a cavitation source to cavitate blood within a blood vessel 105. For instance, actuator 121 may comprise a photothermal source, an acoustic source, a mechanical source, etc. In this example, actuator 121 may cause excitation of blood which is in contact with the vessel wall. Expansion of blood within the vessel will generate a wave field 123 in the vessel wall.

[0046] In other examples, an actuator may induce a displacement field in the sample 105 directly. For example, as illustrated in FIG. IF, a piezoelectric actuator package 124 (e.g., a plastic housing containing an actuator) may be in direct contact with a sample 105 such as on the skin of a human subject. In such examples, the shape of the contact area from the package may be as described with respect to contacts 103, 108, 111, 118 to provide any geometry. For instance, in the illustrated example, piezoelectric actuator package 124 comprises a ring, for example to generate focused elastic waves within the circumference to the ring.

[0047] FIG. 1G illustrates an example system 100 comprising a plurality of actuators 125, 126, 127, 128. In some examples, actuators 125, 126, 127, 128 may be driven at different frequencies to induce a complex wave field in the imaging field of view 104. Imaging at multiple different frequencies simultaneously may speed data acquisition time (reducing the time for a patient to have to remain still or in an uncomfortable position, for example) and may also provide information not readily available via independent scans at different frequencies. For instance, each respective actuator 125, 126, 127, 128 may be driven at one of four respective frequencies (see, e.g., FIG. 16, see also FIG. 15). As another example, each actuator 125, 167, 127, 128 may each be driven independently by a signal comprising multiple frequencies. As another example, actuators 125, 126, 127, 128 may be driven at the same frequency but with different phase offsets or amplitudes. Implementations may employ any combination of wave inducement parameters. [0048] System 100 may include a controller 117 which may include a driver 114. For example, driver 114 may include circuitry to generate a driving signal for the motive energy source, such as analog circuitry, discrete components, application specific integrated circuitry (ASIC), a processor (e.g., a digital signal processor, general purpose processor, embedded microprocessor, etc.) with computer readable medium storing instructions to generate the driving signal, a field-programmable gate array FPGA, etc. The driving signal may determine the magnitude of the energy delivered into the sample over time. In one embodiment, the source 101, 110 may be driven by a single arbitrary frequency (e.g., 2 kHz) or multiple frequencies at the same time. Parameters such as the amplitude, phase, wave shape, duty cycle, etc. of these frequencies may also be tuned. For examples, parameters may be varied energy transfer into the tissue, to match other system parameters, including but not limited to imaging system frame rate, etc., to vary the displacement field to avoid standing waves, etc.

[0049] These signals may generate a variety of elastic waves in the sample, including compression waves, elastic waves, and longitudinal waves. Source 101, 110 may also be driven by a short pulse or many short pulses, one after another. For example, source 101, 110 may be driven by transient frequency pulses (e.g., lasting a second, millisecond, etc. or lasting a certain number of waves) separated by relaxation periods or providing a low frequency signal envelope. In another embodiment, the source may be driven slowly in a static or quasi-static fashion to induce static tissue displacements (e.g., such as via elastic waves with wavelength larger than the tissue sample of interest). In further embodiments, any one or more of these driving signals may also be stitched together consecutively in time or simultaneously to form a custom sequence. In examples comprising multiple displacement energy sources 101, 108 different sources may be driven independently by any of the driving signals described herein. In some examples, driver 114 may output a tunable driving signal, which may be determined/selected via software executed by controller 117.

[0050] System 100 may include an imaging system 106, 113 to detect the displacements that are generated as a result of the energy transfer between source 101, 110 and sample 105. In various embodiments, the imaging modality may be independent of the energy source or sources, driving signals, tissue-coupling, or displacement source geometry. As non-limiting examples, imaging system 106, 113 may include OCT, full-field OCT, line-field OCT, holography, high speed cameras, ultrasound, MRI, etc. [0051] In the illustrated example, imaging system 106, 113 may comprise an OCT imaging system. For example, an OCT imaging system may comprise a light delivery/reception unit 106, such as a galvanometer and scan lens, and an energy generator or sensor 106, such as a laser source (e.g., a vertical cavity surface emitting laser (VCSEL), etc.) and an interferometer. Here, the system may translate the light source in the y-direction after completing a final A- line scan of a particular B-plane. In some examples, such as the one illustrated, imaging system 106, 113 may perform a raster scan. For instance, in the illustrated example, imaging comprises a plurality of axial-aligned (A-line) scans 107 over a plane 104 (e.g., a brightness plane (B-plane)) aligned with the x-z plane.

[0052] System 100 may include any components for displacement detection. For instance, in an example where the imaging system comprises an OCT or other interferometry-based imaging system, displacement measurements may be detected with Doppler-OCT, speckle correlation, or measurements of the phase of the electric field and how it changes between acquisitions. As another example, in embodiments that use a high-speed medical imaging camera, displacements may be detected via digital image correlation. As a further example, embodiments that use ultrasound for displacement detection, displacements may be quantified with cross-correlations of the measured field. Of course, any other suitable imaging modality may be employed and, in some cases, different imaging modalities may be combined. For instance, an implementation might include an OCT device 106/113 and a highspeed camera device.

[0053] System 100 may further comprise a processor 112 to execute instructions stored on a non-transitory computer readable medium to determine a displacement field from the measured field. For example, data from multiple scans of a scan plane location (e.g., plane 104) may be combined to determine the displacement of the wavefield across a time period. As an example, depth data from two consecutive scans may be compared to determine a displacement (e.g., in the z-direction) of the wavefield across the time period.

[0054] In some implementations, the spatial data of the wave field captured by imaging system 106/113 may include distortion, noise, artifacts, etc. For example, in an asynchronous, scanning image capture modality, the wave field may change between scan samples. For instance, if scanning-rate is relatively slow compared to the temporal evolution of the displacement field higher frequency signal distortions may occur. For example, FIG. 2A illustrate an example of output data of a BM-mode OCT image of a displacement field including distortions from A-line scanning in the x-direction and distortions from successive B-line scans in the y-direction.

[0055] For example, as a displacement field, particle displacement S(E, t) of a elastic wave field as a function of position E (x, y, z) at time t may be described as follows (see, FIG. 4A):

Reverberant shear waves are modeled as the superposition of plane waves propagating in many random directions. The particle displacement S(E, t) of the shear wave field as a function of position E (x, y, z) and time t is therefore defined by Eq. (1) where the propagation direction of each wave is represented by n_q and is defined in the spherical coordinate system by angles O and 0. The imaginary constant is j . The polarization n_qi of each wave is normal to the direction of propagation, so it is defined as orthogonal to n_q by a third angle, a. The amplitude and phase of each wave are represented by s_qi, while k is the shear wave number and wo is the shear wave frequency.

[0056] However, when performing a scanning imaging protocol, there is a time delay between each scan capture. For instance, in raster-scanning BM-mode asynchronous imaging, t can not be considered constant within a B-scan when elastic wear speeds are too fast relative to raster-scanning speeds. The non-negligible time delay between each A-line distorts measurements. Implementations may account for such distortion by identifying a carrier frequency term introduced by the time delays. For example, in the OCT example above, for the ith OCT frame, time may be represented as ti + AE/v_SCan where AE is the distance away from the initial A-line and v_SCan is the raster-scanning speed. Accounting for the time delay induced by scanning, Eq. (1) becomes:

The real component of this definition of the displacement compared to previous methods [Eq. (1)] reveals that a difference between synchronous and asynchronous imaging in the term cos (woAE/v_SCan). This term introduces a modulation across each scan by a carrier frequency that is based on the wavefield frequency and the scanning velocity: wo/v_SCan. In an implementation where the wave field comprises a superposition of waves having different frequencies Wi, multiple carrier frequencies Wi/v_SCan may be present in the data. In some examples, processor 112 may amplitude demodulate the spatial data based on a carrier frequency. For instance, FIG. 2B illustrates an example of spatial data after demodulation to remove the carrier frequency. Accordingly, here, distortion remains in the y-direction from delay between B-scans, but distortion within each B-scan (e.g., each constant y slice of x-data) has been removed. While illustrated with respect to a single excitation frequency, a similar procedure may be applied when multiple different frequencies are used simultaneously to excite the wave field. At a common scan frequency, Vscan, each wave frequency w; corresponds to a carrier artifact at Wi/v_Scan. Any techniques to demodulate overlapping multiple amplitude-modulated signals may be applied to remove such artifacts.

[0057] In some examples, after demodulation, the resulting spatially coherent spatial (e.g., displacement) field will be complex-valued. The real component gives the displacement field at one point in time and the phase may be used to produce the entire harmonic elastic wave field. The complex-valued displacement field contains all the information needed reconstruct the entire wave field at any point in time. The resulting (now spatially coherent) displacement field is complex. The real component gives the displacement field at one point in time and the complex phase may be used to produce the entire harmonic shear wave field. The analytical representation of displacement contains all the information needed to generate the entire wave field at any point in time. Multiplying the complex field by a corresponding exponential phase term produces the shear wave field at a different point in time.

[0058] To produce the elastic wave field at a different point in time, e.g. t + At, processor 112 may multiply the complex displacement field by an exponential phase term e^J^^{t >0}, where At, is time offset and phase difference between the elastic wave field recovered at each out- of-plane y-location. In some cases, At may be determined based on imaging parameters. For example, in an OCT example, At may be calculated based the OCT sampling rate and number of samples per frame. For a time-delay of At_y between y-locations (based on the laser repetition rate and sampling parameters), the expected phase shift is e^l&ty^ai°. Processor 112 may then compute an expected phase shift for each y-location. Cumulatively multiplying this phase shift to each xz-displacement field may correct for artifactual phase jumps of the elastic wave field between each y-location. For example, FIG. 2C illustrate the results of results of phase noise correction of the demodulated spatial data illustrated in FIG. 2B. Thus, FIG. 2C illustrates an entire coherent 3D elastic wave speckle field recovered from the original distorted data of FIG. 2A.

[0059] In further examples, processor 112 may perform other steps to remove noise or account for other sources of error/artifacts. For instance, processor 112 may perform surface wave correction, tissue flattening, Gaussian noise reduction, etc.

[0060] For example, in a system 100 including a galvanometer scanner (galvos) 106 for scanning, the galvos might not scan completely orthogonal to each other, which may introduce artifacts. For example, the data may include wave artifacts that tend to propagate spuriously along y due to an additional phase delay. In some examples, processor 112 may perform a calibration procedure to remove such artifacts. Calibration may depend on factors such as a sample being imaged, parameters of a particular galvos, distance to the sample, field of view, etc.

[0061] As another example, processor 112 may apply motion correction to the data. For example, processor 112 may apply motion correction to account for subject movement between scans of different lines, to account for subject movement between scans of the same line, etc. For example, to address motion-induced speckle decorrelation, processor 112 may apply zx bulk motion correction. As another example, to address axial motion between scan repetitions, processor 112 may apply non-rigid axial motion correction on scan data between repetitions (e.g., between B-scan repetitions) in each y-location independently. As another example, to address macroscopic deformation of the elastic wave field due to motion, processor may perform bulk motion correction, such as zx bulk motion correction between B-scan pairs acquired at each y-location. In some cases, processor 112 may perform such motion correction in view of subject topology. For instance, the description of the experiments below include a description of motion correction while accounting for corneal topology. Processor 112 may perform various techniques, such as upsampling, to computer motion detection/correction at a subpixel level. In some examples, processor 112 may perform motion correction to the data prior to generating displacement data (e.g., the processor may determine the displacement field discussed above using motion corrected data). As another examples, processor 112 may apply surface wave correction techniques to the data.

[0062] Returning to wave generation in a biological sample, the described technology may provide tunable wave generation. For example, waves induced in a sample may comprise transverse (e.g., shear) waves, longitudinal waves, combinations thereof. As discussed above, various components may be interchanged or configured to provide a tunable shear wave front generation for biomechanical measurement of biological tissue. As an example, FIGS. 3A-3B illustrate an example imaging system 300 including a wave field generated by indirect wave transmission to a biological sample.

[0063] As illustrated in FIG. 3A, system 300 may comprise an energy source 302 to induce a wave field 310 in a sample 308. Energy source 302 may comprise any suitable source of displacement in a biological sample, such as, for example, a piezoelectric actuator, a sonic or ultrasonic pressure transducer, an acoustic radiator, a laser or other photothermal excitation source, a motor, a servo, an air blower, cavitation source, etc. In this example, energy source

302 is coupled to a contact 303. As discussed above, such as with respect to contacts 103, 108, 111, 118, contact 303 may have any suitable form factor or geometry to provide an origin location for the wave field 310, such as, for example, any number of prongs, lines, blocks, wedges, rings, arcs, or any other geometry. For instance, in the illustrated example, contact

303 comprises blunt tipped prongs 305.

[0064] System 300 may further include a driver 301 to generate a driving signal for energy source 302. For example, the driving signal may determine the wave frequency or frequencies of the wavefield 310. As an example, driver 301 may induce a vibration in contact 304 comprising multiple simultaneous frequencies, such as through a frequency comb, etc. For instance, system 300 may generate waves of multiple distinct frequencies (2 Khz, 4 Khz, 6 Khz, etc.), and then use raster-scanning to individually separate out the displacement field from each one. This may support non-linear frequency dependent measures of shear modulus, or viscoelasticity. Since most tissue is actually viscoelastic (not elastic, as many methods tend to assume), probing multiple frequencies may allow for a relatively more accurate quantitative measure of tissue mechanical properties. As another example, driver 301 may drive source

302 with a frequency synchronous with an imaging frame rate period for controllable amplitude and phase changes as a function of frame rate.

[0065] System 300 further comprises an imaging system 303. For example, imaging system

303 may comprise an OCT imaging system, an ultrasound imaging system, an MRI, or any other suitable imaging system. Here, imaging system 303 is oriented to capture data from sample 308 having a wave field 310 generated via indirect contact by contact 304. In the illustrated example, a wavefront 309 is generated in a neighboring body part 307 to sample 308. Indirectly induced shear waves may travel through multiple boundaries and changes in tissue impendence, leading to a superposition of one or more waves traveling in random directions. As the waves 309 cross the interface 306 between the sample 308 and the neighboring part 307, a wave field 310 is generated comprising a superposition of waves propagating in multiple directions. As described herein, the measurement of the displacement field in sample 308 may be robust to waves traveling in various directions. Accordingly, this may provide for indirect inducement of waves in sample 308 that otherwise would require direct contact, which may be uncomfortable, impractical, etc.

[0066] FIG. 3B illustrates an example of system 300 applied to perform elastography imaging of a subject's cornea 311. In this example, contact 312 comprises a blunt round tip, which is placed against the subject's eyelid 314. A driving signal provided by driver 301 induces a wave field 313 that is coupled to the cornea 311 via the body of the subject's eye 315 providing a direct path to the cornea 311. While not illustrated (for ease of explanation), wave field 313 may comprise a superposition of waves generated via reflection, multiple transmission paths, indirect vibration of the surrounding skull, etc. As illustrated, contact 312 may be manual pressed to the eyelid 314 (e.g., by the subject themselves or by a practitioner), or contact 312 may be positioned via a headset, headrest mounted system, etc.

[0067] In further examples, tissue displacements may be generated in the anterior segment of the eye (including cornea, sclera, limbus, lens) indirectly via energy source coupling around the eye rather the eye itself. Areas around the eye which may be used for this type of indirect coupling of tissue displacements include but are not limited to the upper eyelid, lower eyelid, nose, browbone, temple-region, teeth, ears, and skin. Tissue displacements from any of these regions may also be generate waves anywhere within the eye, beyond the anterior segment, including the posterior eye (retina, optic nerve head). As the described technology is robust to random variations in a wave field, in some cases, a practitioner or user may be able to place a contact 312 on various such locations (e.g., selected for comfort or convenience), without impacting data processing steps.

[0068] As discussed above, systems may perform amplitude demodulation on imaging data to remove signal artifacts from scanning. Implementations may employ any method of amplitude demodulation to remove such artifacts. As an example, FIGS. 4A-F illustrate aspects of artifact removal via Fourier analysis. [0069] FIG. 4A illustrates spherical coordinates and a unitary polarization vector of an example wave. In particular, FIG. 4A illustrates certain geometric parameters discussed above with respect to Eqs. (1) and (2). Here, the right-handed coordinate system is defined such that the depth direction is z (increasing in depth), in-plane lateral dimension is x, and out-of-plane lateral dimension is y. Some implementations, such as phase-sensitive Doppler OCT, may be sensitive to displacements only along the z axis, while others may be sensitive to displacements in multiple directions (e.g., along multiple axes, such as zx, etc.) As illustrated in FIG. 4B, effects of scanning may be modeled as the modulation 402 of the signal of interest 401 by a carrier frequency (e.g., a cosine term 404 or other carrier signal representation). The resultant modulated signal 405 is demodulated 403 to recover the signal of interest 406 (e.g., a spatially-coherent displacement field).

[0070] FIGS. 4C-F illustrate various stages of demodulation 403. When modulated by a cosine, the signal of interest 407 (FIG. 4C) is shifted to both the positive and negative carrier frequency 408, 309 in the Fourier domain (FIG. 4D). To recover the analytical signal of interest, the Fourier transform of the modulated signal may be computed (FIG. 4E) and the spectrum may be down shifted by the carrier frequency such that one side 411 of the desired signal is centered at zero and a bandpass filter may be applied to suppress the second mirrored copy 410. The inverse Fourier transform of the isolated signal 412 returns to the spatial domain and provides a demodulated signal 406.

[0071] In some examples, wave fields are generated such that waves propagate aligned with the scan plane. For example, in FIG. 1A, a wave field might consist of a superposition of waves traveling in positive x direction and the negative x direction. Such aligned waves may propagate along the z axis or may be confined to propagate only in the x axis. As discussed below, this semi-reverberant wave field may support measurement of the field via autocorrelation in the scanning x-axis. Of course, while discussed with respect to shear waves and Doppler OCT, these techniques may be used with any wave field configuration and imaging modality. Additionally, while the below focuses on spatial correlation with respect to the raster-scanning axis x, a similar derivation follows in all dimensions. Modeling shear waves as complex exponential plane waves, the displacement field and velocity field are mathematically identical beyond a scaling factor and n/2 phase shift. Given the definition of the shear wave field in Eq. (1), the normalized ensemble average of the autocorrelation of the displacement yields an analytical expression that is a function of shear wave number k. In other words, the expected value of the shear wave speckle size depends only on shear wavelength for a fully randomized superposition of waves.

[0072] In the fully reverberant case (e.g., a diffuse random wave field), this ensemble average may be computed by an integral representing waves traveling randomly in all directions and polarizations, which corresponds to a uniform spherical integral over all angles of (D, 0 and a. In a typical application, the autocorrelation is computed individually at each time step only with respect to E, so At = 0. This results in a sum of spherical Bessel functions: where S_z is the Doppler phase difference measured displacement along the sensor axis and x is the B-scan axis along which the beam is scanned. While the normalized autocorrelation for both the displacement and velocity fields are identical, a benefit is that displacement measurements do not require temporal coherence. This may drastically reduce frame rate requirements over other approaches that require MB-mode scanning and shear wave synchronization.

[0073] As discussed above, in some cases out-of-plane y-dimension spatial coherence may be recovered via phase shifting B slice data. In some cases, such as where such techniques are impractical or unavailable, autocorrelation along the scanning x-axis may result in identically sized waves traveling at different angles relative to the B-scan plane appearing to have completely different. While in some cases azimuthal averaging addresses this problem by sampling all possible angles, this might not be practical or available when the autocorrelation along y is not accessible. This may be addressed via a semi-reverberant excitation scheme. A semi-reverberant excitation scheme may generate waves whose angle <p with respect to the B-scan plane is known. For instance, waves may be mechanically excited that only travel in-plane within B-scan planes, that travel with constant angles across the B-scan planes, or that have arbitrary angles <p that are known as a function of at least one position coordinate. For instance, <p(x) may represent an arbitrary known oblique angle at different x- coordinates along the B-scan plane.

[0074] This shear wave excitation strategy avoids the bias caused by unexpected waves traveling out of plane by only exciting waves in known directions. The elastography derivation may differ compared to the reverberant elastography derivation, which assumes that waves are traveling in all directions. In the case of in-plane, waves only travel within the plane of the B-scan, which may be described by the constraint 4> = 0. Given a general <p(x) (including <p = c), trigonometric correction may be applied to correct the in-plane B-scan measurements for out-of-plane oblique wave propagation. For example, for a known constant angle cpknown, the measured in-plane wavelength may be corrected by dividing it by cos((pknown). Assuming, for ease of explanation, <p = 0, n_q is now limited to the xy-plane. Recomputing the ensemble average integral and autocorrelation, results in:

The autocorrelation is still a sum of spherical Bessel functions. However, the first term is scaled by half, which results in the central lobe of the Bessel function being narrower. FIG. 5 shows the original and new Bessel function sums on the same axes. A narrower central lobe may reduce autocorrelation window size requirements and moderately improve resolution over a fully reverberant example.

[0075] In some examples, known wave propagation direction may be determined based on placement of the imaging system with respect to the excitation system. For instance, in FIG. 1A, contacts 103, 108 may generate plane waves propagating along the x-axis (at least in the region containing the B-scan plane 104. Thus, aligning the B-scan plane with this axis results in waves traveling in-plane with the B-scan plane 104. Rotating imaging system 106 around the z axis by <p results in waves traveling obliquely to the scan plane 104 and the angle of rotation cp may be used to correct for this. Inclusion of additional contacts, vibration of the contacts at various compound frequencies, arbitrary scan surfaces, etc. may introduce other known directional components, which may be corrected for using standard trigonometric techniques.

[0076] One method to reduce wavelength is increase wave frequency; however higher frequency waves also attenuate more rapidly, making it harder to produce a diffuse wave field. Therefore, for each clinical target, it is critical to understand the underlying tissue architecture and expected range of shear wave speeds in order to design the most effective imaging scheme.

[0077] The coupling of the displacement-generating source is not limited to an individual region. The contact geometry is tunable, enabling one or more points of coupling, which produces a more complex displacement field compared to previous methods. The energy source may also be driven by an arbitrary signal, which further contributes to the complexity of the displacement field. The energy source is also not required to be in direct contact with the sample of interest but may instead be coupled though one or more nearby structures in series or in parallel. While a complex, random-like displacement field is more challenging to measure by conventional tools, it provides the advantage of being more robust to motion, while also being less sensitive to the specific parameters of energy-coupling, such as focal distance and alignment between the energy source and sample.

[0078] FIG. 6 shows an example of a system 600 for elastography in accordance with some embodiments described in the present disclosure. As shown in FIG. 6, a computing device 650 can receive one or more types of data (e.g., imaging data, including raw data and/or preprocessed data) from data source 602. For example, computing device 650 may comprise a control unit or other computer described herein. Data source 602 may comprise an imaging system, such as an OCT imaging system (e.g., a Doppler OCT system, etc.), a high speed camera, an ultrasound imaging system, an MRI system, etc. In some such examples, computing device 650/data source 602 may further include wave field inducement actuators, signal generators, etc. In some embodiments, computing device 650 can execute at least a portion of an elastography process 640 to generate elastography data (e.g., displacement data for a biological sample) from data received from the data source 602.

[0079] Additionally or alternatively, in some embodiments, the computing device 650 can communicate information about data received from the data source 602 to a server 652 over a communication network 654, which can execute at least a portion of the elastography process 640. In such embodiments, the server 652 can return information to the computing device 650 (and/or any other suitable computing device) indicative of an output of the elastography process 640.

[0080] In some embodiments, communication network 654 can be any suitable communication network or combination of communication networks. For example, communication network 654 can include a Wi-Fi network (which can include one or more wireless routers, one or more switches, etc.), a peer-to-peer network (e.g., a Bluetooth network), a cellular network (e.g., a 3G network, a 4G network, etc., complying with any suitable standard, such as CDMA, GSM, LTE, LTE Advanced, WiMAX, etc.), other types of wireless network, a wired network, and so on. In some embodiments, communication network 654 can be a local area network, a wide area network, a public network (e.g., the Internet), a private or semi-private network (e.g., a corporate or university intranet), any other suitable type of network, or any suitable combination of networks. Communications links shown in FIG. 6 can each be any suitable communications link or combination of communications links, such as wired links, fiber optic links, Wi-Fi links, Bluetooth links, cellular links, and so on.

[0081] Referring now to FIG. 7, an example of hardware 700 that can be used to implement data source 602, computing device 650, and server 652 in accordance with some embodiments of the systems and methods described in the present disclosure is shown.

[0082] As shown in FIG. 7, in some embodiments, computing device 650 can include a processor 702, a display 704, one or more inputs 706, one or more communication systems 708, and/or memory 710. In some embodiments, processor 702 can be any suitable hardware processor or combination of processors, such as a central processing unit (CPU), a graphics processing unit (GPU), and so on. In some embodiments, display 704 can include any suitable display devices, such as a liquid crystal display (LCD) screen, a light-emitting diode (LED) display, an organic LED (OLED) display, an electrophoretic display (e.g., an "e-ink" display), a touchscreen, and so on. In some embodiments, inputs 706 can include any suitable input devices and/or sensors that can be used to receive user input, such as a keyboard, a touchscreen, a microphone, and so on.

[0083] In some embodiments, communications systems 708 can include any suitable hardware, firmware, and/or software for communicating information over a data bus, communication network 654 and/or any other suitable communication networks. For example, communications systems 708 can include one or more transceivers, one or more communication chips and/or chip sets, and so on. In a more particular example, communications systems 708 can include hardware, firmware, and/or software that can be used to establish a Wi-Fi connection, a Bluetooth connection, a cellular connection, an Ethernet connection, and so on. In another example, such as where data source 702 and computing device 750 are cohoused, communication systems 708 may include an internal interconnect, such as USB interconnect, general purpose input/out (GPIO) interconnect, etc. [0084] In some embodiments, memory 710 can include any suitable storage device or devices that can be used to store instructions, values, data, or the like, that can be used, for example, by processor 702 to present content using display 704, to communicate with server 652 via communications system(s) 708, and so on. Memory 710 can include any suitable volatile memory, non-volatile memory, storage, or any suitable combination thereof. For example, memory 710 can include random-access memory (RAM), read-only memory (ROM), electrically programmable ROM (EPROM), electrically erasable ROM (EEPROM), other forms of volatile memory, other forms of non-volatile memory, one or more forms of semi-volatile memory, one or more flash drives, one or more hard disks, one or more solid state drives, one or more optical drives, and so on. In some embodiments, memory 710 can have encoded thereon, or otherwise stored therein, a computer program for controlling operation of computing device 650. In such embodiments, processor 702 can execute at least a portion of the computer program to present content (e.g., images, user interfaces, graphics, tables), receive content from server 652, transmit information to server 652, and so on. For example, the processor 702 and the memory 710 can be configured to perform the methods described herein.

[0085] In some embodiments, server 652 can include a processor 712, a display 714, one or more inputs 716, one or more communications systems 718, and/or memory 720. In some embodiments, processor 712 can be any suitable hardware processor or combination of processors, such as a CPU, a GPU, and so on. In some embodiments, display 714 can include any suitable display devices, such as an LCD screen, LED display, OLED display, electrophoretic display, a computer monitor, a touchscreen, a television, and so on. In some embodiments, inputs 716 can include any suitable input devices and/or sensors that can be used to receive user input, such as a keyboard, a mouse, a touchscreen, a microphone, and so on.

[0086] In some embodiments, communications systems 718 can include any suitable hardware, firmware, and/or software for communicating information over communication network 654 and/or any other suitable communication networks. For example, communications systems 718 can include one or more transceivers, one or more communication chips and/or chip sets, and so on. In a more particular example, communications systems 718 can include hardware, firmware, and/or software that can be used to establish a Wi-Fi connection, a Bluetooth connection, a cellular connection, an Ethernet connection, and so on.

[0087] In some embodiments, memory 720 can include any suitable storage device or devices that can be used to store instructions, values, data, or the like, that can be used, for example, by processor 712 to present content using display 714, to communicate with one or more computing devices 650, and so on. Memory 720 can include any suitable volatile memory, non-volatile memory, storage, or any suitable combination thereof. For example, memory 720 can include RAM, ROM, EPROM, EEPROM, other types of volatile memory, other types of non-volatile memory, one or more types of semi-volatile memory, one or more flash drives, one or more hard disks, one or more solid state drives, one or more optical drives, and so on. In some embodiments, memory 720 can have encoded thereon a server program for controlling operation of server 652. In such embodiments, processor 712 can execute at least a portion of the server program to transmit information and/or content (e.g., data, images, a user interface) to one or more computing devices 650, receive information and/or content from one or more computing devices 650, receive instructions from one or more devices (e.g., a personal computer, a laptop computer, a tablet computer, a smartphone), and so on.

[0088] In some embodiments, the server 652 is configured to perform the methods described in the present disclosure. For example, the processor 712 and memory 720 can be configured to perform the methods described herein (e.g., the method of FIG. 4, the method of FIG. 5).

[0089] In some embodiments, data source 602 can include one or more actuators 722, one or more data acquisition systems 724 (e.g., imaging devices), one or more communications systems 726, and/or one or more signal generators 728. In some embodiments, actuators 722 can be any suitable device to induce a wave field in a biological sample, such as a mechanical vibration source (e.g., a piezoelectric actuator, eccentrically weighted motor, MEMS vibrator, etc.) an air source (e.g., a nozzle and air supply, etc.), sound source (e.g., ultrasonic transducer, sonic transducer, etc.), light source (e.g., laser, LED, quantum dot, etc.), etc. In some embodiments, the one or more data acquisition systems 724 comprise imaging device to acquire data from a biological sample. For example, data acquisition systems 724 may comprise galvanometers, scan lenses, interferometers, OCT laser sources, etc. In some embodiments, one or more portions of the data acquisition system(s) 724 can be removable and/or replaceable. In some embodiments, the one or more signal generators 728 can include any suitable devices to generate signal(s) for the one or more actuators 722, such as amplifiers, digital-to-analog converters (DACs), function generator circuitry, etc. As a particular example, computing device 650 may transmit a control signal via communications system 726 to cause signal generator(s) 728 to generate one or more driving signals for actuators 722. For example, as described above, driving signals may include pure tones, complex tones, multiple independent driving signals, etc. [0090] In some embodiments, communications systems 726 can include any suitable hardware, firmware, and/or software for communicating information to computing device 650 (and, in some embodiments, over communication network 654 and/or any other suitable communication networks). For example, communications systems 726 can include one or more transceivers, one or more communication chips and/or chip sets, and so on. In a more particular example, communications systems 726 can include hardware, firmware, and/or software that can be used to establish a wired connection using any suitable port and/or communication standard (e.g., USB, RS-232, ISC, 12C, GPIO, etc.), Wi-Fi connection, a Bluetooth connection, a cellular connection, an Ethernet connection, and so on.

[0091] In some embodiments, any suitable computer-readable media can be used for storing instructions for performing the functions and/or processes described herein. For example, in some embodiments, computer-readable media can be transitory or non- transitory. For example, non-transitory computer-readable media can include media such as magnetic media (e.g., hard disks, floppy disks), optical media (e.g., compact discs, digital video discs, Blu-ray discs), semiconductor media (e.g., RAM, flash memory, EPROM, EEPROM), any suitable media that is not fleeting or devoid of any semblance of permanence during transmission, and/or any suitable tangible media. As another example, transitory computer- readable media can include signals on networks, in wires, conductors, optical fibers, circuits, or any suitable media that is fleeting and devoid of any semblance of permanence during transmission, and/or any suitable intangible media.

[0092] FIGS. 8-17 illustrate various aspects of experiments performed by the inventors and which illustrate various aspects of the technology described herein. These experiments are illustrative of various embodiments described above. Any aspects of the technology described with respect to these experiments may be applied in embodiments described above. Additionally, any suitable alternatives to particular aspects of these experiments may be applied in embodiments.

[0093] FIG. 8 illustrates experimental aspects of experiments performed to validate coherent recovery of the shear wave field while asynchronously imaging. An agar phantom with 1.5 w/v% agar and 0.5 w/v% intralipid was imaged synchronously with 2 Khz piezoelectric excitation (Digi-key 445-181631-ND) and MB-mode scanning to obtain ground truth measurements of the shear wave field. The piezoelectric actuator was in contact with the phantom through 8 independent 3D-printed prongs evenly spaced at a diameter of 20 mm. The data were acquired with a custom OCT system, the OCT source was an optically k-clocked, phase- stable, 100 kHz VCSEL wavelength-swept laser (Thorlabs, SL131090) operating with a center wavelength of 1310 nm. The beam was raster-scanned by galvanometers (Thorlabs, GVS002) through a scan lens (Thorlabs, LSM04) with field of view up to 14.1 mm in both lateral dimensions. The galvanometers were aligned to avoid introducing a phase shift during scanning. The lateral resolution of the system was 10 microns and the axial resolution was 6 microns. The complex fringe was acquired with polarization diverse receivers (Advanced Fiber Resources, China) to improve the SNR of displacement measurements. The incoming signal was then digitized with an acquisition card (AlazarTech, Canada; ATS 9373, 4 GS/s) with input for k-clocking.

[0094] Displacements were computed from the Doppler phase difference relative to the initial A-line at each location. Surface wave correction was applied to compensate for sample surface motion and refractive index difference between the sample and air.

[0095] A synchronized B-scan was obtained with MB-mode scanning (FIG. 8B) by capturing 800 A-lines in a single location before moving onto the next. This was done at 1,024 locations across 5mm. Since the repetition rate of the laser is 100 Khz and the shear waves are 2 Khz, each period of the shear wave field spanned exactly 50 A-lines. From the fully coherent data, down-sampling was performed to emulate the effect of raster-scanning by only keeping one A-line perglobal time step. Then, the experimental demodulation pipeline was applied to fully recover to the original coherent shear wave field. This signal processing pipeline is discussed below.

[0096] In this experiment, a conventional OCT setup with a phase-stable laser (VCSEL) was employed to acquire OCE images (FIG. 8A). The galvanometer (Thorlabs, GVS002) can scan the OCT beam in either synchronous MB-mode or asynchronous BM-mode (FIG.7B) through a scan lens (Thorlabs, LSM04). Example 3D printed plows are shown in (FIG. 8C) with a tissuemimicking agar phantom. The 2 layer agar phantom was composed of a lower stiff layer and upper soft layer (FIG. 8D). For finger imaging, the thumb was placed between the plows with B-scans running across the finger pad (FIG. 8E). After acquiring data before hydration, lotion was applied, and, after 5 minutes, a second dataset on the same finger was acquired.

[0097] With respect to mechanical contrast in a tissue-mimicking phantom, a 2-layer agar phantom was produced with contrasting stiffness to demonstrate asynchronous semi- reverberant elastography in a heterogeneous sample. The lower stiff layer was composed of 1.5 w/v% agar, while the upper soft layer was composed of 0.75 w/v% agar. Both layers had 0.5 w/v% intralipid so that the OCT structural images would look identical despite each layer having different mechanical properties. The data were acquired with the same rasterscanning, k-clocked, phase-stable, 100 kHz VCSEL wavelength-swept laser source (Thorlabs) operating with a center wavelength of 1310 nm. However, instead of synchronizing with the piezoelectric wave source, the data were acquired conventionally by BM-mode rasterscanning. To acquire a volume, 2 consecutive B-scans of 1,792 A-lines each were acquired consecutively in the same y-location before moving to the next out-of-plane y-location. The number of A-lines was chosen to ensure Nyquist sampling across a 14 mm field of view. In total, 100 y- locations over 14 mm were sampled. Given the laser repetition rate of 100 Khz and there being 1,792 A-lines per B-scan, 1981 Hz was chosen for the shearwave excitation frequency to maximize the displacement between B-scans. Shown in FIG. 3C, shear waves were excited in the sample by vertical plows that only generate shear waves in-plane with each B-scan. Shear waves were driven by a function generator through a high performance voltage amplifier (MMech. MX200) and piezoelectric actuator (Digi-key 445-181631-ND).

[0098] With respect to elasticity in vivo before and after hydration, semi-reverberant elastography was used to image skin in vivo on human skin from a finger before and after hydration. First, the skin was imaged prior to hydration with asynchronous semi- reverberant excitation. The plows were positioned on either side of the finger, running parallel to the length of the finger. B-scans were acquired across the finger, orthogonal to the plows. Then, the skin was lathered with thick lotion (Lubriderm, Advanced Therapy) for 5 minutes. After 5 minutes, the same region of skin was imaged (FIG. 3E). Whereas prior methods would have taken up to 6 minutes to acquire this data, the procedure was able produce the elastography results from only 3.58 seconds of imaging. This is a greater than 100-fold improvement in acquisition time, while also being a reasonable amount of time for a human subject to remain still.

[0099] Described with respect to FIG. 9, the following steps describe in detail the signal processing pipeline used in MATLAB to measure shear modulus from an OCT volume. First, the tissue surface was detected (FIG. 9A) and each column was shifted accordingly in order to flatten it (FIG.. 9B). Motivation for tissue flattening primarily arises from the fact that the autocorrelation is computed across the entire B-scan and flattening supports selecting windows where the tissue is relatively homogeneous. For example, flattening allows separation the epidermis from the dermis more effectively, making the assumption of homogeneity more realistic. Additionally, dominant shear waves in strongly layered tissues such the skin are likely to be either Rayleigh waves or lamb waves, both of which travel horizontally along layers.

[0100] After flattening the tissue, the Doppler-OCT phase difference is computed between A-lines at each location by multiplying one frame of the complex tomogram with the conjugate of another frame at a later time step. The phase difference is proportional to displacement by where S_z is the axial displacement, is the central wavelength of light 1300 nm, n is the refractive index of tissue which wasapproximated as 1.4, and ep is the Doppler phase difference. To average out noise, Gaussian filter was applied along z with a kernel size of 4 pixels. Surface wave correction (FIG. 9C) was also applied to compensate for sample surface motion and refractive index difference between the sample and air.

[0101] With a phase-stable laser, there is negligible phase jitter or phase noise, which simplifies processing. Two B-scans are used to determine displacement at each out of plane y-location, so motion artifacts are also relatively minimal. In a volume with two B-scans at 100 different y-locations, there will be 200 total B-scans. After computing the Doppler phase difference, this is reduced to 100 displacement frames. 2D unweighted phase-unwrapping was performed (FIG. 9C) along z to get the unwrapped displacement field. These displacements are still amplitude modulated by the scanning cosine term. To demodulate the displacement field, the ID fast Fourier transform was performed along x to get F_X{S_Z } (FIG. 9E). Here, three dominant peaks are visible: two peaks at the positive and negative carrier frequency in addition to one peak around zero for the DC component.

[0102] Following the transform, the entire spectrum was downshifted by the carrier frequency such that the desired signal is centered at zero. Afterdownshifting, a rectangular window to the signal to retrieve the complex- valued displacement field with a reasonable spatial bandwidth. We define D{F_X{S_Z }; f , fw }, where D represents downshifting, f is the amount by which we want to downshift (equal to the carrierfrequency),and f_w is the width ofthe rectangularwindow. Given shearwave excitation frequency wo and raster-scanning speed v_Scan, f = Wo/vscan. Finally, we compute the inverse Fouriertransform alongxto recoverthefullycoherentcomplex-valuedshearwavedisplacement field, S = F ^D FxfSz }; f , fw }} (FIG. 9F).

[0103] From the complex-valued displacement field S, the autocorrelation along x was computed (FIG. 1G). Note that prior works considered only the real part of the displacement, thus discarding half of the data. By considering the complex-valued field, information was captured from the entire harmonic period of the shear wave field from just one frame. The autocorrelation window spans the entire length of x as well as some rows along z for better averaging. Here, 10 rows were used, which is equivalent to 60 microns. Then, the best fit shear wave number k was found from the sum of square errors relative to Eq. 4. To dynamically adjust for tissue regions with different Bessel function widths, the second zero crossing was found and the fit length was adjusted to only use the first two lobes. Since the autocorrelation spans the entire width of a B-scan, this leaves a ID set of k values at each y- location. From the shear wave number, the shear wave speed c was computed based on the shear w a ve frequency by c = wo/k. For a sample with constant density p, shear modulus G is directly related to shear wave speed by c = VG/p. All together, from an OCT volume with 1,000 pixels in depth, 1, 792 A-lines per B-scan, and 200 B-scans, a 2D shear wave number map was generated with 100 pixels in depth and 100 pixels acrossy.

[0104] To demonstrate recovery of the coherent shear wave field while raster-scanning, an agar phantom was imaged with synchronization for ground truth, then down-sampled the data to emulate scanning (FIG. 10). For the fully synchronized data, 800 A-lines at 1,024 x- locations were acquired. Given the laser repetition rate of 100 Khz and shear wave excitation of 2 Khz, each shear wave period consisted of 50 A-lines. 16 whole periods of the shear wave field were collected to ensure sufficient data in case the synchronization drifted. However, only 50 A-lines were technically required. The raster-scanned data on the other hand only requires 1 A-line per location after down-sampling, resulting in at least a 50 fold theoretically improvement in acquisition time (varying based on laser and shear wave frequency). Finally, the signal was amplitude demodulated to remove the modulation that was introduced by the raster-scanning time delay between A-lines to verify that the result was the same as the synchronized ground truth data. Demodulation was performed it straight-forwardly in the frequency domain and successfully demonstrated recovery of the fully coherent wave field without any need for slow synchronization (FIG. 10C). The modulation effect also conveniently separates out the desired signal from low frequency noise by shifting the displacement field to the carrier frequency. Therefore, when a bandpass filter is applied to recover the displacement field, the low frequency noise is also attenuated naturally as part of the signal processing steps.

[0105] FIG. 11 illustrates elastography results in a tissue-mimicking phantom with mechanical contrast. From the OCT intensity image, there is no major difference by eye between top and bottom layer besides their boundaries. The sample was excited with in-plane 1,981 Hz shear waves through 2 plows on either side of the field of view. Through scanning demodulation, the spatially-coherent displacement field along x can be recovered, from which the autocorrelation can be computed. Applying the rest of the signal processing pipeline, the shear wave number, k and well as shear modulus G may be measured. Note that the data were collected entirely with conventional raster-scanning, and the methods described herein were used to successfully recover the coherent displacement field (FIG. 11B). This enabled us the use of in-plane semi-reverberant Bessel function sums that were derived to fit shear wave number. As expected, since shear wavelength is larger in stiffer tissue, the central lobe of the autocorrelation in the stiff region is wider (FIG. 11D). Since the phantom is homogeneous along x, and the autocorrelation may be computed along the entire x dimension, collapsing each xz-plane of data to a single dimension of k values allows presentation of elastography results along the yz-plane.

[0106] FIGS. 12A-F illustrate elasticity in vivo before and after hydration. Asynchronous semi- reverberant elastography drastically reduces acquisition time and therefore enables practical in vivo imaging of elasticity. The data were collected on skin from a human finger pad before and after hydration with lotion. OCT intensity images of the x-averaged yz- plane show that there is no difference in the structure of the tissue from before or after hydration. The epidermis and dermis layers are visible, supporting the fact that the skin has a layered architecture enabling averaging across x. For the processing, the tissue was flattened finding the surface and shifting the tomogram accordingly. This also allowed averaging across the x direction without more significant mixing of different regions of tissue. The unflattened images are shown in FIG. 12 for visualization.

[0107] Volumes were acquired with 200 B-scans, each with 1,792 A-lines, so the total acquisition time was 3.58 seconds. During acquisition, 1,981 Hz was applied for in-plane shear wave excitation with plows. The data were raster-scanned, with 2 consecutive B-scans at each of 100 y-locations. The 14x14 mm field of view en face OCT intensity image of the surface of the finger is shown in FIG. 12E, from which the finger print can be observed. From each 3D volume there was measured a 2D map of shear wave numbers k. Shear wave numbers are related to shear modulus by tissue density, so assuming a constant tissue density, quantitative shear modulus can also be measured. Finally, the depth dependent shear modulus for dry versus hydrated skin was computed. This is the shear modulus at each depth, averaged across the entire layer. As expected, dry skin is stiffer so it has a higher shear modulus compared to hydrated skin.

[0108] FIG. 13 illustrates an experiment demonstrating shear waves generated from a single contact point.

[0109] FIG. 14 illustrates an example validation of shear wave field measurements with full 3D wave recovery (e.g., across x and y directions). The OCE system's laser was an optically k- clocked, phase stable, 100 kHz VCSEL wavelength-swept laser (SL131090, Thorlabs, USA) operating with a center wavelength of 1310 nm. The collimated input beam had a diameter of 3.6 mm; combined with a 54 mm effective focal length objective scan lens (LSM04, Thorlabs, USA), the system produced an e'² diameter lateral resolution volume of 25 microns. The beam was raster-scanned by galvanometers (GVS002, Thorlabs, USA) through the scan lens with a field of view up to 14.1 mm. The complex fringe was acquired with polarization diverse receivers (Advanced Fiber Resources, China). The incoming signal was digitized with a digitizer with k-clocking (ATS 9373, 4 GS/s, AlazarTech, Canada). For the rabbit corneas, shear waves were coupled into the sample via direct contact with a 3D-printed, 8-pronged excitation ring with diameter of 13 mm.

[0110] Mature New Zealand White rabbit excised eyes (Pel-Freez Biologicals) were used within 30 minutes post-mortem. They were secured in a custom-designed holder, and intraocular pressure (IOP) was maintained using a water column and monitored with a rebound tonometer (IC200, iCare, USA). The corneal epithelium was removed manually, and baseline central corneal thickness (CCT) was measured using ultrasonic pachymetry (Pachmate 2) to serve as an indirect marker of corneal hydration. Saline drops were used to prevent dehydration. The experiment followed a repeating sequence: IOP and elastography measurements were performed, followed by increasing the water column to a new, increased

IOP. [0111] FIG. 14A illustrates an lOP-controlled rabbit cornea OCE experimental setup. FIG. 14B is a photo of an eye globe holder, IOP water column needle, scan lens, and 3D printed prongs. FIG. 14C illustrates OCT reflectance B-scan showing the structure of the cornea and the ocular anatomy below. B-scans are depth resolved images formed on the xz-plane. FIG. 14D illustrates the depth-averaged en face OCT reflectance volume shows the field of view, top of the 3D printed prongs, and the pupil. FIG. 14E illustrates synchronous and asynchronous methods measure the same Doppler OCT shear wave field. Without y- correction, shear waves are distorted along y. Without %-correction, shear waves are also distorted along x.

[0112] Synchronized, phase-locked, MB-mode OCE measurements of lOP-controlled excised rabbit eyes were acquired (FIGS. 14A-D) for comparison. Phase-locked scanning was performed by capturing many consecutive A-lines in a single lateral location before moving the beam. This was done at 100 locations across a 15 mm-wide B-scan. This was repeated at 100 sequential y-locations to produce an entire 3D volume. Since the repetition rate of the laser is 100 kHz and shear waves were 2 kHz, each period of the shear wave field spanned exactly 50 A-lines. However, 100 A-lines were acquired at each location in order to account for scanner repositioning time between lateral sampling locations. The same 3D field of view was also acquired with asynchronous raster-scanning, performing just 2 B-scans of 1,536 A- lines each across lOO y-locations. In total, the synchronized dataset had (z, x, y, t) dimensions of (2048, 100, 100, 100), while asynchronous imaging was only (2048, 1536, 100, 2), representing a 3-fold reduction in acquisition time despite a 15-fold increase in lateral sampling along the x-axis. For the same lateral sampling along x, this would have been a 50- fold reduction in acquisition time, although increased sampling may enable advanced motion correction techniques for imaging in human subjects.

[0113] Mature New Zealand White rabbit cadaver eyes were imaged within 30 minutes post-mortem. IOP was monitored closely with a rebound tonometer (IC200, iCare, USA) and controlled with a water column. Shear waves at 9 frequencies were measured across three lOPs: 8.7, 16, and 30 mmHg. At each IOP, the phase gradient of the autocorrelation of the complex displacement field was used to measure shear wave number k to create an experimental dispersion curve. These curves were used to fit a theoretical dispersion curve from the NITI model (FIG. 15). From the fitted curves, shear moduli of 4.3, 5.5, and 8.6 kPa were measured as IOP increased. [0114] FIGS. 15A-D illustrate phase gradient shear wave number measurements and development of an NITI model. FIG. 15A shows 3D, coherent shear wave fields averaged across the entire depth of the cornea for a range of frequencies in a rabbit cornea ex vivo. Shear wave speckle size decreases as frequency increases. The displacement field is masked by the outline of the cornea. Corresponding phase gradient maps for each displacement field are also shown below. The phase gradient is computed from the autocorrelation of the complex displacement field, which already contains information from the displacement field at all time points. The phase gradients are further masked by the amplitude of the displacement, which is low in areas of speckle nulls. FIG. 15B illustrates from the phase gradient shear wave number k measurements, use of Pitre et al.'s NITI model to fit shear modulus G across lOPs of 8.7, 16, and 30 mmHg as measured by rebound tonometer (IC200, iCare, USA). The effective shear modulus of the cornea increases as IOP increases due to being in a different regime of the nonlinear stress-strain curve of the cornea.

[0115] Asynchronous shear wave imaging recovers the entire coherent 3D wave field. However, to recapitulate quantitative measures of elasticity, one must still extract shear wave number k from the displacement field and use a constitutive model to derive elasticity. The determination of shear wave number k from the phase gradient of the 2D shear wave displacement field was performed because it is impartial to shear wave conditions— phase gradient measurements work in single wave, partia lly-diffuse, and ful ly-diffuse shear wave field conditions. This increases flexibility for methods of shear wave generation in the cornea. [0116] Given a harmonic shear wave within a small window, it was modeled as a plane wave by e^l(kxX+ky^y> . Given this representation, the autocorrelation function of the wave at a distance Ax and Ay is equivalent to multiplication with a shifted complex conjugate,_{or e}-i(fcx(^+Ax)+fc_y(y+Ay)) The argument, or phase, of the autocorrelation function directly results in /c_xAx + k_y y. From this, k = ^k^ + k^ can be computed by taking the magnitude of gradient of the phase of the autocorrelation function with respect to determined values of Ax and Ay, which we set at Ax = Ay = 0. Phase gradient results from a rabbit cornea ex vivo across 9 frequencies are shown in FIG. 15A. As frequency increased, measured shear wave number also increased and shear wave speckle size decreased. At lower frequencies, the shear wavelength was larger than the field of view in some places. The nature of the phase gradient method still allowed recovery shear wave number. [0117] Given the complex-valued and coherent shear wave displacement field, the argument, or phase, of the autocorrelation function directly result in c_xA% + k_y y. From this, /c = was computed by taking the magnitude of gradient of the phase of the autocorrelation function with respect to determined values of A% and Ay, which we set at Ax = Ay = 0. For both the rabbit data and subject data (see FIG. 16), pre-processed displacement tomograms had dimensions of (z, x, y) = (1400, 1408, 72). Prior to computing the autocorrelation function, they down sampled the x dimension to 72 to match the y dimension. The autocorrelation function was then computed with square windows of 6 by 6 pixels. From each window, we computed the argument, or phase, and its gradient magnitude to determine k.

[0118] By repeating this operation at multiple excitation frequencies, they built an experimental dispersion curve /c(ro₀) and used it to determine mechanical properties of the cornea from an appropriate model. Corneal mechanical behavior is dominated by the stroma, composed of vertically stacked lamellae with collagen fibrils, organized with relatively random orientations in a hydrogel matrix, resulting in weak in-plane and strong out-of-plane anisotropies. Advanced models may be used for the cornea, including composite models and fully anisotropic non-linear models. In the experiment, a transverse isotropic model was used to serve as a proof of concept for 3D asynchronous shear wave imaging in humans in vivo. In particular, a constitutive model that includes near incompressibility and assumes negligible longitudinal anisotropy was used. The nearly incompressible transverse isotropic (NITI) model has shown promising results simulating the mechanical response of human and porcine corneas ex vivo and in rabbit corneas in vivo. Following the derivation and using shared MATLAB functions of Pitre et al., guided-mode solutions were determined along with their respective dispersion curves for a planar cornea on top of a fluid water layer. In the NITI model, two shear moduli are distinguished, namely /J. and G, governing in-plane and out-of- plane shear deformations. However, it has been shown that in porcine corneas, the antisymmetric mode solution Ao is mostly insensitive to /z. Consequently, since Ao is the primary mode observed in the experiments, p was set as p = 36G, leaving G as the sole parameter to fit. Other material parameters were assumed constant, and corneal thickness was determined via OCT reflectance images. Here, wave numbers are measured at discrete shear wave frequencies, so the Ao mode dispersion curve was fitted with experimentally determined wave numbers for a finite set of frequencies as opposed to the entire 2D dispersion power spectrum. FIG. 15C shows the experimental points measured at 9 frequencies for the rabbit cornea alongside their respective dispersion curved fits. The fitting procedure yielding the best-fit G values and the 95% confidence intervals is described below. [0119] FIG. 16 illustrates an experiment of corneal imaging with indirect wave field generation. FIG. 16A illustrates a subject operating an actuator as a handheld accessory only in contact with the eyelid. FIG. 16B illustrates OCT reflectance B-scan showing the cornea and surrounding ocular anatomy. FIG. 16C illustrates the depth-averaged en-face OCT reflectance volume showing the field of view, cornea, sclera, eyelids, and eyelashes. The B-scan comes from the location of the dashed line. FIG. 16D illustrates power spectra of the displacement field at each frequency. These spectra have already been frequency-shifted so the amplitude modulation carrier frequency is centered at zero. The brighter center band shows the original DC noise. There are also two copies of the desired displacement signal at the positive and negative carrier frequency. FIG. 16E illustrates motion between scans in human subjects in vivo, showing the desirability of robust motion correction. Bulk motion correction along z and x across all 4 frequencies shows movement up to 30 pm between pairs of B-scans at the same y-location. In the out-of-plane y-dimension, there is movement up to 250 pm. FIG. 16 F illustrates axial motion. Sufficient lateral Nyquist sampling within a B-scan may be used for non-rigid axial motion correction along z. FIG. 15G illustrates demodulated, coherent, en- face, depth averaged, motion corrected displacement fields at 4 frequencies from a human subject in vivo.

[0120] For corneal imaging, the piezo was handheld by each subject via a 3D-printed mount which contacted the eyelid through a single blunt point with a 2 mm radius. Displacements were computed from the Doppler phase differences relative to the initial A-line at each location. They also applied surface wave correction to compensate for sample surface motion and refractive index difference between the sample and air. Shear wave excitation was performed by piezoelectric actuators (445-181631-ND, Digi-key, USA) via a high-performance voltage amplifier (MX200, Micromechatronics Inc., USA).

[0121] A piezoelectric actuator (445-181631-ND, Digi-key, USA) was electrically isolated with electrical tape and mounted in a 3D printed case. The case itself had a single blunt prong on the tip that contacted the eyelid, as depicted in FIG. 16. Since the eyelid was in direct contact with the eye and the piezo was in direct contact with the 3D printed case, shear waves propagated into the cornea. Subjects held the device in their own hand for positioning. For human subject imaging, the piezoelectric device generated shear waves at 944, 1335, 1660, and 2115 Hz. Due to asynchronous imaging, the shear wave displacement fields were amplitude modulated to a carrier frequency. As an example of the carrier frequency shift, the mean power spectrum as a function of depth for each frequency is shown in FIG. 16D for one dataset in one eye from one subject. The mean power spectrum shows sufficient displacement signal at each frequency.

[0122] Phase-locked, synchronous shear wave imaging methods tend to use sparse sampling of the OCT tomogram, typically acquiring just 100 A-lines over many millimeters. This is because shear waves are relatively large (>1 mm), so for completely static tissue, oversampling does not necessarily provide additional information about the wave field. On the other hand, to measure accurate Doppler OCT phase differences in vivo, the same speckle is compared between B-scans. At the employed lateral resolution, motion artifacts of just ~50 pm can completely degrade displacement measurements. In conventional MB-mode imaging, even Nyquist sampling of the OCT speckle laterally, which would increase sampling and imaging time by one order of magnitude (from 10 seconds to up to minutes), does not guarantee that motion will be mitigated. ID M-mode scans lack sufficient tissue contrast to enable cross-correlation approaches to determine motion magnitude and identify the correct M-mode scan for a given physical location in the sample.

[0123] Applying the described technology for shear wave imaging, an entire 3D volume of laterally Nyquist sampled OCT frames can be acquired in just 2.3 seconds over a 13.5 by 13.5 mm lateral field of view. Furthermore, motion-induced speckle decorrelation occurs at the time scale of the two B-scan repetitions that are acquired at each y-location: only 30 ms. For this reason, the fundamental sensitivity to motion is at least an order of magnitude lower than in conventional MB-mode scanning. Moreover, the speckle decorrelation arising from in-plane motion can be corrected by use of Nyquist sampling along x, which supports motion correction techniques that have already seen extensive development for OCTA. Motion detection and correction results from one dataset from one subject are shown as an example in FIGS. 16F and F.

[0124] The motion correction pipeline was a bulk registration computed between each pair of B-scans at individual y-locations. This bulk correction determined one bulk x-shift and z- shift per pair of B-scans. Subpixel image registration was computed by a cross-correlation which was upsampled by a factor of 40 using zero-padding. After bulk motion correction, they also performed a non-rigid axial motion correction between pairs of B-scans with batches of 64 A-lines. Aligning the speckle supports computing the Doppler phase differences these B- scans. In addition, to address macroscopic deformation of the shear wave field due to motion, zx bulk motion was performed correction between B-scan pairs acquired at each y-location. Again, the cross-correlation for subpixel registration was upsampled by a factor of 40. However, since the cornea naturally has a curved surface, straightforward motion correction would tend to flatten desired topology features. To account for the curvature of the cornea, the motion for each B-scan pair was averaged across all excitation frequencies and subtracted from the correction.

[0125] After tomogram reconstruction and motion correction, they detected the surface of the cornea and flattened the tissue accordingly. Flattening the tissue enabled convenient implementation of the used autocorrelation approach since it uses entire x-dimension during demodulation. Flattening ensured each window consists entirely of corneal tissue. Then, they computed the Doppler-measured displacements between the pair of B-scans acquired at each y-location and unwrap the 2D Doppler phase map. Surface wave correction was then performed to compensate for tissue surface motion and refractive index differences between the tissue and air.

[0126] To demodulate the effect of raster-scanning within a B-scan, expected demodulation frequency was calculated. For shear wave excitation frequency ru₀ and raster-scanning speed v_scan, the amplitude modulation term was cos ( o₀ x/v_scan). Then, the data was corrected for small changes in the demodulation pixel due to lateral motion between B-scans in a given pair. Shifting the amplitude modulated spectrum back to DC in the frequency domain is equivalent to applying a phase slope in the spatial domain. This also allowed use of noninteger demodulation pixel values while ensuring that the desired signal is centered on DC. The Fourier transform along x for each z-depth was then computed, and a Blackman bandpass filter centered on DC was applied to remove noise and the mirror copy of the signal at the other carrier frequency. Following, the inverse Fourier transform was performed resulting in a complex-valued displacement field. For any desired phase shift <p, the complexvalued field by e¹^ was multiplied and the real component computed. Now, the shear fields were spatially and temporally coherent along the xz-plane, but out-of-plane y-locations were still misaligned. [0127] To correct for this, the expected time delay (and corresponding phase shift depending on the excitation frequency) was calculated between consecutive y-locations. For a time-delay of At_y between y-locations (based on the laser repetition rate and sampling parameters), the expected phase shift is_e^(lAt-^{y in}-^{0 )}. Cumulatively multiplying this phase shift to each xz-displacement field corrected for artifactual phase jumps of the shear wave field between each y-location (see FIG. 14E). Since the galvos did not scan perfectly orthogonal to each other, waves tended to propagate spuriously along y due to an additional phase delay. This was calibrated with the following procedure: first, rabbit cornea data with the exact same imaging parameters as the human subjects was acquired. The fast and slow axes were switched and MB-mode data was acquired for ground truth. Following, it was determined which additional phase shift between y-locations resulted in identically reconstructed shear wave fields for the 3 datasets. For the system, the calibration was found to be an additional delay of 2.32 A-lines per y-location. After this, the time delay introduced by lateral motion detected between B-scan pairs was computed: lateral motion causes a given cross-sectional displacement to be measured earlier or later than expected by the deterministic time delay and galvanometer calibration. Some spurious phase jumps alongy remained in the recovered shear wave field, presumably due to undetectable out-of-plane lateral motion. To correct for this, the average phase difference between y locations was calculated, followed by cumulative correction along the entire y-axis. The resulting vector unwrapped and fitted with a first-order polynomial. This was subtracted from the cumulative correction, which was applied as the correction. This was fitted to and subtracted from the first-order polynomial of the phase correction in order to avoid subtracting a meaningful phase ramp along y that would be produced by the projection on the y axis of a dominant wave in a single-wave regime.

[0128] FIG. 17 illustrates a study of asynchronous shear wave imaging that was performed on ten healthy volunteers over the age of 18 with no known history of eye disease.: FIG. 17A illustrates experimental dispersion curves (points) and dispersion curves that were fit to the data (lines) across all three subjects are shown side by side. FIG. 17B: shear modulus G measurements for each eye from each acquisition are also shown side-by-side (points). They also combined shear wave number data from all four acquisitions per eye to fit a single shear modulus. FIG. 17C: shear wave SNR across all four excitation frequencies for each dataset shows that indirect shear wave generation through the eyelid is sufficient for generating OCT- measurable displacements.

[0129] The study was conducted at Massachusetts General Hospital with approval from the Institutional Review Board under Protocol 2024P0002050. Each subject went through the informed consent process and was trained on how to operate the handheld piezoelectric device on their eyelid. An effective method for device placement was for the subject to first close their eye, then touch the blunt tip of the case to the outer corner of their eyelid, and finally slowly open their eye. This ensured consistent contact with the soft part of the eye, avoiding the browbone, while also ensuring that subjects would not accidentally touch their sclera or cornea.

[0130] Four measurements were performed on each eye, for a total of eight datasets per subject. Each measurement consisted of four consecutive volumes, during which the piezoelectric device was driven by a different frequency during each volume. Across all ten of the subjects that were recruited and enrolled in the study, we present here data from three subjects. The first four datasets were excluded because imaging parameters (such as shear wave frequency and reference arm polarization state) were improved for latter subjects. In addition, a computer data loss event affected the datasets of an additional subject, leaving just five subjects. Of these five subjects, one did not finish the study. A second subject was not able to use the piezoelectric device to generate shear waves in their cornea, and they could not detect enough shear wave signal for analysis.

[0131] Experimental dispersion curves and the resulting fitted curve for each of the three remaining subjects are shown in FIG. 17A. Shear modulus G across the three subjects ranged from 5 to 50 kPa and was relatively consistent when measured from the same eye. To characterize the quality of indirect-contact shear wave excitation through the eyelid, they also calculated displacement signal-to-noise ratio (SNR) as the ratio between the magnitude of displacement measured in the cornea versus air. Shear wave SNR as a function of frequency across all measurements is shown in FIG. 17C. For the most part, it was consistent per subject and remained above 10 dB. Subject 2 had lower shear wave SNR at higher frequencies, but sufficient to measure consistent shear moduli across 8 datasets. Across all subjects, shear wave SNR decreases as frequency increases. This is expected due to increased attenuation of shear waves at higher frequencies. They were able to repeatedly generate shear waves in humans, even with varying pressure, varying positions on the eyelid, and among different subjects.

[0132] Applying the described technology, the entire temporally and spatially coherent 3D shear wave field can be recovered using a conventional A-line-rate OCT system without the need for phase-locked imaging. This approach reduces sensitivity to motion from minutes or seconds down to milliseconds, supporting phase-sensitive displacement measurements in tissue which rely on repeated measurements of OCT speckle. Combined with indirect-contact shear wave excitation via piezoelectric device, they demonstrated OCE measurements from three healthy human volunteers in vivo enabling ophthalmic elastography.

[0133] One limitation of the NITI model is that it neglects the role of IOP, and the fitting procedure yields an apparent shear modulus G. In fact, it is well known that the cornea stiffens with increasing IOP as collagen fibers uncrimp and begin to stretch. This non-linear behavior is confirmed in a rabbit eye ex vivo (FIG. 15D) wherein the measured corneal elasticity from the same cornea increases with IOP. This could explain differences in shear modulus measurements between the left and right eyes from the same subject who reported no history of eye disease. Since the piezoelectric device is manually pressed against the eye, albeit with protection from the eyelid, it is likely to pre-stress the anterior-segment and could interfere with shear modulus measurements. Another possibility is that these subjects simply had different IOP in each eye, or different corneal elasticity.

[0134] In the foregoing specification, implementations of the disclosure have been described with reference to specific example implementations thereof. It will be evident that various modifications may be made thereto without departing from the broader spirit and scope of implementations of the disclosure as set forth in the following claims. The specification and drawings are, accordingly, to be regarded in an illustrative sense rather than a restrictive sense.

[0135] As used herein in the context of computer implementation, unless otherwise specified or limited, the terms "component," "system," "module," "framework," and the like are intended to encompass part or all of computer-related systems that include hardware, software, a combination of hardware and software, or software in execution. For example, a component may be, but is not limited to being, a processor device, a process being executed (or executable) by a processor device, an object, an executable, a thread of execution, a computer program, or a computer. By way of illustration, both an application running on a computer and the computer can be a component. One or more components (or system, module, and so on) may reside within a process or thread of execution, may be localized on one computer, may be distributed between two or more computers or other processor devices, or may be included within another component (or system, module, and so on).

[0136] In some implementations, devices or systems disclosed herein can be utilized or installed using methods embodying aspects of the disclosure. Correspondingly, description herein of particular features, capabilities, or intended purposes of a device or system is generally intended to inherently include disclosure of a method of using such features for the intended purposes, a method of implementing such capabilities, and a method of installing disclosed (or otherwise known) components to support these purposes or capabilities. Similarly, unless otherwise indicated or limited, discussion herein of any method of manufacturing or using a particular device or system, including installing the device or system, is intended to inherently include disclosure, as embodiments of the disclosure, of the utilized features and implemented capabilities of such device or system.

Claims

CLAIMS What is claimed is:

1. A method, comprising: inducing a wave field in a biological sample, the wave field comprising waves having a wave frequency; scanning the sample at a scanning velocity to perform a plurality of spatial measurements of the sample; determining spatial data from the plurality of spatial measurements; amplitude demodulating the spatial data based on the wave frequency and the scanning velocity; and generating displacement data for the biological sample based on the amplitude demodulated spatial data.

2. The method of claim 1, wherein: scanning the sample comprises capturing a plurality of brightness scans (B-scans) at a scan plane location; determining the spatial data comprises determining a displacement signal for the scan plane location based on the plurality of B-scans; and amplitude demodulating the spatial data comprises amplitude demodulating the displacement signal.

3. The method of claim 2, further comprising: determining a plurality of complex-valued displacement signals for a corresponding plurality of scan plane locations based on a corresponding plurality of demodulated displacement signals for the plurality of scan plane locations; phase aligning the plurality of complex-valued displacement signals based on a time offset between B-line scans of successive scan plane locations.

4. The method of claim 2, wherein the wave field comprises waves having a known propagation direction with respect to the scan plane.

5. The method of claim 4, wherein the wave field comprises first plane waves having a first propagation direction parallel to the scan plane location and second plane waves having a second propagation direction parallel to the scan plane location.

6. The method of claim 4, wherein the wave field comprises waves traveling in a plurality of known propagation directions with respect to the scan plane.

7. The method of claim 1, wherein inducing the wave field comprises inducing a superposition of waves having different propagation directions in the biological sample.

8. The method of claim 7 , wherein the biological sample comprises a first body part and inducing the superposition of waves comprises inducing shear waves in a second body part.

9. The method of claim 7, wherein inducing the superposition of waves comprises inducing shear waves at a plurality of contact points on the biological sample.

10. The method of claim 1, wherein: the wave field further comprises waves having a second wave frequency; and amplitude demodulating the spatial data comprises: generating first amplitude demodulated data based on the first wave frequency and the scanning velocity; and generating second amplitude demodulated data based on the second wave frequency and the scanning velocity.

11. The method of claim 10, wherein: the wave field further comprises waves having a third wave frequency and a fourth wave frequency; and amplitude demodulating the spatial data comprises: generating third amplitude demodulated data based on the third wave frequency and the scanning velocity; and generating fourth amplitude demodulated data based on the fourth wave frequency and the scanning velocity.

12. The method of claim 1, further comprising determining the wave frequency based on the biological sample.

13. The method of claim 1, further comprising selecting a wavefront geometry for the wave field based on the biological sample.

14. A system, comprising: a driver to induce a wave field in a biological sample, the wave field comprising waves having a wave frequency; an imaging system to scan the sample at a scanning velocity to perform a plurality of spatial measurements of the sample; a processor; and a non-transitory computer readable medium storing instructions executable by the processor to: determine spatial data from the plurality of spatial measurements; amplitude demodulate the spatial data based on the wave frequency and the scanning velocity; and generate displacement data for the biological sample based on the amplitude demodulated spatial data.

15. The system of claim 14, wherein: the imaging system is to scan the sample by capturing a plurality of brightness scans (B-scans) at a scan plane location; the instructions to determine the spatial data comprise instructions to determine a displacement signal for the scan plane location based on the plurality of B-scans; and the instructions to amplitude demodulate the spatial data comprise the instructions to amplitude demodulate the displacement signal.

16. The system of claim 15, wherein the non-transitory computer readable medium stores further instructions executable by a processor to: determine a plurality of complex-valued displacement signals for a corresponding plurality of scan plane locations based on a corresponding plurality of demodulated displacement signals for the plurality of scan plane locations; phase align the plurality of complex-valued displacement signals based on a time offset between B-line scans of successive scan plane locations.

17. The system of claim 15, wherein the wave field comprises waves having a known propagation direction with respect to the scan plane.

18. The system of claim 17, wherein the wave field comprises first plane waves having a first propagation direction parallel to the scan plane location and second plane waves having a second propagation direction parallel to the scan plane location.

19. The system of claim 17, wherein the wave field comprises waves traveling in a plurality of known propagation directions with respect to the scan plane.

20. The system of claim 14, wherein the driver is to induce the wave field by inducing a superposition of waves having different propagation directions in the biological sample.

21. The system of claim 20, wherein the biological sample comprises a first body part and the driver is to induce the superposition of waves by inducing shear waves in a second body part.

22. The system of claim 20, wherein the driver is to induce the superposition of waves by inducing shear waves at a plurality of contact points on the biological sample.

23. The system of claim 14, wherein: the wave field further comprises waves having a second wave frequency; and the instructions to amplitude demodulate the spatial data comprise instructions to: generate first amplitude demodulated data based on the first wave frequency and the scanning velocity; and generate second amplitude demodulated data based on the second wave frequency and the scanning velocity.

24. The system of claim 23, wherein: the wave field further comprises waves having a third wave frequency and a fourth wave frequency; and the instructions to amplitude demodulate the spatial data comprise: instructions to generate third amplitude demodulated data based on the third wave frequency and the scanning velocity; and instructions to generate fourth amplitude demodulated data based on the fourth wave frequency and the scanning velocity.

25. The system of claim 14, wherein the instructions are further executable to determine the wave frequency based on the biological sample.

26. A non-transitory computer readable medium storing instructions executable by a computer to: control a driver to induce a wave field in a biological sample, the wave field comprising waves having a wave frequency; receive sample scan data comprising a plurality of spatial measurements of the sample from an imaging system having a scanning velocity; determine spatial data from the plurality of spatial measurements; amplitude demodulate the spatial data based on the wave frequency and the scanning velocity; and generate displacement data for the biological sample based on the amplitude demodulated spatial data.

27. The non-transitory computer readable medium of claim 26, wherein: the sample scan data comprises a plurality of brightness scans (B-scans) at a scan plane location; the instructions to determine the spatial data comprise instructions to determine a displacement signal for the scan plane location based on the plurality of B-scans; and the instructions to amplitude demodulate the spatial data comprise the instructions to amplitude demodulate the displacement signal.

28. The non-transitory computer readable medium of claim 27, storing further instructions executable by a processor to: determine a plurality of complex-valued displacement signals for a corresponding plurality of scan plane locations based on a corresponding plurality of demodulated displacement signals for the plurality of scan plane locations; phase align the plurality of complex-valued displacement signals based on a time offset between B-line scans of successive scan plane locations.

29. The non-transitory computer readable medium of claim 27, wherein the wave field comprises waves having a known propagation direction with respect to the scan plane.

30. The non-transitory computer readable medium of claim 29, wherein the wave field comprises first plane waves having a first propagation direction parallel to the scan plane location and second plane waves having a second propagation direction parallel to the scan plane location.

31. The non-transitory computer readable medium of claim 29, wherein the wave field comprises waves traveling in a plurality of known propagation directions with respect to the scan plane.

32. The non-transitory computer readable medium of claim 26, wherein the instructions to control the driver to induce the wave field comprise instructions to control the driver to induce a superposition of waves having different propagation directions in the biological sample.

33. The non-transitory computer readable medium of claim 32, wherein the biological sample comprises a first body part and the driver is to induce the superposition of waves by inducing shear waves in a second body part.

34. The non-transitory computer readable medium of claim 32, wherein the driver is to induce the superposition of waves by inducing shear waves at a plurality of contact points on the biological sample.

35. The non-transitory computer readable medium of claim 26, wherein: the wave field further comprises waves having a second wave frequency; and the instructions to amplitude demodulate the spatial data comprise instructions to: generate first amplitude demodulated data based on the first wave frequency and the scanning velocity; and generate second amplitude demodulated data based on the second wave frequency and the scanning velocity.

36. The non-transitory computer readable medium of claim 35, wherein: the wave field further comprises waves having a third wave frequency and a fourth wave frequency; and the instructions to amplitude demodulate the spatial data comprise: instructions to generate third amplitude demodulated data based on the third wave frequency and the scanning velocity; and instructions to generate fourth amplitude demodulated data based on the fourth wave frequency and the scanning velocity.

37. The non-transitory computer readable medium of claim 26, wherein the instructions are further executable to determine the wave frequency based on the biological sample.

-M-