Motion Estimation in Ultrasonic Images
In order to provide a suitable test of the reliability of this algorithm, experiments were performed using the synthesized images and the real ultrasonic images of muscle contractions and tendon motion. All of these experiments were performed on a Pentium 500MHz machine.Ⅲ
Accuracy Measurement
For synthesis images where the correct motion flows are known, the angular error measurement reported by Barron, Fleet, and Beauchemin [24] is used to evaluate the results. The error between the unit vectors of the estimated flow Ve and the correct flow Vc in the 3D space is:
.
Another choice of error measure is the magnitude of the 3-D difference vector, but it tends to have a higher value for large motion areas than for small motion areas. The angular error was chosen as the criterion for measurement accuracy.
In real motion cases the correct flow field is unavailable. The angular error is not applicable for accuracy measurement in real image sequences. In real motion cases, the peak signal-to-noise ratio (PSNR) is used as the accuracy measure. PSNR is based on the mean squared displaced pixel difference (dpd). It measures the correspondence of pixels described by the flow vectors between two frames. PSNR, in dB, is formulated as:
[ ]
target locations in the image plane.Synthesized Image Generation
Since a correct flow field is unavailable in real image sequences, tests with synthesis data sets are essential for performance evaluation. This measure is less meaningful if the properties of the synthesized data are very different from the properties in the real images. To generate synthesized
images with the correct flow field, a feature-line metamorphosis method [23] was adopted to produce the desired displacements between the pixel pairs in the two corresponding sonograms, which are then analyzed for measurement accuracy.
By giving a feature line pair, PQ and P'Q' as shown in Fig. 23 (a), pixel X with respect to
PQ is equivalent in geometrical relationship to pixel X’ with respect to P'Q'. We can establish the local coordinate systems in the source and reference images. The coordinate vector for each point is represented as (s,t) and (s’,t’) in the source image and the reference image respectively as in (32) and (33). In the source image, PR is the projection of PX on PQ, and s is the ratio of length of PR and PQ; and t is the projection distance from point X to PQ. To preserve the same geometrical relationship between point X and PQ for point X’ and P'Q', the coordinate vector of X with respect to PQ is (s,t) and the coordinate vector of X’ with respect to P'Q', denoted (s’,t’), is defined by the equation in (34).
2 where Perpendicular(⋅) returns the vector perpendicular to the input vector.
For cases with multiple feature lines, as shown in Fig. 23 (b), each feature line pair generates a corresponding point Xi’ for point X. The influences from each feature line pair are combined to obtain the corresponding point Xi’ with respect to X using the inverse distance-weighted interpolation, producing the corresponding point X’ using equations given in (2) and (3).
(a) (b)
Fig. 23. Feature line metamorphosis. (a) is the single feature line case, (b) is the multiple feature lines case
To generate the simulated ultrasonic image pair, the feature line pairs are assigned as shown in Fig. 24 (a). The arrows in Fig. 24 (a) represent the general motion defined by these feature line pairs. A backward mapping of the feature line morphing process is applied to generate the reference image as in Fig. 24 (b). A forward mapping of the feature line morphing process is
applied to generate the correct flow field as in Fig. 24 (c). Thus according to the desired motion, the simulated image pair and correct flow field can be generated. Since the reference image and the flow field were generated from a real case image, the simulated data set is closer to the real case than those generated with artificial patterns.
(a) (b) (c) Fig. 24. Generated data sets. (a) Source image with given feature line pairs,
(b) generated reference image, and (c) generated flow field from the given feature line pairs.
Block size selection
In the block matching technique, the size of the block is an important issue. A large block size holds more information but loses its local properties. A small block size maintains its local properties but may cause the aperture problem and fall into a local minimum. Computing time is another important consideration. A large block size always requires more time. An accuracy-efficiency curve (A-E curve) [25] was adopted to measure the performance of various block sizes. An exhaustive search block-matching algorithm was tested here. In our experiments with 232×240 ultrasonic images, the A-E curve in Fig. 25 indicates that a block sizes with a 15×15 or 21×21 pixel width produces optimal performance. Here, a block size with a 21×21 pixel width was chosen for further experiments.
21X21 41X41
51X51
31X31 25X25 15X15
11X11 7X7 5X5
0 2000 4000 6000 8000 10000 12000 14000 16000
0 5 10 15 20 25 30
Accuracy - Mean Anglar Error (degree) Efficiency - Computation Time (second)
Fig. 25. Accuracy-Efficiency Curve for block matching algorithm in different block sizes
Results from the synthesized images
In this section, the generated synthesis image pair is used to test the performance of the proposed method. In this simulation study, the difference between the original image and the motion compensation with correct flow is the optimal situation. In that situation, the optimal PSNR=36.48dB.
Since the proposed method is based on feature metamorphosis, the results are dependent upon the performance of the feature points. With feature points that are highly structured, a better correspondence and angular error can be achieved and thought of as an optimal result. In this case, the mean angular error=4.91(degree) was calculated only for the feature points.
In Table , the experimental results from the proposed method are compared with the Ⅱ exhaustive search block matching algorithm and the hierarchical method [15]. The computation time, PSNR, and the mean angular error are taken as the performance criteria. In these synthesis experiments, we can see that the proposed method performs better than the other methods in computation time, PSNR and angular error. The exhaustive search block matching method requires lots of time to produce a less accurate result and performs better in PSNR than the mean angular error since it does not take the smoothness constraint into account and only minimizes the intensity error to get the result. The hierarchical method reduces the computation time and performs better than the exhaustive block matching method in mean angular error. The proposed method can achieve a much better result using a similar computation time as the hierarchical method. Even without the refinement, the proposed method can still achieve a better result than the hierarchical method while saving computation time.
From Table , we can see the angular error before and after the refinement. In the proposed Ⅱ method, we can see that the angular error is increased by the effect of the non-feature points after the feature metamorphosis. Through the adjustment of the energy-based refinement, we can see that the result is improved significantly. The refinement process is not only helpful for the proposed method, but also helpful for other methods. As we can see from Table , the angular Ⅱ errors drop rapidly after applying the refinement process to both the hierarchical method and the proposed method.
TABLE Ⅱ
COMPARISON BETWEEN THE EXHAUSTIVE SEARCH METHOD, THE HIERARCHICAL METHOD AND THE PROPOSED METHOD. Computation
Time (s) PSNR (dB) Mean Angular Error (degree) Exhaustive Search Block Matching 4394 31.23 9.08
Hierarchical method 27.73 28.79 8.41
Hierarchical method with refinement 46.91 31.57 6.91 without
refinement 8.74 30.25 8.18
Feature metamorphosis
method with refinement 29.53 34.24 6.28
Results in real ultrasonic images
In order to evaluate this algorithm using living tissue, a muscle contraction experiment (Fig.
26) and a tendon motion experiment (Fig. 27) were performed. The muscle contraction was obtained by scanning along the longitudinal direction of a forearm in the proximal region while moving the middle finger and keeping the other fingers still. The tendon motion was obtained by scanning around the wrist area along the same direction using a similar finger motion. Since the correct flow is not available in the case of an actual image, we treated the PSNR as the measurement criterion. The motion compensated difference image was used to illustrate the flow estimation match error.
Two frames of the muscle contraction sequence with feature points and the corresponding flow upon it are illustrated in Fig. 26 (a) and (b). The dense motion vector field after feature metamorphosis is plotted in Fig. 26 (c), and the motion compensated difference image is shown in Fig. 26 (d) shows the displacement pixel difference. After the refinement process, we can obtain the dense motion vector field plotted in Fig. 26 (e), and the motion compensated difference image is shown in Fig. 26 (f). Although the actual displacements are not available, the estimated motion field is in agreement with the anatomical changes observed by clinical experts.
Figures 27 (a) (b) illustrate two frames of the tendon motion sequence with feature points and the corresponding flow. The dense motion vector field after feature metamorphosis is plotted in Fig. 27 (c) and the motion compensated difference image is shown in Fig. 27 (d). The dense motion vector field after the refinement process is plotted in Fig. 27 (e) and the resulting motion compensated difference image is shown in Fig. 27 (f). From the motion vector field, we can differentiate among the tendon and bone groups where the tendon group illustrates the movement from the anatomical changes and the bone group remains still.
Table shows the experimental PSNR results. Comparing the results in tⅢ he generated case, the PSNR results are acceptable since the motion in the actual case is much more complex than in the generated case. These results were also verified visually by doctors to be consistent with their clinical expectations.
(a) (b)
(c) (d)
(e) (f)
Fig. 26. A case measuring muscle motion. (a) is the source image, (b) is the reference image, (c) is the morph flow field, (d) is the motion compensated difference image of (c), (e) is the flow field after refinement, (f) is the motion
compensated difference image of (e)
(a) (b)
(c) (d)
(e) (f)
Fig. 27. A case measuring tendon motion. (a) is the source image, (b) is the reference image, (c) is the morph flow field, (d) is the motion compensated difference image of (c), (e) is the flow field after refinement, (f) is the motion
compensated difference image of (e)
TABLE Ⅲ
EXPERIMENT RESULTS OF ULTRASONIC IMAGES IN REAL CASES PSNR (dB)
The proposed method
PSNR (dB)
The hierarchical method
Muscle Case 27.64 24.34
Tendon Case 29.32 24.48
Motion Estimation in Tagged MR Images
To evaluate the performances of the proposed method, both the motion-synthesis phantom data sets and the clinical data sets of moving LV myocardium were applied to examine the accuracy. The grid pattern is generated by two orthogonal sets of tagging planes, which are spaced by a fixed 6mm interval. The obtained images (256 x 256 matrices, gray level digital data, 20 frames for each sequence) were analyzed on a Pentium III 1GMhz PC.
Motion-synthesis phantom experiments
In order to evaluate the accuracy of the tag tracking methods, a simulator based on a 13-parameter kinematic model of Arts et al. [55] has been proposed in [56] to simulate a sequence of tagged MR images. By setting the parameters of the model, the ground truth displacement field can also be generated with the simulated image sequence. The ground truth displacement field can be compared with the estimated displacement field for accuracy measurement.
In order to generate synthetic image sequences closed to the real situation, a motion-synthesis algorithm is proposed based on the phantom data sets. As shown in Fig. 28, the static phantom was first imaged using MR tagging schemes. Then, an affine transformation was applied to the phantom sequence to simulate the phantom deformation. The tag-tracking algorithm is then applied to the motion-synthesis sequence, and the ground truth displacement field can also be calculated with the affine parameters as the reference for evaluating motion estimation error.
Fig. 28. Illustration of generation of motion-synthesis images.
A good estimation on the strain field depends on the accuracy of the estimated displacement field. To evaluate the estimated displacement field, the PSNR (peak signal-to-noise ratio) is often used as the accuracy measure. PSNR is based on the mean squared displaced pixel difference (dpd). It measures the correspondence of pixels described by the displacement vectors between
the two image frames. The PSNR, in dB, is formulated as:
where Q is the total number of pixels, Λ denotes a set of target locations in the image plane, and dpd is the intensity difference between the corresponding pixels.
Another common measure criterion is the angular error measure, which is reported by Barron, Fleet, and Beauchemin [57]. The error between the unit vectors of the estimated flow Ve
and the correct flow Vc is given as: Fig. 29 shows an example of motion-synthesis phantom experiment where the phantom was stretched in the 45 degree direction. Angular error measurement is applied to evaluate the motion estimation. The value of PSNR and the optimal PSNR are also computed. The optimal PSNR indicates the optimal accuracy that the grey level deviation is computed with ground truth motion compensation. Fig. 29 (a) is the original phantom image at time=0, and Fig. 29 (b) is the motion-synthesis image at time=1 with estimated tag lines plotted over the images. The estimated displacement vectors are plotted on Fig. 29 (c), and the ground truth displacement vectors are plotted on Fig. 29 (d) where the color classification map is shown in Fig. 29 (e). Results of experiments on several different affine parameter tests are shown in Table IV. In this table, the angular errors of the proposed method are all smaller than four degrees. Also, the values of PSNR are spoiled less than 0.5dB when comparing with the optimal PSNR. The proposed algorithm performs very well for both measures in angular error and PSNR.
(a) (b) (c) (d) (e)
Fig. 29. Experimental results on motion-synthesis phantom images. (a) Tag lines in original phantom image P0. (b) Tag lines in motion-synthesis image P’1. (c) The estimated displacement field. (d) The ground truth displacement field generated from the affine parameters. (e) The color classification model where hue indicates the direction and
saturation indicates the magnitude of the displacement.
Table IV. ACCURACY EVALUATION IN MOTION-SYNTHESIS PHANTOM IMAGES Affine Parameters
(a0, a1, a2, a3)
Angular Error (degree)
PSNR (dB) of the Estimated Motion
PSNR (dB) of the Ground-Truth Motion
(1.05, 0, 0, 1.05) 2.21 18.86 19.27
(0.95, 0, 0, 0.95) 1.92 18.04 18.09
(1, 0.05, 0.05, 1) 2.84 18.49 18.69
(1, 0.05, -0.05, 1) 2.64 18.46 18.72
(1.05, 0.03, 0.03, 1.05) 3.06 18.93 19.39
(0.95, -0.03, -0.03, 0.95) 2.04 17.83 17.95
(1.05, 0.03, -0.03, 1.05) 2.05 18.97 19.05
(1.05, 0.03, -0.03, 0.95) 2.56 18.20 18.46
Results on clinical experiments
In the clinical experiments with the tagging schemes, several patients were scanned by targeting the LV in the base, middle, and apex slices. Displacement and deformation fields are computed to analyze on these clinical data sets. 3D motion analysis can also be extended from the proposed method easily. This can be achieved by integrating the longitudinal flow vectors from extra scanned long-axis view images. However, it is time consuming in scanning the dense 3D data set in both short-axis and long-axis views. To ease the load of patients, we scanned the LV only in base, middle, and apex slices.
To evaluate the accuracy of the clinical data sets, the MR delayed enhancement images [58, 59] are adopted as the reference. The delayed enhancement imaging is a useful technique in the assessment of the viability of the myocardium. The delayed enhancement sequence suppresses the signals from normal myocardium while demonstrating increased signals in the infracted areas of the myocardium where pooling of gadolinium contrast agents occurs. The clinical outcomes are examined by comparing the tag analysis results with the delayed enhancement images.
Fig. 30 shows an example of pathological case. The tag lines, which were initialized by the projection profiles of the edge image, were referred as the undeformed tag lines (as shown in Fig.
30 (a)). Fig. 30 (b) shows the tracked tag lines in the end systole period. Fig. 30 (c) shows the estimated displacement field, and Fig. 30 (d) illustrates the estimated principle strain map. From the displacement field, we can see that the myocardium contracts and rotates counterclockwise.
However, we can hardly indicate the undeformed regions due to the global translation. On the contrary, we cannot observe the motion from the strain map but we can inspect the deformation in each region at each time frame. In the strain map, we can find that the upper-left regions with low strain values are less contractile during the whole systolic period. These regions, which stop contracting while the rest of the myocardium contract normally, might be injured due to reduced blood flow. Comparing with the corresponding MR delayed enhancement image as shown in Fig.
30 (e), we can also find that the upper-left regions where are intensity enhanced are injured
regions. The strain map indicates the abnormal regions, which agree with the results of the delayed enhancement MR image.
(a) (b) (c)
(d) (e)
Fig. 30. (a) Initial tag lines as undeformed tag lines. (b) Tracked tag lines at end-systole period. (c) Estimated displacement field. (d) Estimated principle strain map and (e) the delayed enhancement MR image, where the plotted
arrows indicate the injured region of the myocardium.
Fig. 31 shows another experimental result of a patient case. The top, middle, and bottom rows show the results of the base, middle, and apex slices, respectively. The left column shows the tracked tag lines at end-systole. The middle column shows the estimated displacement fields.
The right column shows the estimated principle strain maps. In the strain maps, we can observe that the upper-left regions reveal low strain values. In this case, the upper-left regions of myocardium show less contractibility and are identified as the myocardial ischemia.
Fig. 31. Experimental results of a patient. From top raw to bottom raw are base, middle, and apex slices. Left column:
tracked tag lines at end-systole period. Middle column: estimated displacement fields. Right column: estimated principle strain map. Arrows indicate the abnormal regions.