System and method for wave velocity determination

Abstract

Systems and methods for measuring material properties of a medium are described. Operations include receiving data from a receiver such as an ultrasound receiver, the data including a representation of particle motion over time in a spatial dimension, generating a time delay map from the data, wherein a first axis of the time delay map corresponds to each of a plurality of points in the spatial dimension as a measurement point and wherein a second axis of the time delay map corresponds to each of the plurality of points in the first spatial dimension as a reference point; determining a time delay surface based on a maximum value of a cross-correlation function for the time delay map, wherein the cross-correlation function is configured to determine a relationship between a plurality of measurement points and a plurality of reference points, and determining a wave velocity based on the time delay surface.

Description

Docket No.2023-293 (630666.01478) SYSTEM AND METHOD FOR WAVE VELOCITY DETERMINATION STATEMENT OF FEDERALLY SPONSORED RESEARCH

[0001] This invention was made with government support under HL145268 awarded by the National Institutes of Health. The government has certain rights in the invention. BACKGROUND

[0002] This disclosure relates to the measurement of wave velocity for evaluation of mechanical properties. More particularly, this disclosure relates to a method that uses correlation of all available signal pair combinations to obtain robust measurements of wave velocity. In some examples, this method could be applied for vascular applications as well as shear wave elastography.

[0003] Ultrasonic shear wave elastography (SWE) has been applied in clinical practice to obtain the absolute value of stiffness of human tissues. Shear wave speed (SWS) maps or a single value of SWS are provided. Pulse Wave Imaging (PWI) can estimate local vascular stiffness by mapping the motion of the vessel wall caused by the natural pulsation of vessels. These methods benefit from an accurate measurement of the propagation speed of waves. The vessel wall motion can be complex because of, for example, turbulent hemodynamics. In comparative systems and methods, the decorrelation between signals can degrade the reliability of the cross-correlation. SUMMARY OF THE DISCLOSURE

[0004] The present disclosure addresses the aforementioned drawbacks by providing systems and methods which utilize the full information of spatiotemporal maps and cross- correlation.

[0005] The foregoing and other aspects and advantages of the present disclosure will appear from the following description. In the description, reference is made to the accompanying drawings that form a part hereof, and in which there is shown by way of illustration one or more embodiments. These embodiments do not necessarily represent the full scope of the invention, however, and reference is therefore made to the claims and herein for interpreting the scope of the invention. 1 QB\630666.01478\90616201.2Docket No.2023-293 (630666.01478) BRIEF DESCRIPTION OF THE DRAWINGS

[0006] FIG.1 illustrates B-mode of in vivo common carotid artery. Top wall and bottom walls are illustrated.

[0007] FIG.2 illustrates example spatio-temporal maps of wall motion. (a) illustrates the velocity of a wall (b) illustrates the acceleration of the wall.

[0008] FIG. 3 illustrates an example decorrelation between a reference signal and a studied signal.

[0009] FIG. 4 illustrates an example time delay map between the reference lines and the subsequent lines.

[0010] FIG. 5A illustrates an example of a map of the maximum value of the normalized cross-correlation coefficient between the reference lines and the subsequent lines.

[0011] FIG. 5B illustrates an example of the weight matrix using correlation coefficients. (a) illustrates the maximum value of the normalized cross-correlation, (b) illustrates pixels of (a) in which values are less than a threshold (i.e., 0.6) are set to be zero, and (c) illustrates binary weight after thresholding.

[0012] FIGS. 6A-6B illustrate example surfaces of time delays and linear regression planes.

[0013] FIG. 7 illustrates examples of (a) two-dimensional full correlation, and (b) multi-dimensional full correlation.

[0014] FIG.8 illustrates in vivo data of six patients.14 measurements were made at the outflow vein, 5 cm away from the arteriovenous fistula (AVF) for each subject. (a) illustrates Median pulse wave velocity (PWV) and (b) illustrates IQR (interquartile range) of PWV.

[0015] FIG. 9 illustrates an example SWS map for a homogeneous phantom (E = 45 kPa). (a) illustrates the tangent plane method and (b) illustrates multi-dimensional full cross- correlation (FCor).

[0016] FIG.10 illustrates an example workflow for estimating shear wave speed (SWS) according to embodiments described in the present disclosure.

[0017] FIG.11 illustrates an example of an ultrasound system.

[0018] FIG. 12 is a flowchart of an example method for estimating wave velocity according to embodiments described in the present disclosure. 2 QB\630666.01478\90616201.2Docket No.2023-293 (630666.01478)

[0019] FIG.13 illustrates an example workflow for estimating wave velocity according to embodiments described in the present disclosure. DETAILED DESCRIPTION

[0020] As noted above, ultrasonic shear wave elastography (SWE) has been widely applied in clinical practice to obtain the absolute value of stiffness of human tissues. A SWS (shear wave speed) map or a single value of SWS are provided by this modality. Pulse Wave Imaging (PWI) can estimate local vascular stiffness by mapping the motion of the vessel wall caused by the natural pulsation of vessels. These methods rely on the accurate measurement of the propagation speed of waves.

[0021] The displacement data along the depth dimension can be obtained using a 1D autocorrelation method from in-phase / quadrature (IQ) data or using cross-correlation from radiofrequency data. In a comparative example, the interframe axial displacements were estimated from two consecutive images and this is often called the particle velocity data. The top and bottom vessel walls in a longitudinal view are manually traced (FIG. 1). In the figure, yellow lines represent the vessels walls. The walls are tracked across different frames.

[0022] The waves at the walls can be caused by ARF (acoustic radiation force) for SWE or the natural pulsation of vessels for PWI. Considering PWI as an example, the propagation of the pulse wave is visualized by plotting the particle velocity of the walls coded in color, over time, along the path of propagation in the longitudinal direction as shown in FIG. 2(a). These maps are called spatio-temporal maps. FIG.2(a) corresponds to the axial wall velocity temporal waveforms at each lateral position v(i), i = 1, 2, …, N (N = 128)). In the case of PWI, the temporal derivative of the velocity (i.e., wall acceleration) may be used as shown in FIG.2(b). In the illustrated example, the axial wall acceleration was derived by performing temporal differentiation, a(i) = dv(i) / dt.

[0023] The wave velocity (i.e., pulse wave velocity (PWV) for PWI or group / local velocity for SWE) can be calculated by time delay (Δt) between two lateral positions (their coordinates are xiand xj). The wave velocity, V is given by ^^ ൌ^^^ െ ^^^(1) Δ^^

[0024] The time delay can be estimated using cross-correlation. Cross-correlation is an operation that quantifies the similarity between two signals, f(t) and g(t), as a function of the delay of one relative to the other, τ. It is defined mathematically as follows: 3 QB\630666.01478\90616201.2Docket No.2023-293 (630666.01478) ^^ ^^ ^^ ^^^ ^^^ ൌ∑௧ ^^^ ^^ ^ ^^^ ^^^ ^^^ (2)^∑௧ ^^^ ^^ ^ ^^^ଶ ∑௧ ^^^ ^^^ଶ ^^⁄ ଶ

[0025] The thetime delay at which the are the similarity of f(t) and g(t).

[0026] In a comparative method, the cross-correlation functions between one reference line (a(refi), a blue line in FIG. 2(b) as an example) and all the subsequent lines of the spatio- temporal map were calculated to obtain the time delay between them. The time delays were plotted against the distance traveled by the pulse wave. A linear regression fit was applied on the time-distance plot. The wave velocity was calculated as the reciprocal of the linear regression. This will be referred to herein as SRef, short for single reference. In another comparative method, the reference lines were extended to N positions, calculating the cross- correlation function between the acceleration (or velocity) profile at the first position and the acceleration (or velocity) profiles at all the other N − 1 positions. This process was repeated using each of the N lines of the spatio-temporal map successively as the reference signal. N plots of the time delay versus distance were created. From each plot, the wave velocity was calculated, and the N wave velocity estimates are averaged to yield the final estimate value.

[0027] The wall motion can be complex because of the turbulent hemodynamics in the lumen of the vessel. Thus, in the comparative methods, the decorrelation (FIG.3) between the two signals (a(refi) and a(j)) degrades the reliability of the cross-correlation. The present disclosure describes systems and methods, referred to herein as full cross-correlation (FCor), to utilize the full information of the spatiotemporal maps and cross-correlation.

[0028] By providing for systems and methods that are more robust to turbulent hemodynamics, that provide for more reliable cross-correlation, and that generate more accurate and reliable wave velocities, the present disclosure effects improvements in the technological fields ultrasound imaging, including sub-fields of ultrasonic shear wave elastography and pulse wave imaging.

[0029] Advantageously, the disclosed systems and methods implement a full correlation-based method for measuring pulse wave propagation, which utilizes all lateral positions as reference signals. A time delay map with all combinations of reference signals and signals at all lateral positions is calculated and correlation coefficients are used as weights. A plane is fit to the time delay data using weighted least squares, for example, which provides pulse wave velocity values. The disclosed systems and methods provide an advantage over 4 QB\630666.01478\90616201.2Docket No.2023-293 (630666.01478) previous techniques for measuring pulse wave velocity because good reference signals can be selectively utilized based on correlation coefficients; thus, the disclosed systems and methods are more robust to noise at some lateral positions.

[0030] The time delays may be calculated between a(i) and a(j) (i, j =1, 2,…, 128) by finding the time lag of the maximum of the cross-correlation as shown in FIG.4. Note that all lateral positions (denoted as x-axis) contribute to wave velocity calculation as reference lines (denoted as xr-axis). The maximum value of the normalized cross-correlation function characterizes the similarity of a reference signal and a studied signal (FIG.5A). Combinations for which correlation coefficients are less than a threshold (i.e., 0.5 in this example) are set to zero in the figure, and they are excluded for wave velocity calculation. FIG.6A shows the time- delay surface plot of FIG. 4 after excluding points with their correlation coefficients are less than 0.5.

[0031] A tangent plane of this time-delay surface at a center position is given by ^^ ^^^ ^^ ^^ ൌ^^^ ^ ^ ^(3)^^ ^^ ^^ ^^^ ^^^^

[0032] where C is the

[0033] The data point Δt(i,j) on the time-delay surface can be modeled when considering the measurement follows:Δ^^^^^, ^^^ൌ Δ ^^̅^^^, ^^^^ ^^(4)

[0034] where Δt̅ (i,j) is theε is a measurement error term. Assuming the true time delay satisfies the equation of the tangent plane in a homogeneous media, the following expression was obtained: Δ^^̅^ ^^, ^^^ ൌ^^ ^^ ^^ ^^ ^1 ^^ (5)^^ ^^ ^^^ ^^^ ^^^ ^^^ ^ ^^ ൌ ^ ^^ ^^^ ^ ⋅ ^^

[0035] where xi and xirare the coordinates in x and xraxes and B is defined by ^^^ ் (6) ^^ ൌ ^^^ ^^^ ^^ ^^ ^^ ^^ ^^^൨

[0036] To find B, the error e between the observed and true time-delay in a least squares sense was minimized as follows: ேேேேଶ ^(7)Docket No.2023-293 (630666.01478)

[0037] where N is the number of lateral positions. In matrix form, this is given by: ^^ ൌ | ^^െ ^^ ^^|ଶ(8)

[0038] where |•| is the Euclidean norm, X is a N2× 2 matrix whose one row is [xi xjr], ^^ … ^^ … ்^^ ൌ ^^ ^ ^^^ … ^^^ … ^^^ … ^^^ … ^^ே … ^^ே … ^^ே (9)^^^ ^ ^ ^ ^ ^ ^ ^ ^^ … ^^^ … ^^ே … ^^^ … ^^^ … ^^ே … ^^^ … ^^^ … ^^ே ൨

[0039] and Eqn.(8) is a least-squares problem, and its closed form solution is given by ^^ ൌ ൫X்X൯ି^X்Y(10)

[0040] This full correlation method can be extended to increase robustness to outlier data in in vivo data, which may be caused by physiological motion, low displacement signal- to-noise ratio, and / or spatial inhomogeneities in tissue. In some implementations, the time delay of FIG. (4) may be weighted because different data points may not be equally reliable. A weighted least-squares method is the generalization of the ordinary least-square method that weights each measurement data unequally. The error e of Eqn. (8) may be modified to incorporate weights of value Wij as follows: ேே(1 ^^ ൌ ^^ ^^^^^ ^^^ ^^, ^^^ െ ^^^^^^^1) 1^^ ⋅ ^^^ଶ

[0041] In^ ଶ ฬଶ^(12) ^^ ൌ W ^^െ ^^ ^^^ฬ

[0042] where W is a diagonal matrix (N2× N2) whose element is Wij. The matrix B to minimize Eqn. (12) is given by ^^ ൌ ൫X்WX൯ି^X்WY(13)

[0043] The maximum value of the normalized cross-correlation function (i.e., correlation coefficient) of Eqn. (2) characterizes the similarity of the two studied signals. Combinations for which correlation coefficients are shown in FIG.5B(a), which can be utilized as a weight matrix. To ensure that some outliers (e.g., those whose correlation coefficients are relatively small) do not contribute to the TP estimation, elements in the weight matrix which are less than a threshold (i.e., 0.6 in this example) can be set to zero as shown in FIG. 5B(b). 6 QB\630666.01478\90616201.2Docket No.2023-293 (630666.01478) FIG. 5B(c) shows a binary weight after thresholding. FIG. 6B shows the time-delay surface plot of FIG.4 after excluding points whose correlation coefficients are less than 0.6.

[0044] In either case, the wave velocity is estimated as the inverse of the slope of the plane with respect to the x-axis, which is given by ^^ ൌ 1 / ฬ^^ ^^ (14)^^ ^^ฬ

[0045] Eqn. (11) was derived for the case where motion is two-dimensional (i.e., motion in x-t plane) as shown in FIG.7(a). This can be extended to process three-dimensional data as shown in FIG.7(b), where M positions are across the z-axis. When considering motion in the lateral (x) and axial directions (z), the wave velocity can be obtained solving the Eikonal equation. The arrival time satisfies the Eikonal equation if it is Lipschitz continuous: |∇ ^^^ ^^, ^^^| ൌ1 (15)^^

[0046] where T(x,z) is thethat is ∇^^ ൌ ൬^^ ^^ ^^ ^^ (16)^^ ^^ ,^^ ^^^

[0047] For the formulation ofbe considered as follows: ^x^,Δ ^^^^, ^xଶ,Δ ^^ଶ^, … , ൫x^,Δ ^^^൯, … , ^x^,Δ ^^^^(17)

[0048] whereEach input, xp is a column vector of length 4: x^ ൌ ^^^^^^^^ ^^^^^^^^்(18)

[0049] where x and xrand zrdenote depth positions and reference depth positions, respectively.

[0050] The Δtp(p = 1, …, P) points constitute a multi-dimensional surface of time- delays, and its tangent hyperplane at a center position is given by ^^ ൌ^^ ^^^^ ^^^ ^^ ^^^^^^ ^^^^ ^^^ ^^ ^^^^ ^^(19)

[0051] where C7 QB\630666.01478\90616201.2Docket No.2023-293 (630666.01478)

[0052] To fit the hyperplane to data points, represented as the time intercept b0 and slope parameters b1, b2, b3, and b4, the following linear relationship holds between the input and output data. ^^^^ ^^^^^^^ ^^^^^^ଶ^ ^^^^^ଷ^ ^^^^^^ସ^ Δ ^^^, ^^ ൌ 1, … , ^^ (20)

[0053] where ^^^^^^^^^ଶ^^ଷ^^ସ^ ൌ^ ^^^^ ^^^^ ^^ ^^ ^^ ^^ ^^ (21) ^^ ^^^^ ^^^^^ ^^^^ ^^^൨

[0054] The above can be written more compactly in a matrix form ^^^^^ ^ Δ ^^^, ^^ ൌ 1, … , ^^ (22)

[0055] where ^^^ ൌ ^1 ^^^^^^^ ^^^^^^^^ (23)

[0056] and^^ ൌ ^ ^^^^^^^^ଶ^^ଷ^^ସ^்(24)

[0057] The error e between the observed and hyperplane model in a least squares sense may be described as follows: ^ (2 ^൫ଶ5) ^^ ൌ ^^^ ^^ െ Y൯

[0058] where Y is a P × 1

[0059] The closed form solution of B to minimize e is given by ^^ ൌ൫X்X൯ି^X்Y (26)

[0060] where and X is a P × 5 vector whose pthelement is Xp. Eqns. (25) and (26) can also be modified to incorporate weights, as discussed above. From B and Eqn. (12), the wave velocity is given by: ^^ ൌ11 ൌ(27)

[0061] PWI andVerasonics Inc., Kirkland, WA, USA). To capture the wave propagation, PWC (plane wave compounding) was performed. The tracking beam had a center frequency of 5 MHz and its pulse repetition frequencies (PRF) for PWI and SWE were 1 kHz and 12.5 kHz, respectively. 8 QB\630666.01478\90616201.2Docket No.2023-293 (630666.01478) PWC excited all transducer elements with three angles [-4°, 0°, 4°] without transmit apodization. The receive beamformer (F / 1.5 with rectangular apodization) on the Verasonics system was used and beamformed IQ was processed to estimate motion. The frame rate of PWC was increased to the PRF using the time-aligned PWC (TAPWC) method on the acquired IQ data.

[0062] Most end-stage renal disease patients receive hemodialysis (HD) treatment, which requires gaining vascular access to a HD machine preferably with an arteriovenous fistula (AVF). Measuring the failure- to-mature of AVFs is critical to guide AVF interventional treatment. Six subjects have been enrolled in an in vivo study and 14 measurements were made at the outflow vein, 5 cm away from the AVF for each subject. IQR (interquartile range) was measured to compare results from FCor and SRef. Median PWV and IQR are shown in FIGS. 8(a,b). Diamonds and squares represent the top and bottom walls, respectively. Red and blue colors represent the FCorr and SRef results, respectively. FCorr provides comparable or smaller IQRs than SRef, indicating an increase in robustness. The use of the full data set can provide a more robust PWV measurement.

[0063] For a comparison study, an elasticity quality assurance phantom (Model 039) with nominal Young’s moduli of E = 45 kPa, manufactured by CIRS, Inc. (Norfolk, VA, USA) was scanned. One focused push beam (center frequency = 4.1 MHz, F / # = 1.5, focal depth = 20 mm) was transmitted for a duration of 400 μs.

[0064] The tangent plane method (TPM) estimates the wave velocity resorting to the Eikonal equation and calculates the tangent plane of the time-delay surface. This method can be considered to use a single reference when calculating the time-delays. A moving window was used for performing cross-correlations between the reference pixel and all other pixels within the window. The reference pixel was the center of the window, and the window size was 11 pixels.

[0065] FIG. 9(a,b) uses TPM and multi-dimensional full correlation, respectively. Within the red box of FIG.9(a), multi-dimensional full correlation shows less noise than TPM in the rectangle from z = 30-40 mm as well as near the top of the image z = 0-5 mm.

[0066] For another comparison study, a robust 2D shear wave speed calculation known as the tangent plane method (TPM) was adapted using the methods described in the present disclosure. In TPM, the inverse Eikonal equation is solved, and a tangent plane of the arrival time surface is computed to minimize a least squares error. Here, TPM is improved using the calculation of cross-correlation between all combinations of pixels in a window according to 9 QB\630666.01478\90616201.2Docket No.2023-293 (630666.01478) embodiments described in the present disclosure. Additionally, the normalized cross- correlation (NCC) coefficients are used as weights for the weighted least squares method when estimating the hyper-tangent plane.

[0067] The calculation of local shear wave speed in a window (11×11 pixels) with a center at (m,n) using a time-of-flight (TOF) estimation is shown in FIG. 10A. Motion was obtained using a 1D autocorrelation method from IQ data of consecutive frames. A directional filter to separate two directions of shear wave propagation ( ^^ = 0 or 180). TOF, ∆ ^^^ ^^, ^^^ (the travel time from the point (0,0) to ^ ^^^, ^^^^), was obtained by finding the lag, ^^, to maximize the cross-correlation coefficient between the point (0,0) to ^ ^^^, ^^^^. The tangent plane to fit this local TOF data points was estimated using a least square method. Considering a window, the data is three-dimensional with the size of 103×11×11 (103 is the number of frames of motion data). This data was converted to 2D data of the size of 103×121, which is shown in FIG.10B. As an illustration, motion at (m,n) and (i,j) pixel are plotted in FIG.10C. In previous studies, the cross-correlation functions between one reference point at (m,n) and all the neighboring pixels were calculated to obtain the time delays between them. Here, the disclosed systems and methods calculate TOF are between all combination of vel(k) and vel(g) (k, g = 1, 2, …, 121) as shown in FIG. 10D. Note that all pixels contribute to shear wave speed calculation as reference points. Hyper-tangent plane is estimated to fit TOF data points. For a more robust estimation, the NCC coefficients can be utilized as weights for a weighted least square method. Note that NCC characterizes the similarity of the two studied signals (FIG. 10E). Combinations for which NCC is less than a threshold are set to zero in FIG. 10F, and they can be excluded for shear wave speed calculation. Shear wave speed can then ne estimated as the inverse of the slope of the hyper-tangent plane.

[0068] The method could be expanded to a third spatial dimension and a temporal dimension to obtain the TOF ∆ ^^^ ^^, ^^, ^^^ (the travel time from the point (0,0) to ^ ^^^, ^^^, ^^^^). A hyperplane could be fit to the data for evaluating the wave velocity.

[0069] FIG.11 illustrates an example of an ultrasound system 100 that can implement the methods described in the present disclosure. The ultrasound system 100 includes a transducer array 102 that includes a plurality of separately driven transducer elements 104. The transducer array 102 can include any suitable ultrasound transducer array, including linear arrays, curved arrays, phased arrays, and so on. Similarly, the transducer array 102 can include a 1D transducer, a 1.5D transducer, a 1.75D transducer, a 2D transducer, a 3D transducer, and so on. 10 QB\630666.01478\90616201.2Docket No.2023-293 (630666.01478)

[0070] When energized by a transmitter 106, a given transducer element 104 produces a burst of ultrasonic energy. The ultrasonic energy reflected back to the transducer array 102 (e.g., an echo) from the object or subject under study is converted to an electrical signal (e.g., an echo signal) by each transducer element 104 and can be applied separately to a receiver 108 through a set of switches 110. The transmitter 106, receiver 108, and switches 110 are operated under the control of a controller 112, which may include one or more processors. As one example, the controller 112 can include a computer system.

[0071] The transmitter 106 can be programmed to transmit unfocused or focused ultrasound waves. In some configurations, the transmitter 106 can also be programmed to transmit diverged waves, spherical waves, cylindrical waves, plane waves, or combinations thereof. Furthermore, the transmitter 106 can be programmed to transmit spatially or temporally encoded pulses.

[0072] The receiver 108 can be programmed to implement a suitable detection sequence for the imaging task at hand. In some embodiments, the detection sequence can include one or more of line-by-line scanning, compounding plane wave imaging, synthetic aperture imaging, and compounding diverging beam imaging.

[0073] In some configurations, the transmitter 106 and the receiver 108 can be programmed to implement a high frame rate. For instance, a frame rate associated with an acquisition pulse repetition frequency (“PRF”) of at least 100 Hz can be implemented. In some configurations, the ultrasound system 100 can sample and store at least one hundred ensembles of echo signals in the temporal direction.

[0074] The controller 112 can be programmed to design an imaging sequence to capture images that may be used to implement the techniques described in the present disclosure, or as otherwise known in the art. In some embodiments, the controller 112 receives user inputs defining various factors used in the design of the imaging sequence.

[0075] A scan can be performed by setting the switches 110 to their transmit position, thereby directing the transmitter 106 to be turned on momentarily to energize transducer elements 104 during a single transmission event according to the designed imaging sequence. The switches 110 can then be set to their receive position and the subsequent echo signals produced by the transducer elements 104 in response to one or more detected echoes are measured and applied to the receiver 108. The separate echo signals from the transducer elements 104 can be combined in the receiver 108 to produce a single echo signal. 11 QB\630666.01478\90616201.2Docket No.2023-293 (630666.01478)

[0076] While FIG. 11 illustrates the ultrasound system 100 as including both a transmitter 106 and a receiver 108, in some embodiments the ultrasound system 100 may omit the transmitter 106 or the transmitter 106 may not be used or activated. In such implementations, rather than detecting waves applied from an external source, the receiver 108 may detect waves from endogenous sources (e.g., a heartbeat, a pressure pulse, etc.).

[0077] The echo signals are communicated to a processing unit 114, which may be implemented by a hardware processor and memory, to process echo signals or images generated from echo signals. As an example, the processing unit 114 can determine or estimate wave velocity using the methods described in the present disclosure. However, in some implementations the full cross-correlation may be performed by a controller or other processing unit that is external to the ultrasound system 100. Images produced from the echo signals by the processing unit 114 can be displayed on a display system 116. In implementations in which an external device is used (e.g., to perform the full cross-correlation), the images produced from the echo signals may be transmitted to the external device for display in addition to the processing.

[0078] The processing unit 114 may include a plurality of electronic processors each having one or more processing cores, and / or may include a single electronic processor having one or more processing cores. In implementations where multiple electronic processors are present, all of the electronic processors may form a part of the processing unit 114, all of the electronic processors may be part of one or more external devices, or the electronic processors may be distributed among the processing unit 114 and one or more external devices. The operations described herein may be performed by a single electronic processor or processing core, or may be divided among several electronic processors and / or processing cores in any combination.

[0079] FIG.12 illustrates an example method, and FIG.13 illustrates a corresponding example workflow, for estimating wave velocity in accordance with embodiments described in the present disclosure.

[0080] The method includes receiving data from a receiver, as indicated at step 1202. The receiver may be a transducer of an ultrasound system, for example. The receive is configured to acquire data about a wave field in a medium (e.g., a tissue including a vessel wall) over a period of time. The acquired data include a representation of particle motion over time in a first spatial dimension. When the medium is a vessel wall, the first spatial dimension may be a longitudinal dimension of the vessel wall. Additionally or alternatively, the data may also 12 QB\630666.01478\90616201.2Docket No.2023-293 (630666.01478) include a representation of particle motion over time in a second spatial dimension. When the medium is a vessel wall, the second spatial dimension may be a radial dimension of the vessel wall.

[0081] From the acquired data, a time delay map is generated, as indicated at step 1204. A first axis of the time delay map corresponds to each of a plurality of points in the first spatial dimension as a measurement point. A second axis of the time delay map corresponds to each of the plurality of points in the first spatial dimension as a reference point. When the data include a representation of particle motion over time in a second spatial dimension, an additional time delay map can be generated for each of a plurality of points in the second spatial dimension.

[0082] A time delay surface is then generated from the time delay map, as indicated at step 1206. In some implementations, the time delay surface may be a hyperplane. The time delay surface may be generated based on a maximum value of a cross-correlation function for the time delay map. The cross-correlation function can be configured to determine a relationship between a plurality of measurement points and a plurality of reference points. Determining the time delay surface can include discarding points on the time delay map having a value that is lower than a predetermined threshold. Additionally or alternatively, determining the time delay surface can include minimizing an error function. The error function may correspond to a difference between an observed time delay and a true time delay for each measurement point – reference point pair.

[0083] A wave velocity is then determined based on a derivative of the time delay surface, as indicated at step 1208.

[0084] From the wave velocity, one or more mechanical or material properties can be generated, as indicated at step 1210. The material properties may include Young’s modulus, shear modulus, bulk modulus, or the like. In some implementations, the material properties may include a local vascular stiffness of a vessel wall, which is determined based on the wave velocity.

[0085] The wave velocity and / or the material properties can then be output using a computer system, as indicated at step 1212. For instance, the wave velocity and / or material properties can be output as part of a report generated by the computer system. The report may include textual information, quantitative information (e.g., calculated values of wave velocity, material properties, etc.), and the like. Additionally or alternatively, the wave velocity can be output as a wave velocity map that indicates a spatial distribution of wave velocity values in 13 QB\630666.01478\90616201.2Docket No.2023-293 (630666.01478) the subject from which the data were acquired. Similarly, the material properties can be output as material property maps that indicate the spatial distribution of the material properties in the subject from which the data were acquired.

[0086] The present disclosure has described one or more preferred embodiments, and it should be appreciated that many equivalents, alternatives, variations, and modifications, aside from those expressly stated, are possible and within the scope of the invention. 14 QB\630666.01478\90616201.2

Claims

Docket No.2023-293 (630666.01478) CLAIMS What is claimed is:

1. A system for measuring material properties of a medium, the system comprising: a receiver configured to acquire data about a wave field in the medium over a period of time; and an electronic processor operatively connected to the receiver, the electronic processor configured to: receive the data from the receiver, the data including a representation of particle motion over time in a first spatial dimension, generate a time delay map from the data, wherein a first axis of the time delay map corresponds to each of a plurality of points in the first spatial dimension as a measurement point and wherein a second axis of the time delay map corresponds to each of the plurality of points in the first spatial dimension as a reference point, determine a time delay surface based on a maximum value of a cross-correlation function for the time delay map, wherein the cross-correlation function is configured to determine a relationship between a plurality of measurement points and a plurality of reference points, and determine a wave velocity based on a derivative of the time delay surface.

2. The system of claim 1, wherein to determine the time delay surface, the electronic processor is configured to discard points on the time delay map having a value that is lower than a predetermined threshold.

3. The system according to claim 1 or claim 2, wherein to determine the time delay surface, the electronic processor is configured to minimize an error function, the error function corresponding to a difference between an observed time delay and a true time delay for each measurement point – reference point pair.

4. The system of any one of claims 1 to 3, wherein the data further includes a representation of particle motion over time in a second spatial dimension, and 15 QB\630666.01478\90616201.2Docket No.2023-293 (630666.01478) the electronic processor is configured to generate the time delay map for each of a plurality of points in the second spatial dimension.

5. The system of claim 4, wherein the time delay surface is a hyperplane.

6. The system of any one of claims 1 to 5, wherein the medium is a tissue including a vessel wall.

7. The system of claim 6, wherein the electronic processor is further configured to determine a local vascular stiffness of the vessel wall based on the wave velocity.

8. The system of claim 6 or claim 7, wherein the first spatial dimension is a longitudinal dimension of the vessel wall.

9. The system of any one of claims 1 to 8, further comprising an ultrasound transmitter configured to produce the wave field in the medium, wherein the receiver is an ultrasound receiver.

10. A method of measuring material properties of a medium, the method comprising: receiving data from a receiver, the receiver being configured to acquire data about a wave field in the medium over a period of time, the data including a representation of particle motion over time in a first spatial dimension; generating a time delay map from the data, wherein a first axis of the time delay map corresponds to each of a plurality of points in the first spatial dimension as a measurement point and wherein a second axis of the time delay map corresponds to each of the plurality of points in the first spatial dimension as a reference point; determining a time delay surface based on a maximum value of a cross-correlation function for the time delay map, wherein the cross-correlation function is configured to determine a relationship between a plurality of measurement points and a plurality of reference points; and determining a wave velocity based on a derivative of the time delay surface. 16 QB\630666.01478\90616201.2Docket No.2023-293 (630666.01478) 11. The method of claim 10, wherein determining the time delay surface includes discarding points on the time delay map having a value that is lower than a predetermined threshold.

12. The method according to claim 10 or claim 11, determining the time delay surface includes minimizing an error function, wherein the error function corresponds to a difference between an observed time delay and a true time delay for each measurement point – reference point pair.

13. The method of any one of claims 10 to 12, wherein the data further includes a representation of particle motion over time in a second spatial dimension, and the method comprises generating the time delay map for each of a plurality of points in the second spatial dimension.

14. The method of claim 13, wherein the time delay surface is a hyperplane.

15. The method of any one of claims 10 to 14, wherein the medium is a tissue including a vessel wall.

16. The method of claim 15, further comprising determining a local vascular stiffness of the vessel wall based on the wave velocity.

17. The method of claim 15 or claim 16, wherein the first spatial dimension is a longitudinal dimension of the vessel wall.

18. The method of any one of claims 10 to 17, wherein the receiver is an ultrasound receiver.

19. The method of any one of claims 10 to 18, further comprising performing a temporal differentiation on the data to determine an acceleration of the medium. 17 QB\630666.01478\90616201.2Docket No.2023-293 (630666.01478) 20. A non-transitory computer-readable medium storing instructions that, when executed by at least one electronic processor of a system, cause the system to perform operations including the method of any one of claims 10 to 19. 18 QB\630666.01478\90616201.2

Images