METHOD AND APPARATUS TO MEASURE TISSUE DISPLACEMENT AND STRAIN

Abstract

A method for reconstructing displacement and strain maps of human tissue of a subject in vivo or other objects includes applying a mechanical deformation for a target area, imaging tissues while deformation is applied on the target area, measuring the axial and lateral displacements and axial, lateral, and shear strains of tissue in the target area, differentiating between tissues of different stiffness and strain responses, and labeling tissues of different stiffness and strain responses. In some embodiments, displacement and strain imaging are performed using a high-resolution, high-speed, highly-accuracy hierarchy recursive displacement tracking from conventional ultrasound brightness mode images. Particularly, this invention relates to the reconstruction of displacement and strain images from conventional ultrasound B-mode images acquired using any ultrasound machine of different manufacturers without the need to use carrier signals.

Claims

1. A method for reconstructing displacement and strain maps of objects under stimulus comprising: Applying a mechanical deformation for a target area of a body organ in a subject; imaging tissues while deformation is applied on the target area; measuring the axial and lateral displacements and axial, lateral, and shear strains of tissue in the targeted organ within human body in vivo due to the mechanical deformation using high-resolution, high-speed, highly-accuracy hierarchy recursive displacement tracking (HR-DT); differentiating between tissues of different stiffness and strain responses; and labeling tissues of different stiffness and strain responses.

2. The method of claims 1, wherein the variation of strain comprises: varying the strain by at least 0.5%.

3. The method of claim 1, wherein input images are acquired directly from ultrasound scanners using DICOM, Digital Video Interface (DVI), High-Definition Multimedia Interface (HDMI), video graphics adapter (VGA), S-Video, composite video out, BNC, or the universal serial bus (USB).

4. The method of claim 1, wherein input images are sent via DICOM, Ethernet, WI-FI, Bluetooth or other data exchange method.

5. The method of claim 1, wherein the mechanical deformation can be applied from outside or inside the body comprises deforming the target area.

6. The method of claim 1, wherein the deformation in tissue can be due to a thermal stimulus applied comprises changing the temperature of the target area.

7. The method of claim 1, wherein the mechanical deformation is applied by a transducer coupled to the ultrasound scanner.

8. The method of claim 1 can be applied to measure displacements and/or strains in objects using images from nuclear medicine cameras, optical cameras, x-ray based imaging, and magnetic resonance imaging.

9. A system for measuring the displacement and strain maps of objects under stimulus, the system comprising: a stimulating device to apply a mechanical deformation for a target area of a body organ in a subject; an ultrasound device to image tissues while deformation is applied on the target area using an ultrasound scanner; a device connected to ultrasound machine for measuring the axial and lateral displacements and axial, lateral, and shear strains of tissue in the targeted organ within human body in vivo due to the mechanical deformation using high-resolution, high-speed, highly-accuracy hierarchy recursive displacement tracking (HR-DT); a device to differentiate between tissues of different stiffness and strain responses; and a labelling device for tissues of different stiffness and strain responses.

10. The system in claim 9 uses conventional gray scale images of B-mode images, cine-loop and videos acquired directly from ultrasound machines via DICOM, Digital Video Interface (DVI), High-Definition Multimedia Interface (HDMI), video graphics adapter (VGA), S-Video, composite video out, BNC, or the universal serial bus (USB).

11. The system in claim 9 uses conventional gray scale images of B-mode images, cine-loop and videos sent via the DICOM, Ethernet, WI-FI, or Bluetooth.

12. The system in claim 9 uses conventional gray scale images of B-mode images, cine-loop and videos transferred via a data storage device.

13. The system of claim 9, wherein the displacement and strain change comprises a a mechanical deformation.

14. The system of claim 9, wherein the displacement and strain change comprises a thermal stimulus due to a temperature change.

15. The system of claim 9, is coupled to and integrated with the ultrasound imaging device.

16. The system of claim 9, is not coupled to nor integrated with the imaging device, wherein images are transferred remotely via network or data storage devices.

17. The system of claim 9, can be used to measure displacements and/or strains in objects using images from nuclear medicine cameras, optical cameras, x-ray based imaging, and magnetic resonance imaging.

18. The system of claim 9, can has his processing unit inside or outside the system like on an internal/external processor, cloud computing system, personal computer, laptop or smartphone or any other processing unit that can be supported with a display or already having one.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0011] A clear understanding of the key features of the invention summarized above may be had by reference to the appended drawings, which illustrate the method and system of the invention, although it will be understood that such drawings depict preferred embodiments of the invention and, therefore, are not to be considered as limiting its scope with regard to other embodiments which the invention is capable of contemplating. Accordingly:

[0012] FIG. 1 is an illustration of the system schematic of the ultrasound elastography add-on module system for providing non-invasive displacement and strain measurements.

[0013] FIG. 2 is an illustration of the first part of the displacement estimation algorithm starting from acquiring the ultrasound images to initial displacement values and determining skipping frames factor.

[0014] FIG. 3 is an illustration for the feature extraction and frame selection part in the algorithm.

[0015] FIG. 4 is an illustration for last part of the algorithm where elastogram is generated.

[0016] FIGS. 5A-5E illustrates (A) is a B-mode gray scale ultrasound image of breast phantom including 2 masses with a dashed line passing through the inclusion (red) and solid line on background material (blue), (B) axial displacement image estimated by HR-DT, (C) graph showing gradient of axial displacement along axial depth at the two lines (red and blue) shown (A), (D) axial strain by the HR-DT technique, and (E) axial strain by DSPview code [2] that used RF signals. Both (D) and (E) show the two masses in low strain values (blue). Ultrasound B-mode images were acquired using the Sonix RP machine by Ultrasonix (Canada).

[0017] FIGS. 6A-6D illustrates (A) B-mode gray scale ultrasound image of breast phantom including 2 masses B-mode, (B) axial displacement map, (C) axial strain map, (D) superimposed pseudo-strain image. Ultrasound B-mode images were acquired using Acuson machine by Siemens (Germany).

[0018] FIGS. 7A-7D illustrates (A) B-mode gray scale ultrasound image of breast phantom including 2 masses B-mode, (B) axial displacement map, (C) axial strain map, (D) superimposed pseudo-strain image. Ultrasound B-mode images were acquired using Digison machine by IBE (Egypt).

DETAILED DESCRIPTION OF THE INVENTION

[0019] Referring to FIG. 1 images are captured from the ultrasound machine (102) using a hardware interface (104) such as video capture device then they are processed by our method in the processing unit (106) where displacements and strains between different frames are measured and displayed on the display unit (108). Images can be transferred to the processing unit by means of network connections (110) through Ethernet or WI-FI and/or external data storage (112).

[0020] FIG. 2 refers to the first part of the displacement estimation algorithm where ultrasound images, cine-loop, or videos the region of interest (ROI) is automatically segmented at which there are imaging data (scan lines) acquired by the transducer. In the first block (202) the imaging depth is entered or determined automatically for processing successive or selected frames (204) acquired after time t. Template and search window size (206) are determined automatically by taking a ratio from the imaging frame size (such as 10-20% for 1st HR and 5-10% for the 2nd FIR) and or it can be adjusted manually by the user and also it's considered to not violate the minimum size of template based on speckle size and this is why we all user interfere. The template window is a vector of length equals to 2 or 3 elements and determines the template window size as well as the levels of the hierarchy used. Overlap percentage between template windows is determined to enhance the spectral resolution of the estimated displacement field between chosen imaging frames. There are two objective functions that can be used for the block matching. These functions are computationally efficient and/or accurate including first the Zero-mean normalized cross correlation (ZMCC)

[00001]$\begin{array}{cc}R\left(.Math..Math.x,.Math..Math.y\right)=\frac{{.Math.}_{x=\frac{\mathrm{tx}}{2}}^{\frac{\mathrm{tx}}{2}-1}.Math.{.Math.}_{y=\frac{-\mathrm{ty}}{2}}^{\frac{\mathrm{ty}}{2}-1}.Math.\left(f\left(x+.Math..Math.x,y+.Math..Math.y\right)-\stackrel{\_}{f}\right)\left(g\left(x,y\right)-\stackrel{\_}{g}\right)}{\sqrt[2]{\begin{array}{c}\left({.Math.}_{x=\frac{\mathrm{tx}}{2}}^{\frac{\mathrm{tx}}{2}-1}.Math.{.Math.}_{y=\frac{-\mathrm{ty}}{2}}^{\frac{\mathrm{ty}}{2}-1}.Math.f\left(x+.Math..Math.x,y+.Math..Math.y\right)-\stackrel{\_}{f}\right)\\ \left({.Math.}_{x=\frac{\mathrm{tx}}{2}}^{\frac{\mathrm{tx}}{2}-1}.Math.{.Math.}_{y=\frac{-\mathrm{ty}}{2}}^{\frac{\mathrm{ty}}{2}-1}.Math.\left(g\left(x,y\right)-\stackrel{\_}{g}\right)\right)\end{array}}}& \left(1\right)\end{array}$

[0021] where (x, y) is point of interest, f is the template area around the point of interest in the reference frame at t, g is the search area around the point of interest in the compressed frame at t+1, g is the mean of the search window, f is the mean of the template window, ty and tx are the axial and lateral template window size, respectively. The second function is the sum of absolute difference (SAD)

[00002]$\begin{array}{cc}R\left(.Math..Math.x,.Math..Math.y\right)={.Math.}_{x=\frac{\mathrm{tx}}{2}}^{\frac{\mathrm{tx}}{2}-1}.Math.{.Math.}_{y=\frac{-\mathrm{ty}}{2}}^{\frac{\mathrm{ty}}{2}-1}.Math..Math.f\left(x+.Math..Math.x,y+.Math..Math.y\right)-g\left(x,y\right).Math.& \left(2\right)\end{array}$

[0022] Enhanced Hierarchy displacement estimator

[0023] The frame selection described in block 204 is trivial such as frame t and t+1. However, the selection of frames used in elastography calculations depends on the frame rate used while acquiring the B-mode ultrasound images and the compression rate done by the user's hand. Frames are included in Elastography calculations whenever the strain between them is within a valid range to show a high quality Elastogram, i.e. if the average strain between two successive frames (t and t+1) is outside this range, frame t+1 is excluded and will be replace by t+2. If the strain between frame t and t+2 is between that range, frame t+2 will be included, and will be excluded otherwise, and so on. This technique assumes a maximum compression of 3% axial strain is applied with typically 1 Hz compression-relaxation cycles as example illustrated in [1], also skilled in the art can expand or shrink the range according to the requirements. The first run of the algorithm is reserved for determining the frame skipping rate and this is done by considering the average displacement and average strain between every 2 successive frames i.e. t & t+1 and then determine the skipping factor that is based on axial displacement and strain between frames to ensure us having a pure axial motion with efficient strain and this is what takes place in block 212.

[0024] The size and location of searching window relative to the template window have critical roles in determining the speed and accuracy of the algorithm and that is where we apply our constrains (axial, lateral, and hierarchical). At the 1.sup.st iteration of the hierarchical model (level 1) we choose search window size so that

sx=tx+offsetx

sy=ty+offsety (3)

where sx is the lateral dimension of the search window, sy is the axial dimension of the search window, tx is the lateral dimension of the template window, and ty is the axial dimension of the template window, offset x and offset y are chosen to form the search area. For example, these offsets can be 22 pixels in the axial direction and 12 pixels in the lateral direction assuming and with a displacement greater than 11 pixels in the axial direction and 6 pixels in the lateral dimension, we will get a bad quality strain image This has speed up our algorithm since at the 1.sup.st level of the hierarchy we usually choose a larger template window and it is frequent to choose the search window1.52 of the template window and this is computationally expensive specially for a large template window size.

[0025] Estimating the displacement field between two frames is performed using block matching techniques. As mentioned before, we can use mainly 2 equation (SAD & ZMCC). SAD is a fast and accurate function that we usually use for block matching so most of results are recorded using SAD function. In another embodiment, ZMCC can be used. In order to improve the processing speed and render a high quality axial strain image at almost real time processing, we adopt a smart algorithm to optimize the trade off between accuracy and speed for estimating high quality axial strain images. The 3 constrains (axial, lateral and hierarchical) are included in the hierarchy recursive (HR) speckle tracking. Axial and lateral constrains have been used before [2] but on estimating axial strain image using RF data not B-mode gray scale ultrasound images like this application. Here we are introducing a hybrid or merged technique involving the power of HR in estimating an accurate and high resolution displacement image along the high computational power added with the axial and lateral constrains based on our application. In one embodiment, we have observed that the 2nd levels of HR tracking lead to a high quality strain image at a reasonable time. At first level we may choose a large template window enough to divide our reference frame (pre-compression) into 1015 overlapping template windows at most, then we start to search for these template windows inside the second frame (post-compression) in a search window with the criteria discussed and we search row by row. After finishing the 1.sup.st row, we apply the lateral constrain by assuming that compression is mostly axial so the motion in the lateral direction is just due to poisson effect as described in equation (4).

.sub.lateral=.sub.axial (4)

[0026] This lateral strain is due to very small displacements, so we just determine the displacement of the pixel from the above row as pixel (x,y) and the point in the next row is (x,y+1) and let us assume that point (x,y) moved with displacement (dx,dy) where dx is the lateral displacement and dy is the axial displacement. Then by applying the lateral constrain on the searching window, the search window in the second frame will be centered at (x+dx,y+1+dy) with size equals to x+4 pixels the template window lateral size, so that it allows only searching x+2 pixels to the right and x+2 pixels to the left. Extending this criteria to the axial constrain, the search window will be equal to the template window axial size x+5, and here it starts to examine if this is a compression cycle or relaxation cycle (motion is upward or downward) and we search only in the direction of motion. Applying these 2 constrains has minimized the required computations as below

k=mn+1 (5)

where k is the number of operations, m the size of searching window, and n the size of template window. Using equation (5), if a template window size of 3232 is used along with a search window size with just 1.5 of template window size (4848),then we have k=17 axial operations and 17 lateral operations with total 289 operations for each tracking point. After applying axial and lateral displacement constrains for the same example we have template window with 3232 and search window with 3637, so we will have a total of 20 operations for each tracking point, i.e. the computational complexity reduced with about 14 times than the exhaustive search. After estimating the displacement of each row, we apply a 2D median filter to smooth and avoid any spurious displacement vectors.

[0027] After finishing the first level of enhanced HR tracking, displacement estimations between two frames will be obtained with high SNR and low resolution. This displacement field will be used to decide whether these frames underwent enough and valid compression to cause good contrast between tissues or not and that is done by extracting many features that have been studied and reported [4]-[6] besides what we have assumed based on Poission ratio that the later displacement should be symmetric around the axis of the axial compressions and others based on our observation like the entropy of the axial displacement and it's relation with the computed axial strain. To demonstrate that some features may affect the quality of the axial strain image and also features that we have found usefull according to what we have observed from WEKA. To take such decision, we have built a database containing many elastography videos that have been acquired using ultrasound breast phantom (Model 059, CIRS Inc., Norfolk, Va.) and many ultrasound machines of different manufacturers including, and not limited to, IBE, Siemens, Samsung, Philips, GE, Toshiba, Mindray, and Ultrasonix and then we have applied our algorithm over all frames with frame skipping factor equals one i.e. frames to be processed are at t & t+1 and selected good frames (qualitatively) manually (good frames are frames that produce high quality strain image) to extract features and then we have built a data set containing good and bad frames to start feature selection using information gain algorithm included in WEKA. After selecting frames we have trained many classifier with different number of features to see how accurate we can classify good frames from bad frames, we have used KNN, SVM, Nave bayes that are already implemented in WEKA.

[0028] For the second HR level, we apply procedures similar to the first level but with some hierarchy constrains. There is no need to estimate all the first row as we have done in the first level using exhaustive search, because we already know a priori info about its displacement field. We have a higher resolution 2nd level grid to track. First we obtain an initial estimate for the displacement at each point of interest in the 1.sup.st row of the new grid using the estimated displacement field of level 1 HR, and then the searching area is shifted with a step equals to estimated displacement. The size of this searching area becomes smaller than the one used in the first HR level and we use the small offset used in the axial and lateral constrains based on our calculations of the different displacements and strains from the 1.sup.st HR level that's implicitly a characteristic feature for the person applying this compression. For the subpixel accuracy in both levels, skilled in the art can use any algorithm based on his accuracy and computation speed requirements, we have used 3 point gaussian fit but still there are plenty of ways like 4-point gaussian estimator, interpolating the original images, center of mass, etc. Also, in addition to our constraints that avoids false peaks in the correlation field, we have used multi-correlation validation technique that is based on applying correlation equation several times on the same searching area but with little shifted or resized template window and multiply all correlation field with each others considering the difference between the template windows location and this can provide us with a correlation field with suppressed noise and amplitude true peak compared with other correlation fields.

[0029] After the second HR level is completed, a high resolution and accurate displacement field is reconstructed which can be used directly to estimate a high quality strain image. In one embodiment, we have used different techniques for strain estimation including the least-Square method [7], Savitzky-Golay filter [8], direct displacement gradient. All strain estimation techniques have produced good results since it depends mainly on the quality of the displacement estimation (strain is the gradient of the displacement).

[0030] Accumulation of axial strain to produce a higher quality image and increase the SNR. Accumulation depends mainly on the way we are computing strain if it is a lagrange strain then its characterized by higher contrast to noise ratio (CNR) at accumulation and if an Euler strain then its characterized by higher SNR that depends mainly on either we are tracking an object or a pixel. Last step, and this is also based on the physical assumptions that allow us to enhance our results, is applying adaptive normalization to the resulting strain image at which we move with a very wide window that normalize all tissue on the same horizontal level and this is to avoid stress leakage with the depth, since this would lead us to fake decisions thinking that there are hard tissue at high depth.

[0031] Also to be stated that the 2.sup.nd Hierarchical level can totally different than the first one according to matching equation (ZMCC, SAD, SSD), also it can be different algorithm/technique such as using optical flow for the 2.sup.nd Hierarchical level or using different different form or data to be tracked (register pre-compression image using the estimated displacement field). Also, we can use more number of Hierarchical level and do at each level something different than the previous one, with considering the results estimated from the previous one to enhance and speed up the estimation in the current level and so on until maximum accuracy is achieved at the last level of hierarchy.

[0032] While the present invention has been described in terms of particular embodiments and applications, in both summarized and detailed forms, it is not intended that these descriptions in any way limit its scope to any such embodiments and applications, and it will be understood that many substitutions, changes and variations in the described embodiments, applications and details of the method and system illustrated herein and of their operation can be made by those skilled in the art without departing from the spirit of this invention. Rather, the scope of the invention is defined by the following claims. We therefore claim as our invention all that comes within the scope and spirit of these claims.