 
 
IMPROVING TISSUE ELASTICITY IMAGING USING A KALMAN  
FILTER-BASED NON-RIGID MOTION TRACKING ALGORITHM  
 
 
 
 
 
By 
 
Susheel Reddy Vadde 
 
 
 
Submitted in Partial Fulfillment of the Requirements for 
 
The Degree of  
 
Master of 
 
Computing and Information Systems  
 
 
 
 
 
 
 
 
 
 
 
 
YOUNGSTOWN STATE UNIVERSITY 
 
May, 2011 
 
 
 
 
ii 
 
IMPROVING TISSUE ELASTICITY IMAGING USING A KALMAN  
FILTER-BASED NON-RIGID MOTION TRACKING ALGORITHM  
 
Susheel Vadde 
 
 
I hereby release this thesis to the public.  I understand that this thesis will be made 
available from the OhioLINK ETD Center and the Maag Library Circulation Desk for 
public access.  I also authorize the University or other individuals to make copies of this 
thesis as needed for scholarly research. 
 
 
Signature: 
     
  Susheel Reddy Vadde, Student  Date 
 
 
 
 
Approvals: 
     
  Dr. Yong Zhang, Thesis Advisor  Date 
 
 
 
 
     
  Dr. John Sullins, Committee Member Date 
 
 
   
 
     
  Dr. Kriss A. Schueller, Committee Member Date 
 
 
 
 
     
  Peter J. Kasvinsky, Dean of School of Graduate Studies and Research Date 
 
 
 
 
iii 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
© 
 
Susheel Vadde 
 
2011 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
iv 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
DEDICATION 
I dedicate this thesis to my parents who have encouraged me at all phases of my life  
and for their immense love on me. 
 
 
 
 
 
 
 
 
 
 
 
v 
 
ABSTRACT 
 
Imaging biomechanical properties of biological tissues under deformation 
has strong implications to many medical applications such as cancer detection and 
surgery planning. Quantifying the elasticity of soft tissue using an optical sensor is 
particularly attractive because it is non-invasive and easy to operate. However, the 
current computing method is plagued by the existence of noises in the two-frame 
optical flow solutions. In this thesis, a Kalman filter based tracking algorithm is 
examined, aiming to improve the quality of cumulative motion over a long video 
sequence. The proposed method is robust and is capable of handling non-rigid 
motion that is typical of soft tissue. Experiments of using videos of four rat tissue 
specimen subject to a biomechanical tensile test indicates that the proposed tracking 
method is very promising in generating a smooth, accurate, and continuous multi-
frame motion field. This type of multi-frame motion data not only allows us to 
compute a more accurate individual strain elastogram, but also provide valuable 
information for calibrating a series of relative strain images over the entire 
deformation process.    
  
 
 
 
 
 
 
 
vi 
 
ACKNOWLEDGEMENTS 
   
I would like to thank Dr. Zhang for giving me this opportunity to work on the 
biomedical imaging project and helping me to complete this thesis successfully. He has 
provided me tremendous encouragement and knowledge of how to approach a 
challenging problem by thinking beyond the box. 
I would like to thank Dr. Sullins for being a great advisor during my entire 
graduate study at YSU and for guiding me to choose courses of my interest. I also thank 
Dr. Schueller for taking time from his busy schedule to serve in the thesis committee.  
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
vii 
 
TABLE OF CONTENTS 
ABSTRACT  ....................................................................................................................... V 
ACKNOWLEDGEMENTS .............................................................................................. VI 
LIST OF FIGURES ........................................................................................................ VIII 
LIST OF TABLES  ........................................................................................................... IX 
1. INTRODUCTION ...........................................................................................................1 
2. RELATED WORK ..........................................................................................................2 
3. EXPERIMENT DESIGN .................................................................................................4 
               3.1 MECHANICAL TESTING ............................................................................4 
               3.2 VIDEO ACQUISITION AND PREPROCESSING  ......................................5 
4. COMPUTING METHODS  ............................................................................................8 
               4.1 TOW-FRAME MOTION BY OPTICAL FLOW ALGORITHM .................8 
               4.2 MULTI-FRAME MOTION BY ADDITION AND KALMAN FILTER  ...10 
               4.3 STRAIN IMAGING BY WEIGHTED GRADIENT FILTER  ....................13 
               4.4 COMPUTING TEST SETUP AND PARAMETERS  .................................14 
5. RESULTS AND DISCUSSIONS  .................................................................................19 
               5.1 MOTION TRACKING  ................................................................................19 
               5.2 STRAIN IMAGING  ....................................................................................23 
               5.3 EXISTING PROBLEMS  .............................................................................25 
6. CONCLUSIONS ...........................................................................................................26 
7. REFERENCES  .............................................................................................................27 
  
 
 
viii 
 
LIST OF FIGURES 
 
FIGURE 1. EXPERIMENT SETUP  .................................................................................. 4 
FIGURE 2. SAMPLES OF THE CROPPED FRAMES. ................................................... 7 
FIGURE 3. OPTICAL FLOW IMAGES .......................................................................... 20 
FIGURE 4. MULTI-FRAME MOTION FIELDS  ........................................................... 22 
FIGURE 5. THE STRAIN IMAGES ................................................................................ 24 
FIGURE 6. EXAMPLES OF STRAIN ELASTOGRAMS   ............................................ 25 
 
 
 
 
 
 
 
 
 
 
 
 
 
ix 
 
LIST OF TABLES 
 
TABLE 1. SPECIMEN VIDEOTAPING INFORMATION  .............................................. 6 
TABLE 2. VIDEO INFORMATION WITH TOTAL/BREAK TIME AND FRAME   ..... 6 
TABLE 3. KALMAN FILTER BASED TRACKING ALGORITHM  ............................ 12 
 
 
 
 
 
 
 
 
1 
 
1. INTRODUCTION  
 
Measuring elastic properties of biological tissues has found applications in a wide 
range of domains, particularly in cancer diagnosis and treatment planning. Various sensor 
modalities have been studied to provide quantitative measure of tissue abnormalities, 
including the ultrasound imaging, the magnetic resonance imaging, and the optical 
coherence tomography. Optical elastography is a relatively new modality that infers an 
object’s biomechanical property variation from the observed tissue surface deformation 
captured in an optical video sequence. Recently, optical elastography has been used in 
skin diseases analysis and surgery assessment. 
To generate an elastogram or strain image, an optical elastographic system must 
consist of three basic components: (i) A sensor that captures the tissue deformation; (ii) 
An algorithm that tracks the non-rigid tissue motion across multiple video frames; (iii) A 
forward or inverse algorithm that computes an absolute elasticity or relative strain image. 
One of the main difficulties in obtaining a high-quality elastogram is that the two-
frame motion (optical flow solution in this study) is often noisy, which has a severe 
adverse impact on elastogram quality, because the flow noise could be greatly amplified 
by a gradient based algorithm. The main contribution of this thesis is to utilize a Kalman 
filter to improve the multi-frame motion tracking and hence the elastogram quality.    
The proposed method was evaluated using the videos of four rat tissue samples 
underwent a tensile test. Experimental results suggest that, in comparison to a simple 
addition method, the Kalman filter based method can produce a smooth and dense motion 
field over a much longer frame length that is essential for elastogram calibration.  
 
2 
 
2. RELATED WORK  
 
There is a large body of literature on imaging tissue elasticity. In this section, a 
brief review of the state of the art of elastography technology is given, particularly the 
papers that are most relevant to this study. More comprehensive surveys and in-depth 
discussions can be found in [1, 2]. 
Since the earlier work of elastography [3, 4, 5], a large advance has been made in 
the algorithm development, sensor design, system integration, and clinical assessment, 
with a particular focus on early cancer detection and tissue characterization of internal 
organs, such as prostate, liver, breast, kidney, and heart [6, 7, 8, 9, 10]. The majority of 
the studies used either ultrasonic sensors or magnetic resonance imaging techniques. 
These types of modalities require direct contact of imaging sensors or excitation devices 
with a patient’s body to generate and gather the tissue deformation data. The force 
applied to the tissue could be in the form of dynamic mechanical waves or quasi-static 
loads. However, they may not be appropriate for certain clinical situations, for example, 
the skin patients who suffer burn damage or are in the post-surgery recovery period. 
Optical coherence tomography (OCT) has also been utilized to measure soft tissue 
biomechanical properties [11, 12]. An OCT based elastographic modality has the 
advantage that it has much higher resolution than the MRI based methods and less 
artifacts than the ultrasound based elastography. However, OCT based method has a 
drawback that its penetration depth is very limited and therefore its main application is in 
skin tissue characterization. Another problem of an OCT based method is that it only 
operates on a small surface area and also requires external excitation source to generate 
tissue deformation.  
 
3 
 
It is not until recently that researchers start to investigate the feasibility of in vivo 
imaging of skin elasticity using a regular optical camcorder [13, 14]. It has been 
demonstrated that vital information of skin tissue elasticity can be obtained with a high 
definition optical video acquired at a distance. Promising results have been observed in 
the cases of hypertrophic lichen planus, seborrheic keratosis, and psoriasis vulgaris. The 
optical strain images reveal the presence of significant tissue abnormalities associated 
with the underlying skin diseases. A similar technique based on finite element modeling 
and computer vision technique has also been applied to burn scar assessment [15].  
One of the major obstacles in optical imaging of tissue elasticity is that the motion 
data computed by an optical flow algorithm is often corrupted by noise, because the 
algorithm is very sensitive to the illumination variation that is difficult to avoid during 
video recording, even in a well-control laboratory setting. One possible solution is to 
compute a cumulative multi-frame motion field in the framework of Kalman filter.  
Kalman Filter is a mathematical tool for estimating a system’s state parameters 
given the system’s dynamic model, control inputs and noisy measurements [16]. Kalman 
filter has been successfully utilized in many tracking systems such as ballistic missiles 
and satellite navigation devices. In tissue property imaging [17], a method was developed 
using kalman filter and pulse echo ultrasound for detecting harmonic motion induced by 
ultrasonic radiation force. The method estimates the amplitude and phase of the motion at 
a tissue region with high sensitivity. In this study, the Kalman filter is mainly used to 
generate a more accurate multi-frame motion field by consolidating the noisy optical flow 
solutions. 
 
 
4 
 
3. EXPERIMENT DESIGN 
 
3.1 MECHANICAL TESTING 
 
Specimens of rat tissue were subject to a normal tensile test that determines the 
biomechanical properties. At first, a tissue sample was taken from the abdominal area of a 
rat and was cut into a dumbbell shape and preserved in the saline solution. An Instron 
Tensiometer (5500R) with a 100N load cell capable of 0.25% accuracy was used to 
stretch the specimen vertically. To avoid slipping effect, two grips were used to tighten 
the specimen to the pulling ends. The grips were also covered with aluminum foil to 
reduce potential contamination (see Figure 1).  
The location of the bottom grip was fixed and the top grip was pulled upward at a 
slow and constant speed. Duing a tensile test, a specimen exhibitded a series of sequential 
and gradual deformation patterns that can be characterized into three reprenstiative 
stages: Linear defromation, plastice deformation, and pre-rupture deformation.  
   
(a) Linear Deformation (b) Plastic Deformation. (c) Pre-rupture. 
        
Figure 1: The illustration of the tensile experiment setup and the rat tissue deformation. 
The bottom grip was fixed while the top grip was pulled upward steadily until the 
specimen was broken. Each specimen underwent an entire deformation process that 
includes the linear, plastic, and pre-rupture stages. 
 
 
5 
 
3.2 VIDEO ACQUISITION AND PREPROCESSING 
For each specimen, its entire deformation cycle was videotaped using a Sony 
HDR-SR1 camcorder with a capture speed of 30 frames per second. Recording condition 
includes a normal indoor lighting illumination and a fixed background setting in a 
mechanical engineering laboratory, with a constant room temperature of 25
o 
C. The 
camcorder was placed at a location that was about 1.5 meter from the tensiometer. 
However, the focus distance of the camcorder was adjusted slightly depending upon the 
size of a particular specimen and the study area of interest. For example, videos of certain 
specimens were recorded with a zoom-in mode to capture more details of tissue 
deformation. Once the camcorder was focused, tension force was applied to the specimen 
slowly until it broke. This procedure ensured that a smooth tissue motion was obtained to 
facilitate the subsequent optical flow and strain computations. Of the 41 specimens that 
were recorded, 4 representative ones were selected in the study of motion tracking and 
strain imaging. Table 1 summarizes the specimen and recording information. 
To perform motion tracking and strain computation, a video has to be pre-
processed on the basis of individual frames. At first, all videos were split into thousands 
of frames using the Sony Vegas Pro-8 software package. The process can use either the 
start/ending times or start/ending absolute frame numbers to specify a particular video 
segment. The output images were in the default JPEG format. The image aspect ratio 
(width/height ratio) was set to 1.0 square to ensure that there was no distortion involved 
during the frame extraction process. Table 2 lists the time/frame statistics of the 4 videos, 
and Figure 2 shows a few cropped frame samples. 
 
 
6 
 
Table 1:  Specimen Videotaping Information 
 
Specimen  
ID 
Recording 
Date 
Treatment 
Matrix 
Zoom-in 
View 
Break 
Location 
5A6I 9/20/2010 Yes No Top 
5A6S 9/20/2010 Yes No Above Matrix 
6A4S 9/21/2010 Yes Yes Top 
6A4I 9/21/2010 Yes No Above Matrix 
 
 
 
Table 2: Video Information with Total/Break Times and Frames 
Specimen  
ID 
Total # of 
Frames 
Break  
Frame 
Total  
Time 
Break  
Time 
5A6I 3479 3202 1min 56sec 1min 22sec 
5A6S 3420 3058 1min 54sec 1min 48sec 
6A4I 3914 3902 2min 20sec 2min 12sec 
6A4S 3074 2458 1min 42sec 1min 26sec 
 
 
 
7 
 
 
Frame 3450 Frame 3418 Frame 3902 Frame 3042 
 
Frame 8 Frame 12 Frame 5 Frame 1 
5A6I 5A6S 6A1I 6A4S 
 
Figure 2: Samples of the cropped frames that were used in computing optical flow and 
strain images. Note that the background does not have a significant effect on optical flow 
results because of the relatively slow tissue motion. 
 
 
 
8 
 
4. COMPUTING METHODS  
 
4.1 TWO-FRAME MOTION BY OPTICAL FLOW ALGORITHM 
 
Optical flow represents an apparent velocity field that describes an object’s motion 
observed in the image plane [18]. An optical flow field is obtained by solving a governing 
partial differential equation that is constrained by two fundamental assumptions: (i) The 
brightness of a material point remains unchanged between two image frames; (ii) All 
material points in a relatively small image neighborhood move at the same or similar 
speed and direction. 
Considering the image brightness as a function of spatial coordinates and time, 
E(x,y,t), then the brightness constancy between two images can be expressed as a total 
derivative using  the well-known conservation principle: 
 
                                     
.0  
dt
dE
                                                                (1) 
 
For the majority of motion tracking problems, the brightness function can be 
treated as differentiable in both spatial and temporal domains and the equation (1) can be 
expanded to: 
 
      
 ((), (),)
0.
dE x t y t t E dx E dy dE
dt x dt y dt dt
 w w
      
 w w
                                 (2) 
 
In an image plane, a motion vector has two components: U = [u, v]
T
 , with u = 
dx/dt for the horizontal motion and v = dy/dt for the vertical motion.  Using E
t 
to denote 
the temporal derivatives of the brightness function, the brightness constancy equation can 
be written in a more compact form: 
 
 
9 
 
                                
.0)(     �
t
T
EUE
                                                    (3) 
 
It should be noted that the PDE as expressed in equation (3) cannot be solved 
without ambiguity due to the lack of higher level knowledge of how an object moves in 
the 3D space, i.e., the aperture dilemma. To solve this problem, the brightness constancy 
equation must be regularized with a smoothness prior term: 
 
   
 
22222
(,) ( ) ( ),
x yt xyxy
obj uv Eu Ev E u u v v O        
                   (4) 
 
Where E
x
, E
y
, u
x
, u
y
, v
x
 and v
y
 represent the partial derivatives of the corresponding 
variables, and   is the Lagrangian multiplier (regularization coefficient). An optical flow 
solution (u, v) between two images can then be computed by minimizing the objective 
function (4). It is clear that this is a constrained solution that represents a compromise 
between the observed motion and the smoothness penalty. 
In many situations, the two fundamental assumptions underlying the optical flow 
equation are often violated, especially in the presence of moving physical boundaries and 
depth discontinuities such as the one between a tissue specimen and the background 
(Figure 2). To deal with this problem, a robust algorithm is utilized that includes the 
piecewise smoothing and multi-resolution methods [19]. The breach of brightness 
constancy and spatial smoothness conditions is handled by an adaptive smoothing 
strategy that adjusts the penalty weight depending upon the local discontinuity property. 
The actual algorithm implementation can be found in [20]. In all of the experiments of 
this study, the optical flow parameters were pre-set based on the recommended values. 
But the iteration number was increased to improve the convergence performance because 
of the image size. 
 
 
 
10 
 
4.2 MULTI-FRAME MOTION BY ADDITION AND KALMAN FILTER 
 
In optical flow computation, two adjacent frames were used to derive the pixel level 
motion field in order to satisfy the small displacement requirement of the governing 
equation. But the resulting two-frame flow was too small to be used in strain imaging, 
because the flow noise could be largely amplified by a gradient strain filter. Therefore, a 
tracking method is needed that calculates a larger cumulative motion by integrating 
multiple two-frame flows over a longer sequence. 
Two tracking methods were considered: (a) A simple addition method that adds up 
flow results at the fixed pixel locations; (b) A robust method that is implemented in the 
framework a Kalman filter.   
In the first method, the horizontal and vertical components of the motion vector 
computed from multiple frame pairs are added at each pixel:  
112 1( , ) ( , ) ( , )  ... ( , )ui i n ui i ui i ui n i n              
    (5) 
112 1( ,  ) ( ,  ) ( ,  )  ... ( ,  )vi i n vi i vi i vi n i n              
 
   (6) 
Where u(i, i+n) and v(i, i+n) denote the cumulative motion between Frame i and Frame 
i+n. The cumulative motion vector will then be used to compute strain elastograms. It is 
difficult to estimate an optimal value of n because it depends on both an object’s motion 
speed and the camera capture speed. In this thesis, an empirical value of 15 was used 
based on the results of a number of evaluation tests. 
The addition method has the advantage that it is easy to implement and can process 
a large number of frames fairly efficiently. However, the method is based on the 
assumption that the motion of a physical point is relatively small, which is only valid for 
 
11 
 
a short sequence. For example, with the videos of rat tissue, 10-50 frames usually 
produce good tracking results. But for a longer sequence (>100 frames), the method will 
introduce significant artifacts and lead to poor strain images. 
The Kalman filter was initially derived from linear dynamic system over discrete 
time domain. The state of system is defined as a vector of real numbers. A Kalman filter 
calculates at each discrete time interval to generate a new state from the present state. 
Linear operators are applied for each state with some noise measurements included.   
In this study, a Kalman filter was used to integrate noisy two-frame optical flows. 
The method has the following advantages: (i) The measurement uncertainties (optical 
flow noise) is handled in a recursive prediction-correction manner that ensures an optimal 
tracking; (ii) The algorithm can work with a large number of points to produce a dense 
motion field, which is necessary for generating a strain image with a resolution close to 
the original input images. 
The algorithm will produce an optimal estimate of a state vector by continuously 
updating the predictions of a process mode. It is optimal in the sense that the error 
covariance is minimized. A state vector (s
t
) is defined as: s
t
 = [x
t
, y
t
, u
t
, v
t
]
T
, where (x
t
, y
t
) 
and (u
t
, v
t
) are the true position and the true motion of a point at time t. A measurement 
vector (m
t
) is defined as: m
t
 = [x*
t
, y*
t
, u*
t
, v*
t
]
T
, where  (x*
t
, y*
t
) is the position 
estimates and (u*t, v*t) is an optical flow solution. The multi-frame tracking is 
implemented with a recursive filter consisting of a process model and a measurement 
model. The process model describes the evolution of the state vector: s
t
 = As
t-1 
+ q
t-1
, p(q) 
~ N (0,Q), where A is the state transition matrix, and q is a zero-mean Gaussian noise 
with a covariance matrix Q. The measurement model relates the state vector to the noisy 
 
12 
 
flow vector: m
t
  = Hs
t
 + r
t
, p(r) ~ N (0,R), where H is the identity measurement matrix, 
and r is the measurement noise with a zero-mean Gaussian distribution and a covariance 
matrix R. 
The tracking alternates between two stages: (i) In prediction stage, the state vector 
at the time t is estimated using the values obtained at the previous time t-1 by: �
t/t-1
 = A�
t-1
. 
As the measurement is not used at this stage, �
t/t-1 
can be viewed as the a priori estimate. 
The a priori error covariance matrix is computed by: P
t/t-1 
= AP
t/t-1
A
T
 + Q
t-1
, which 
quantifies the estimation uncertainty; (ii) In the correction stage, the Kalman gain matrix 
is computed by: K
t
 = P
t/t-1
H
T
(HP
t/t-1
H
T
 + R
t
)
-1
, which is a blending factor that minimizes 
the a posteriori error covariance. With a blending factor, the a posteriori estimate of the 
state vector is updated by the measurement: �
t
 = �
t/t-1 
+ K
t 
(m
t
  - H�
t/t-1
). The a posteriori 
error covariance matrix is also updated by: P
t
 = (I - K
t
H)P
t/t-1
. The overall Kalman filter 
based tracking procedure is summarizing in Table 3. 
 
Table 3.  Kalman filter based multi-frame tracking algorithm. 
Stage Equation Interpretation 
Prediction 
�
t/t-1
= A�
t-1
Predict the state (a priori estimate) 
P
t/t-1
= AP
t/t-1
A
T
 + Q
t-1
Predict the error covariance (a priori estimate) 
Correction 
K
t
= P
t/t-1
H
T
(HP
t/t-1
H
T
 + R
t
)
-1
Compute the Kalman gain matrix 
�
t
= �
t/t-1
 + K
t
(m
t
  - H�
t/t-1
) Update the state vector with measurement (a posteriori estimate) 
P
t
= (I - K
t
H)P
t/t-1
Update the error covariance (a posteriori estimate) 
 
 
 
13 
 
4.3 STRAIN IMAGING BY WEIGHTED GRADIENT FILTER 
 
 Given a multi-frame motion field, a strain elastogram can be computed using a 
Cauchy strain tensor (0 ): 
                   
1
()
2
1
()
2
xx yx
xy yy
uuv
x
vu v
xy x
 H H
 H
 H H
 � � w w w
     
  
 w w  w
      
 w w  w
 ,                           (7) 
where (0
xx
, 0
yy
) are normal strains and (0
xy
, 0
yx
) are shear strains. Since the motion vector is 
presented as optical flow images, spatial derivatives of the motion vector is conveniently 
obtained through a convolution operator: 
* ,       * ,   
xx
uv
Su Sv
 w w
    
                                (8) 
* ,       * ,   
yy
uv
Su Sv
yy
 w w
    
                                   (9) 
where S
x
 and S
y
 are the weighted Sobel filters. Using the spatial derivatives of the motion 
vector, a strain image can be computed for each strain tensor component (0
xx
, 0
yy
, 0
xy
, 0
yx
). 
If both horizontal and vertical strains are large, a magnitude elastogram can be obtained:   
 .      
2222
mxxxyyxy
 H  H H H H     
                                             (10)         
 In this study, a specimen was mainly pulled vertically. As a result, the vertical 
strain elastogram (0
yy
) carries much more information about a sample’s biomechanical 
property variation. For a better visualization, a strain image was also normalized to a grey 
scale of [0, 255], with brighter/darker pixels for higher/lower strains values.   
 
 
14 
 
4.4 COMPUTING TEST SETUP AND PARAMETERS 
 
 A test was carried out in six steps: (1) Frame extraction; (2) Image conversion; (3) 
Image cropping; (4) Computing optical flows; (5) Multi-frame motion tracking; (6) 
Computing strain elastograms. The most computation intensive steps of optical flow and 
motion tracking were implemented with C/C++ and the high-level frame scheduling was 
managed using Perl.  All tests were conducted in a Cygwin or Linux environment. 
 
(1) Frame Extraction 
1. Absolute frames of a video sequence were extracted using Sony Vegas Pro 8. 
2. Parameters: HDV 1080-60i (1440 x 1080, 29.97 fps, pixel aspect ratio: 1.0). 
3. The average time of extracting an entire video is about 90 minutes. 
(2) Image Conversion 
1. The extracted image frames are in the JPEG format. 
2. JPEG images were converted to PGM format using a Perl script. 
 
 
 
 
 
 
 
 
 
 
 
 
 
#!/usr/bin/perl 
$src = "./6A_1I"; 
$dst = "./6A_1I_pgm"; 
opendir(DIR, "$src"); 
@allfiles = readdir(DIR); 
$i = 0; 
foreach $file (@allfiles) { 
@base = split (/_/, $file); 
if($file =~ /jpeg/) { 
system("convert  $src/$file  $dst/$base[0]$i.pgm"); 
print  "converting $file to  $base[0]$i.pgm \n"; 
$i = $i + 1; 
} 
} 
 
15 
 
(3) Image Cropping 
1. Purpose:  remove the background that is irrelevant to tissue strain imaging.   
2. Inputs:  the coordinates of two controlling corner points. 
 
 
 
 
 
 
 
 
 
 
(4) Computing Optical Flow 
1. This is the most time consuming step. 
2. Inputs:  cropped images indexed by two integer numbers. 
3. Key parameters: Iteration number = 20, Stages = 5, Frame interval = 1. 
4. Outputs: Image/text files of motion vector in horizontal/vertical directions.  
Specimen  No of Frames 
No. of Frames 
Tracked 
Computing Time 
5A6I 3479 800 14 hrs 20 min 
5A6S 3420 800 14 hrs 7 min 
6A4I 3712 800 14 hrs 37 min 
6A4S 3074 800 14 hrs 25min 
 
 
!/local/bin/perl 
#--------------------------- 
$inp = "."; 
$out = "../6A_1I_pgm_cropped/1_1000/"; 
#--------------------------- 
$start = 1; 
$end   = $start + 1000; 
#--------------Crop sub image --------------------------- 
for($i = $start; $i < $end; $i++){ 
    $name = "Image$i"; 
    system("../../../Rat_Code_Crop/crop  $name.pgm  708 386 788 918 "); 
    print "cropping $name \n"; 
} 
system("mv  cropped*.pgm $out"); 
print "Finish moving images. \n"; 
 
16 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
#!/local/bin/perl 
 
$rows = 791; 
$cols = 171; 
 
$interval = 1; 
$skip     = 0; 
 
$start = 1; 
$end   = $start + 20; 
#-------------------------------------------------- 
$src = "./6A_1I_pgm_cropped/1_20"; 
$dst = "./6A_1I_pgm_cropped_flow/1_20"; 
#-------------------------------------------------- 
for($i = $start; $i < $end; $i++){ 
    $j = $i + $interval; 
    system("../../Rat_Code_Flow/src/gnc $i $j 4 1  $src/cropped_Image  $dst/flow-
$i-$j- -l1 10.0 -l2 1.0 -S1 10.0 -S2 1.0 -s1 10.0 -s2 0.05 -stages 5 -nx $cols -ny 
$rows -f 0.5 -F 0 -w 1.995 -iters 20 -end '.pgm' -skip 15");  
    print  "optical flow:  $i  $j  \n";  
    $i = $i + $skip; 
} 
#----------- Delete unused files ------------------------------ 
system("rm $dst/*-0.*"); 
system("rm $dst/*-1.*"); 
system("rm $dst/*-2.*"); 
system("rm $dst/*-3.*"); 
system("rm $dst/*-data-*"); 
system("rm $dst/*-spatial-*"); 
system("rm $dst/*-stable-*"); 
system("rm $dst/*.hdr*");       
 
#----------- Move Text files ------------------------------ 
system("mkdir $dst/text"); 
system("mkdir $dst/u_flow"); 
system("mkdir $dst/v_flow");   
 
system("mv $dst/*.txt     $dst/text/"); 
system("mv $dst/*-u-*.pgm $dst/u_flow/"); 
system("mv $dst/*-v-*.pgm $dst/v_flow/");   
 
17 
 
(5)  Motion Tracking 
1. Both Addition method and Kalman filter method were used. 
2. Inputs: multiple optical flow fields to be tracked. 
3. Parameters: boundary width (BD), smoothing window, starting/ending frames.  
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
#!/local/bin/perl 
 
$rows = 791; 
$cols = 111; 
$skip = 0; 
$BD                 = 40; 
$nSmooth        = 1; 
$smooth_halfWin = 2; 
 
#-------------------------------------------------- 
$src = "./6A_1I_pgm_cropped_flow/1_2000/text"; 
$dst = "./6A_1I_pgm_cropped_flow_track/1_2000"; 
system("mkdir $dst"); 
#-------------------------------------------------- 
 
$startID = 500; 
$endID   = $start + 995; 
 
#------------------------------         exe  src  dst  startID  endID  rows  cols  BD  
nSmooth  smooth_halfWin stage 
system("../../Rat_Code_TRACK_Flow/TRACKFLOW $src $dst $startID 
$endID $rows $cols $BD $nSmooth $smooth_halfWin 4    "); 
 
#----------- Move files --------------------------- 
system("mkdir $dst/track_text"); 
system("mkdir $dst/track_u"); 
system("mkdir $dst/track_v"); 
 
system("mv $dst/*.txt    $dst/track_text/"); 
system("mv $dst/*-u*.pgm $dst/track_u/"); 
system("mv $dst/*-v*.pgm $dst/track_v/");   
#--------------------------------------------------  
 
 
18 
 
(6)  Computing Strain Images 
1. The strain image was computed using a gradient filter. 
2. Inputs: multi-frame motion field. 
3. Parameters: FDorSobel: 1 for finite difference method and 2 for Sobel filter. 
4. Outputs: a strain elastogram. 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
#!/local/bin/perl 
$rows = 791; 
$cols = 171; 
$exe    = "../../Rat_Code_Strain/STRAIN"; 
$inpdir = "./5A_6S_pgm_cropped_flow/3000_4000/addflow/addflow_text"; 
$outdir = "./5A_6S_pgm_cropped_strain"; 
system("mkdir $outdir"); 
$nSmooMot  = 0;  
$nSmooStr  = 0; 
$ConWidth  = 1; 
$halfWin   = 1; 
$FBD       = 50;    
$SBDTOP    = 210; $SBDBOT    = 20; 
$SBDLEFT   = 20; $SBDRIGHT  = 15; 
$FDorSobel = 1; 
#-------------------------------------------------- 
opendir(DIR, "$inpdir"); 
@allfiles = readdir(DIR); 
  
foreach $file (@allfiles) {  
     if($file =~ /-U/) { 
        @token = split(/\./, $file); 
        $base1 = $token[0];  $base2 = $token[0]; $base2 =~ s/-U/-V/; 
        #        exe  base1  base2  inpdir  outdir  rows  cols  nSmooMot  nSmooStr  
ConWidth  halfWin  FBD  SBDTOP  SBDBOT  SBDLEFT  SBDRIGHT  
FDorSobel 
        system("$exe $base1 $base2 $inpdir $outdir $rows $cols $nSmooMot 
$nSmooStr $ConWidth $halfWin $FBD $SBDTOP $SBDBOT $SBDLEFT 
$SBDRIGHT $FDorSobel"); 
        print  "$exe $base1 $base2 $inpdir $outdir $rows $cols $nSmooMot 
$nSmooStr $ConWidth $halfWin $FBD $SBDTOP $SBDBOT $SBDLEFT 
$SBDRIGHT $FDorSobel \n\n";  
     } 
} 
 
19 
 
5. RESULTS AND DISCUSSIONS 
 
5.1 MOTION TRACKING 
 
The two-frame motion fields computed with the optical flow algorithm are quite 
noisy (see Figure 3). In the vertical flow images, darker pixels indicate the upward 
motion and brighter pixels indicate downward motion. In the horizontal flow images, 
darker pixels suggest the motion to the left and brighter pixels are for the motion to the 
right. It should be noted that the motion values at all pixels in an image are normalized to 
a grayscale of 0 to 255. Therefore, a flow image is of relative meanings only. In other 
words, a vertical flow image and a horizontal flow images cannot be compared on the 
basis of their absolute pixel intensity values. In fact, since the tissue deformation in a 
tensile test is dominated in the vertical direction, the magnitude of vertical motion is 
much larger than that of its horizontal counterpart. 
Several observations can be made from the two-frame optical flow results: (i) The 
motion vector has been successfully computed (both horizontal and vertical components), 
although data quality is compromised by flow noise; (ii) Since the difference between the 
horizontal and vertical components could be of several orders of magnitude, the strain 
images will be primarily influenced by the vertical component; (iii) Computing an optical 
flow using the two adjacent frames satisfies the small displacement condition of the 
governing equation. However, the motion values are in the sub-pixel level (< 0.4 pixel), 
which makes it vulnerable to random noise. One possible solution is to use larger frame 
interval without violating the constraint condition, which will be investigated in future 
studies; (iv) Including the background in the cropped images does not seem to have a 
significant impact on the flow results, probably due to the small two-frame deformation. 
 
 
20 
 
 
Frame 1022 Frame 1023 Frame 2045 Frame 2046 
 
 
 
Horizontal Motion 
(Frames 1022-1023) 
Vertical Motion 
(Frames 1022-1023) 
Horizontal Motion 
(Frames 2045-2046) 
Vertical Motion 
(Frames 2045-2046) 
 
Figure 3: Sample frames from the videos of specimen 6A1I and the corresponding two-
frame optical flow results. The geometry of the specimen can be roughly identified in the 
flow images that are generally noisy.  
 
 
21 
 
Figure 4 shows the multi-frame tracking results using the addition method and the 
Kalman filter based method. The darker pixels indicate upward or leftward motion and 
brighter pixels indicate downward or rightward motion, respectively, with the intensity 
values proportional to the motion magnitude. 
For the addition method (Figure 4, top row), 15 frames were added to generate a 
multi-frame motion. It is apparent that the motion quality is much improved as evidenced 
by the smoother image intensity, in comparison to the two-frame flow images. If a flow 
image can be modeled as motion signal plus the white noise, adding multiple flow images 
at the fixed pixels clearly averaged out the random noises. Because the two-frame flow 
results show much larger noises at the boundary, 30-50 boundary pixels were cut off to 
minimize the noise effect (the black pixels). 
 Large improvements were also observed in the multi-frame motion images 
computed using the Kalman filter algorithm (Figure 4, bottom row).  Given a short video 
segment (frames 1-16 for the addition method and frames 1-33 for Kalman filter method), 
both approaches produced good results. However, only the Kalman filter method is 
capable of handling a long sequence that captured the entire tissue deformation (linear 
stage, plastic stage, and the pre-rupture stage). For example, the cumulative motion of 
more than 3,000 frames (6A1I) has been tracked with a good quality. This type of 
tracking over a large number frames is very useful in matching the stress-strain history, 
which can be further utilized in strain calibration. Some of the feature points were lost 
during the Kalman filter tracking. This problem can partially be solved by recovering the 
lost points (e.g., assigning the average motion values of its neighboring points).     
 
 
22 
 
 
Horizontal Motion 
(Frames 1-16) 
Vertical Motion 
(Frames 1-16) 
Horizontal Motion 
(Frames 1001-1016) 
Vertical Motion 
(Frames 1001-1016) 
 
 
 
Horizontal Motion 
(Frames 1-33) 
Vertical Motion 
(Frames 1-33) 
Horizontal Motion 
(Frames 1000-1300) 
Vertical Motion 
(Frames 1000-1300) 
 
Figure 4: Multi-frame motion fields computed using a simple addition method (top row) 
and a Kalman filter based robust tracking algorithm (bottom row). 
 
 
23 
 
5.2 STRAIN IMAGING 
 
 Given the multi-frame motion fields tracked by either the addition method or the 
Kalman filter method, strain elastograms can be obtained by applying the weighted 
gradient filter to the motion data (see Figure 5).    
 For a short sequence (15 frames and 33 frames), both tracking methods resulted in 
good strain images of similar patterns that characterize the tissue matrix. For example, 
the muscle structure (layers of fibers) is clearly visible, which is indicative of the 
existence of material property heterogeneity that is difficult to see in the raw images. Due 
to the lack of quantitative measurement of collagen and elastin fibers of the specimen, it 
is not clear yet whether the observed strain patterns can be correlated to the property 
variations in different deformation stages (linear, plastic, and pre-rupture). In future 
investigations, it is highly desirable to establish a quantitative correlation based on the 
tests of a large number samples, so that a strain elastogram can be interrelated with a 
sound biophysical support, at least on the level of collagen and elastin fibers that play the 
major role in determining tissue elasticity. 
 The strain images of a long sequence (300 frames) tracked by the Kalman filter is 
somewhat different from that of a short segment (15 frames) by the addition method, with 
the former showing more disparity in strain values at the top and middle sections, which 
is more descriptive of the tissue deformation history over time. It should be emphasized 
that both types of strain elastograms provide valuable information of tissue property 
changes associated with the pathological process and hence have great potentials in 
revealing diseased induced tissue abnormalities. 
 
24 
 
  
Horizontal Strain 
(Frames 1-16) 
Vertical Strain 
(Frames 1-16) 
Horizontal Strain 
(Frames 1001-1016) 
Vertical Strain  
(Frames 1001-1016) 
 
 
 
 
 
Horizontal Strain 
(Frames 1-33) 
Vertical Strain 
(Frames 1-33) 
Horizontal Strain 
(Frames 1000-1300) 
Vertical Strain  
(Frames 1000-1300) 
 
Figure 5: The strain images derived from the multi-frame motion tracked by the addition 
method (top row) and a Kalman filter algorithm (bottom row). All images were computed 
using the videos of 6A1S. 
 
 
25 
 
5.3 EXISTING PROBLEMS 
 
Since optical elastography is still in its early developing stage, there are many 
technical issues that need to be addressed before it can be used in a clinical trial and 
eventually real medical applications. A large number of evaluation tests have to be 
conducted using various tissue samples.  In this study, it was found that, for certain video 
frames, the strain elastograms using the tracking data of the addition method had artifacts 
in the forms of salient irregular curves (Figure 6, image on the right). It is not clear what 
causes the artifacts that might be related to the smoothing operators or boundary effect. 
Another problem is that the feature points could be lost during Kalman filter tracking as 
seen in Figure 4.  All of those issues need to be addressed in future investigations. 
 
Vertical Strain 
(Frames 21-36) 
Vertical Strain 
(Frames 41-56) 
 
 
Figure 6: Examples of strain elastograms with and without artifacts using the addition 
method (the vertical strain component of frames 41-56 from specimen 6A4S). 
 
 
 
26 
 
6. CONCLUSIONS 
This thesis investigates the feasibility of using a Kalman filter based approach to 
improve the quality of motion tracking over a long video sequence and the subsequent 
strain elastogram computation. The experiments were carried out using the optical videos 
of rat tissue that were subject to a normal biomechanical tensile test. Given the noisy two-
frame optical flow solutions, the proposed Kalman filter based tracking algorithm proves 
to be a very promising approach in minimizing the impact of flow noise and providing a 
stable and accurate estimation of tissue deformation. In comparison to the addition 
method, the Kalman filter method has the advantage that it is capable of tracking a very 
long video sequence without causing a severe deviation. More importantly, the strain 
elastograms derived from the long tracked multi-frame motion field has significant 
implications to both basic medical imaging research and clinical applications.  
  
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
27 
 
7. REFERENCES 
[1] Ophir J, Alam SK, Garra BS, Kallel F, Konofagou EE, Krouskop T, Merritt 
CRB, Righetti R, Souchon R, Srinivasan S, and Varghese T, Elastography: 
Imaging the elastic properties of soft tissues with ultrasound, Journal of Medical 
Ultrasonics, 2002; 29 (Winter): 155-171. 
[2] Parker KJ, Doyley MM, and Rubens DJ, Imaging the elastic properties of tissue: 
the 20 year perspective, Physics in Medicine and Biology, 2011; 56: R1-R29. 
[3] Ophir J, Cespedes I, Ponnekanti H, Yazdi Y, and Li X, Elastography: A 
quantitative method for measuring the elasticity of biological tissues, Ultrasonic 
Imaging, 1991; 13: 111-134. 
[4] Parker KJ, Huang SR, Musulin RA, and Lerner RM, Tissue response to 
mechanical vibrations for “sonoelasticity imaging”, Ultrasound in Medicine and 
Biology, 1990; 16: 241-246. 
[5] Muthupillai R, Lomas DJ, Rossman PJ, Greenleaf JF, Manduca A, and Ehman 
RL, Magnetic resonance elastography by direct visualization of propagating 
acoustic strain waves, Science, 1995; 269: 1854-1857. 
[6] Krouskop TA, Wheeler TM, Kallel F, Garra BS, and Hall TJ, The elastic moduli 
of breast and prostate tissues under compression, Ultrasonic Imaging, 1998; 20: 
260-274. 
[7] Righetti R, Kallel F, Stafford RJ, Price RE, Krouskop TA, Hazle JD, and Ophir 
J, Elastographic characterization of HIFU-induced lesions in canine livers, 
Ultrasound in Medicine & Biology, 1999; 25: 1099-1113. 
 
28 
 
[8] Garra BS, Céspedes I, Ophir J, Spratt S, Zuurbier RA, Magnant CM, and 
Pennanen MF, Elastography of Breast lesions: initial clinical results, Radiology, 
1997; 202: 79-86. 
[9] Weitzel WF, Kim K, Rubin JM, Xie H, and O’Donnell M: Renal advances in 
ultrasound elasticity imaging: measuring the compliance of arteries and kidneys 
in end-stage renal disease, Blood Purification, 2005; 23:10-17. 
[10] Konofagou EE, D’hooge J, and Ophir J, Myocardial elastography, A feasibility 
study in vivo, Ultrasound in Medicine and Biology, 2002; 28: 475-482.   
[11] Ko HJ, Tan W, Stack R. Boppart SA, “Optical Coherence Elastography of 
Engineered and Developing Tissue”, Tissue Engineering, 2006; 12: 63-73. 
[12] Sun C, Standish B and Yang VXD, Optical coherence elastography: current 
status and future applications, J. Biomed. Opt. 16, 043001 2011. 
[13] Zhang Y., Brodell R. T., Mostow E. N., Vinyard C. J., and Marie H., In vivo skin 
elastography with high definition optical videos, Skin Research and Technology, 
2009; 15: 271-282. 
[14] Zhang Y, Sullins JR, Goldgof DB, Manohar V, Computing strain elastogram’s of 
skin using an optical flow based method, Proceed. of 5th Inter. Conf. on 
Ultrasonic Measurement and Imaging of Tissue Elasticity, Snowbird, Utah, 
October, 2006. 
[15] Zhang Y, Goldgof DB, Sarkar S, Tsap LV. A modeling approach for burn scar 
assessment using natural features and elastic property, IEEE Transactions on 
Medical Imaging, 2004; 10:1325-1329. 
 
29 
 
[16] Kalman RE, A New Approach to Linear Filtering and Prediction Problems, 
Transactions of the ASME - Journal of Basic Engineering, 1960; 82: 35-45.  
[17] Yi Zheng, Shigao Chen, Wei Tan, Kinnick R. Greenleaf JF,  Detection of tissue 
harmonic motion induced by ultrasonic radiation force using pulse-echo 
ultrasound and kalman filter, IEEE Transactions on Ultrasonics, Ferroelectrics 
and Frequency Control,  2007: 54: 290-300. 
[18] B. K. P. Horn and B. G. Schunck, Determining optical flow, AI Memo 572, 
Massachusetts Institute of Tech., 1980. 
[19] M. J. Black, P. and Anandan, The robust estimation of multiple motions: 
parametric and piecewise-smooth flow fields, Computer Vision and Image 
Understanding, 1996; 63: 75-104. 
[20] http://www.cs.brown.edu/~black/