Medical Image Registration and Visualization on Tumor Growth with Time Series

Kuocheng Wang; Thenkurussi Kesavadas

doi:10.47739/1005

Medical Image Registration and Visualization on Tumor Growth with Time Series

Research Article | Open Access | Volume 1 | Issue 1

Article DOI : https://doi.org/10.47739/1005

Kuocheng Wang Thenkurussi Kesavadas

^1. Department of Industrial and System Engineering, University of Illinois at UrbanaChampaign, USA

+ Show More - Show Less

Corresponding Authors

Kuocheng Wang, Department of Industrial and System Engineering, University of Illinois at Urbana-Champaign, Urbana, IL, USA Tel: 217-377-9982;

Abstract

Data registration is a common process in medical image analysis. The goal of data registration is to solve the transformation problem with multiple images’ alignment. Conventionally, diagnosing the tumors periodically requires understanding the growth and spread of tumor which is performed by doctors by visual inspections of multiple MRI scans taken over different stages in time series. Due to the misalignment of patient’s posture, comparison of these multiple MRI scans is tedious. This problem is addressed often using image registration of non-rigid body. However, this can be slow and hard to implement. On the other hand, rigid body registration is sometimes faster and easier to implement. The downside is that rigid body registration doesn’t usually take deformation into consideration. In this paper, two rigid body registration methods were explored, which are the BFP method and the PA method. Those results are later compared with the ICP method.

Index Terms

Rigid body registration; Tumor registration; BFP; PA.

Citation

Wang K, Kesavadas T (2016) Medical Image Registration and Visualization on Tumor Growth with Time Series. Comput Sci Eng 1(1): 1005.

ABBREVIATIONS

ICP: Iterative Closest Point; VTK: Virtual Tool Kit; MRS: Magnetic Resonance Spectroscopy; PET: Positron Emission Tomography; STL: Stereo Lithography; MRI: Magnetic Resonance Image; BFP: Best Fit Plane; PA: Principal Axes; FEA: Finite Element Analysis

INTRODUCTION

Image registration is a process of finding one common coordinate for different images. There are two common methods for image registration, one is deformable registration, and the other is rigid registration. Deformable registration has been used largely for tumor registration. Brock et al. used the FEA method for liver and tumor registration [1]. The downside is that it takes substantial processing time. Lamecker and Pennec [2] used a cost function to minimize the uncertainty of correspondences and unreliability of image intensity. This approach ignores the information around the lesion so it loses the importation information of the lesion. Jenkinson and Smith [3] implemented affine transformation (a linear mapping that can map points from one coordinate to another) of brain images.

The method tried the global optimization registration and the author claimed that it will be more likely to take less than 1 hour to run the program. Kaus et al. [4] , developed a surfacebased registration method and implemented on human organs. The author extracted control points and achieved a running time within a few seconds. However, manually selecting control points can be tedious. Mohamed [5] built a statistical model to constrain the registration.

Cuadra et al., employed Maxwell’s demons [6] registration algorithm with lesion growth model [7]. One of the drawbacks is that the seed location requires expert’s manual choice. Gendrinet et al. [8], investigated the feasibility of real-time organ motion monitoring with rigid body registration. Mang et al. [9], did an extensive study on the consistency of rigid body registration and suggested using normalized mutual information for image registration between time points. Brett et al. [10], combined both affine transformation and non-rigid registration in a single cost function and tried different threshold for the cost function.

In this paper, we explored two different rigid body registration approaches on tumor registration. One is BFP Registration and the other is PA Registration. Based on our knowledge, no one has ever implemented these two methods on tumor registration. The transformation matrix is found by finding the transformation of the BFP, the PA between the source image and the target image, respectively. Results are later compared with the best known approach, ICP.

In order to visualize the aligned tumors and compare them, VTK is used with different colors representing different time series of tumor. Doctors are able to compare between these tumors more easily.

The rest of paper is structured as follows: Section 3 introduces the method and the data we used for registration. Section 4 discusses and compares the results among the 3 methods. Section 5 lists some future work we plan to do.

MATERIALS AND METHODS

In this section, we provide a detailed description of the method we proposed and the data we used to generate the result.
Methods

BFP: Curve fitting is a usual technique engineers or scientists like to use when describing the behavior of the data from either simulation or experiment.

The idea we came up was to use BFP to fit the data, and then found the transformation matrix to transform the plane to the template plane. BFP belongs to polynomial (an expression of more than two algebraic terms, especially the sum of several terms that contain different powers of the same variable (s)) fit, in our case, it is a first order linear fit. Consider the plane as a set {Pi } where i goes from 1,…n, n is the number of data points. Pi is a 3*1 vector represented in 3D space as $\left [ X_{i}Y_{i}Z_{i} \right ]^{T}$ . The equation of a BFP is

z = ax + by +c (1)

The coefficient {a, b, c} of the equation can be computed with 3 non-linear points. Since we have more than 3 points, finding {a, b, c} essentially becomes an optimization problem. That is, to find a {a, b, c} that minimizes the error between the original points and the points projected onto the plane. The error between the original data set and the BFP is

$E^{2}=\sum ^{n}_{i=1}\left \| Z_{i-} \left ( ax_{i} +by_{i}+c\right )\right \|^{2}$ (2)

Denoting X,Y,Z as 3 vectors, $X=\left [ x_{1}x_{2}...x_{n} \right ]^{T}, Y=\left [ y_{1}y_{2}...y_{n} \right ]^{T}, Z=\left [ z_{1}z_{2}...z_{n} \right ]^{T}$ .In order to solve this problem, let’s rewrite the equation as

$\left [ X Y 1 \right ]\begin{bmatrix} a \\b \\ c \end{bmatrix}=Z$ (3)

The equation now becomes a matrix multiplication format, Ax = Z. a ,b, c is solved by taking the pseudo-inverse of A.

$X=\left ( A^{T}A \right )^{-1}A^{T}Z$ (4)

As long as AT is full rank, a, b, c can be uniquely found.

The transformation matrix calculation is based on [11]. We assume a point-wise correspondence between the 2 BFP. The reason is that the rotation due to alignment is reasonably close and it will not exceed 90° .

PA Registration: Moment of inertia is commonly used in physics or other engineering fields, referring to its ability of resistance to the change of its motion. In computer vision, an image moment is some weighted average of pixels’ intensities, or a function, as described in section 3.1.2.1. The orientation of an image can be extracted from the image moment, which is referred to as PA. After getting the PA for both source image and template image, the transformation that is used to match the source image’s PA to the target’s PA can be used to transform the sourced image to the target image.

Moment of inertia: Based on [12,13], the 3D ordinary moments, mijk , is defined as follows

$m_{ijk=\int_{-\infty }^{\infty }\int_{-\infty }^{\infty }\int_{-\infty }^{\infty }x^{i}y^{j}z^{k}f\left ( x,y,z \right )dxdydz }$ (5)

Where i, j, k are the integers and i+j+k is the order of the moment, f is a binary function. In most of laser scans or medical image scans, the 3D object is approximated by a tessellation of basic shapes, such as triangle or tetrahedral. Hence the discredited formulation is more useful. Let A be the set that contains all the triangles, $A_{1},A_{2},A_{3},....A_{N}$ . N is the number of element. f(x, y, z) is defined as

$f(x,y,z)=\left\{\begin{matrix} 1, (x,y,z)\epsilon A & \\ 0, elsewhere& \end{matrix}\right.$ (6)

Considering only the tetrahedral case, the volume integral can be rewritten as a surface integral of each single element.

$m_{ijk}=\sum_{t=1}^{N}\left ( m_{ijk} \right )_{t}=\sum_{t=1}^{N}\int \int _{A_{t}}x^{i}y^{j}[z\left ( x,y \right )]^{k}ds$ (7)

mijk is the ordinary moment for each individual element. Moment contains some information on image properties.

This is the definition of first and second order moment [14].

$m_{000}=m_{010}=m_{001}=0$ (Standardized position)

$m_{110}=m_{101}=m_{011}=0$ (Standardized orientation)

The moment we used for extracting rotation information is the secondary moment, which means i+j+k=2. It can be used to calculate principal axis, as we will see later.

PA: PA can be used to describe the orientation of an object. PA works well even when the shape of the object is not symmetric.

Denote $\left \{ u_{1},u_{2,}u_{3} \right \}$ as the eigenvector of moment of inertia. When the 3D object is centered at the origin, the PA is defined as the eigenvector of the inertia matrix [13,15].

$I=\begin{bmatrix} I_{xx} -I_{xy} -I_{xz} & & \\-I_{xy} I_{yy} -I_{yz} & & \\ -I_{zx} -I_{zy} -I_{zz} & & \end{bmatrix}$

where

$I_{xx}=m_{020}+m_{002},I_{yy}=m_{200}+m_{002},I_{zz}=m_{200}+m_{020}$ (8)

$I_{xy}=I_{yx}=m_{110},I_{xz}=I_{zx}=m_{101},I_{yz}=I_{zy}=m_{011}$ (9)

Ambiguity Elimination: Based on the previous section definition, { $u_{1},u_{2},u_{3}$ } are the eigenvectors (a vector which, when operated on by a given operator, gives a scalar multiple of itself) of I matrix. Assume that {u1 , u2 , u3 } is the normalized eigenvector, the sign of ui cannot be uniquely determined. To see why,

$Iu_{i}=\lambda u_{i},I\left ( -u_{i} \right )=\lambda \left ( -u_{i} \right )$ (10)

There are eight groups of eigenvectors in total,{±u1, ±u2, ±u3}.This suggests that there are eight orientations that can describe the 3D object.

The ambiguity in this paper was eliminated based on Galvez and Canton [14]. First, the requirement of right-handed system will eliminate four combinations of PA. Then a heuristic procedure was developed to get rid of the rest three combinations.

Transformation matrix extraction: After extracting PA from both images, we need to find a way to transform sensed image’s axes to template’s axes. Given two sets of axes $X_{1}$ and X2, suppose they are at the same origin, the rotation matrix that transforms 1 set of axis to the other is just a linear transformation:

$X_{1}=RX_{2}$ (11)

Where $R,X_{1},X_{2}$ are 3*3 matrices. Because $X_{1},X_{2}$ are principal axes, they are full rank and the inverse of either matrix exists, R matrix can be found by:

$X_{1},X_{2}^{-1}=R$ (12)

BFP: The following are the basic steps of BFP method.

Step 1: Load the STL file of the organ.

Step 2: Calculate the centroids of two data sets, p and p0.

Step 3: Find the BFP, more specifically, find 4 vertices used to form the BFP for two data sets.

Step 4: Calculate the rotation matrix R and translation matrix T.

Step 5: Transform the source data to target use the transformation found with the BFP.

PA: This is the steps for PA method.

Step 1: Import the STL file into a CAD software.

Step 2: Find the PA by using the values from the moment of inertia of the area at the centroid.

Step 3: Calculate the centroid of source image and target image sets.

Step 4: Translate all the tumor to the origin.

Step 5: Calculate the rotation matrix R with equation 12.

Step 6: Rotation the source image with rotation matrix R.

Tumor data representation

There are different modalities that can be used for collecting data of a tumor from a person, like CT scan, MRS and PET. The data we got was STL files on different time series, which belongs to OSF Saint Francis Medical Center. The data were initially taken from MRI scans. The doctors then did a 3D reconstruction using a medical image software called Osirix. This would give a 3D block of data, which contained lung, skeleton, bones, fat, heart, catheter and muscle. To extract a 3D organ or a tumor from the block, the doctor chose a threshold, which was used to separate different parts from the data. The data was saved as separate STL files.

The purpose of registration is to transform all other tumors from different time series to the first time series. In other words, the first time series is the target. Figure 1 shows the original data. Different time series of the tumor suggest that the tumor is shrinking.

Figure 1: Shows the original data. Different time series of the tumor suggest that the tumor is shrinking.

RESULTS AND DISCUSSION

Figure (2-4) shows PA Registration, BFP Registration and ICP, respectively.

Figure (2) shows the result of PA Registration performs better than BFP for translation. The rotation between two time series is different because the surrounding clusters do not match. This can be due to the reason that we sub sampled the tumor in order to feed into solid modeling software. The second moment of inertia is thus inaccurate. PA based method is popular because even the shape of the object is unknown, PA describe the original object well. Finding the transformation matrix is constant time so is getting the principal axis from a CAD software (as long as the number of vertices are small), so fast computation can be performed.

BFP is a simplistic way of finding the features of data. The underlying assumption is that the shape of the image from different time series does not change much. Ideally, if the segmentation from MRI image is good and the shape is well preserved, the BFP should be a good representation of the spatial location of the data throughout different time series. The run time for calculating the BFP is O (n). The run time for calculating transformation matrix is constant time because only 4 vertices are required for this calculation.

Figure 2: Shows the result of PA Registration performs better than BFP for translation

Figure (3) is the result for BFP Registration. Registration results show that the translation does not work well. It is hard to judge rotation because the shape of tumor changed. There are many parts in the tumor, while matching the core part of tumor is easy; matching all of the connected clusters is hard. As we can see, none of the above images match perfectly. The advantage of using BFP is its easy implementation and speed especially when shape of the two objects is close; it is a good representation of the data set features. However there are two main disadvantages. One is that when the tumor has significant deformation, the BFP may not represent the data well. The second is that when two hearts are upside down from each other, the one to one correspondence between 2 BFP changes. At this moment, there is no way to check whether the correspondence is correct.

Figure 3: Shows PA Registration, BFP Registration and ICP, respectively.

Figure (4) shows the result for ICP algorithm. ICP algorithm minimizes the total Euclidean distance between each point in the template and other time series. The result shows that each cluster of tumor is reasonably close to the template. However, ICP takes more than 4 hours to run with 20 iterations.

Figure 4: Shows PA Registration, BFP Registration and ICP, respectively.

CONCLUSION

In this paper, we studied three methods on tumor registration. Comparing the three methods, ICP seem to produce the best results. But ICP takes too much computing time. The workstation we used has Intel Core i7-5820K and nvidia geforce gtx 980 graphics card. With this configuration, PAR and BFT can be implemented a few minutes but for the case we studied, ICP took 4 hours to complete.

Several research challenges will be addressed in the future. First, instead of relying on a CAD software to calculate the second moment of inertia, we plan to calculate it internally through a new algorithm. It is worth comparing the results by taking into account all the points. Second, we will compare three approaches with a larger dataset from more patients. Third, we plan to test a new nonlinear registration method.

ACKNOWLEDGEMENTS

We would like to thank OSF Jump Simulation Center, specifically, Dr.Matthew Bramlet, Lela.G. Dimonte and Brent Cross generously provide us with all the image data we used for simulation.

JSM Computer Science and Engineering

JSM Computer Science and Engineering