Transcript DissertationProposal

PRESENTATION
MRI-BASED ATTENUATION
CORRECTION FOR PET/MRI AND
MEDICAL IMAGE REGISTRATION
METHODS ON GPUs
Dissertation Proposal
by
Wei-Hung Cheng
September, 9, 2009
Introduction
Image acquisition systems:
– Magnetic Resonance Imaging (MRI)
– X-ray Computed Tomography (CT)
– Ultrasound
– X-ray mammography
Nuclear medicine imaging techniques:
– Positron Emission Tomography (PET)
– Single Photon Emission Computed Tomography (SPECT)
– cardiovascular imaging
– bone scanning
Medical image analysis methodologies :
– image segmentation
– image registration
– motion tracking and change detection from image sequences
– the measurement of anatomical and physiological parameters from images.
Medical image detector research on multimodality imaging:
– PET/CT methodology
– SPECT/CT methodology
– PET/MRI methodology
Our research focuses on:
– MRI-based attenuation correction on PET/MRI methodology
– The CUDA programming model with a GPUs platform for medical image
registration.
Background
PET/MRI Methodology
Attenuation Correction (AC) Processes
Segmentation Approach
Registration Approach
Combined Approach
Image Registration on GPUs Platforms
CUDA Programming Model
PET/MRI Methodology
The status of estimating PET
attenuation maps to form MRI
Attenuation correction
– attenuation maps: assign to
the volume of the respective
attenuation coefficients
Standard method for finding the
attenuation map is to perform
transmission scans with CT
measurements
MRI (AC) → CT(AC) → PET
(A)Non-attenuation corrected images and (B)
attenuation corrected images in a patient with liver
metastasis
Attenuation Correction (AC) Processes
PET AC processes in 2 stages:
1 – Finding the attenuation map
2 – Applies the attenuation map to PET emission data
The sinogram SINac is corrected for photon attenuation which can be
determined by point wise multiplication of every entry in ACF with the
corresponding entry in the sinogram SINdet which are then measured by the
PET detectors.
Equation:
SINac = SINdet × ACF
The AC factor (ACF) of a line of response (LOR) is the measured signal I
to the signal I0 which can be measure in the absence of attenuation. µPET is
the linear attenuation coefficient for E = 511-keV photons.
Equation:
ACFLOR = I0 / I = exp [ ∫LOR µPET × dx]
Segmentation Approach
Locating and mapping the primary attenuation structures in the scan object.
In general, it can be achieved in two stages:
Segmentation stage – segments the scan object into area of tissues and
organs which have different attenuation properties.
Mapping stage – maps the segmented tissues and organs to corresponding
linear attenuation coefficients at 511 keV.
PET Image
Registration
Transformation
from MRI to
PET image
Improved
Reconstruction
Algorithm
MR Image
Generation of
Attenuation Map
Attenuation Map
Darko Zikic proposed :
Pre-processing the image to apply mutual information algorithm (MI) for
registration of the MR and PET image
– The MI algorithm failed because some features in MR images are
missing in the PET images such as all tissue does not show activity.
– Pre-processing method: an intensity based seeded region growing
algorithm to track the segment region containing all voxels and a
simple binary threshold filter to removing background noise.
– Several alternative registration algorithms : orthogonal projections
algorithm, line matching algorithm, and marker-based registration.
The generation of the attenuation map is to distinguish different tissue
types, such as bone, bone marrow, muscle, fat, lungs, and air present in the
MR images.
Two approaches for the generation of AC such as intensity based
segmentation and joint histogram clustering
Intensity based segmentation:
– iso-density algorithm : segments the intensity histogram into the number
of classes (air, muscle, fat, and bone-marrow).
– a multiple pass binary threshold filter algorithm: gives the user more
control of the process with more flexibility.
Exemplary results of the iso-density segmentation algorithm on MR
images of a human foot: Top: PD sequence, 4 classes; Middle: T1
weighted, 4 classes; Bottom: T2 weighted, 6 classes
Joint histogram clustering: combines two different MR sequences’ intensity
values into one joint histogram.
Four parts of algorithm:
Co-registration – two different MR images of the same scan object
Generation of the Joint Histogram – same clustering algorithm
Clustering of the Joint Histogram – K-means, fuzzy means, mixture
of Gaussians, self-organizing feature maps, and hierarchical
clustering.
The output image production – after the result of clustering the joint
histogram, cluster values are assigned to the output image.
Results of the joint histogram clustering algorithm after manual
adjustment. (PD and T1 weighted MR images)
Registration Approach
MRI/CT atlas registration was proposed by Kops, Herzog, and Hofmann.
It used MRI with atlas registration to determine the attenuation map to
predict a pseudo-CT.
Required an atlas MRI image and a corresponding CT image on the same
subject.
First Step – aligning MRI image with the new subject’s MRI image by a
nonrigid registration algorithm which can be computed automatically for
deformation.
Second Step – applying the deformation from the first step to the atlas CT
image that could get the desired result.
Principle of atlas-based MR-AC. The atlas consists of a matching MR-CT image volume that can be
generated from a patient. An atlas of attenuation values at 511 keV is generated from matching CT images.
The atlas MR image (top left) is coregistered to the MR image volume of a specific patient (bottom left).
This transformation is applied to the corresponding CT atlas, thus generating an attenuation image (i.e.
pseudo-CT image) that approximately matches the patient anatomy
Combined Approach
Hofmann and Steinke proposed combining local pattern recognition with
atlas registration.
First – they used pattern recognition with Gaussian process by adding the
registered coordinate of a training point as input and transform two
different images into a single coordinate system.
Second – they included prior knowledge of the atlas registration by setting
the mean function to the average value of the registered CT images.
Atlas Database of
corresponding MRCT Pairs
MRi
MR Patient
Image MRn
For every database MRi
calculate registration to the
patient MRn and apply to
MRi and CTi
CTi
PET
Detector
Sinogram
SINdet
MRregi
CTregi
Find all neighboring MR
patches in the registered
database
Perform Gaussian Process
Regression on patch and
position
This part is on a PET/CT scanner
SINac =
SINdet / ACF
Filtered
BackProjection
Attenuation
Corrected PET
Image
Overview of steps involved in Hofmann’s method for obtaining attenuation corrected PET image, based on
PET detector sinogram and MR image. Re-sampling of MRI, PET, and CT to required resolution is
performed wherever necessary.
Image Registration on GPUs Platforms
The computing capacities of
graphics processing units (GPUs)
have improved exponentially in
the recent decade.
NVIDIA released a CUDA
programming model for GPUs.
The CUDA programming
environment applies the parallel
processing capabilities of the
GPUs to medical image
processing research.
CUDA Programming Model
A parallel programming model and software environment designed to
handle parallel computing tasks.
Similar to the traditional single instruction, multiple data (SIMD) parallel
model.
Major abstractions:
– a hierarchy of thread groups
– shared memories
– barrier synchronization
It provide a programming model for data parallelism, thread parallelism,
and task parallelism.
Computation paradigm of the CUDA
A program is divided into blocks.
– A block is a group of threads mapped to a
single multiprocessor by the programmer to
share the memory.
Host
Device
Grid 1
Kernel
1
The data is also divided amongst all threads
in a SIMD fashion by the programmer.
Block
(0, 0)
Block
(1, 0)
Block
(0, 1)
Block
(1, 1)
Grid 2
Kernel
2
All threads are organized into warps.
– Each warp is a group of 32 parallel scale
threads, which can run concurrently on the
multi-processors.
Block (1, 1)
(0,0,1) (1,0,1) (2,0,1) (3,0,1)
Thread Thread Thread Thread
(0,0,0) (1,0,0) (2,0,0) (3,0,0)
Thread Thread Thread Thread
(0,1,0) (1,1,0) (2,1,0) (3,1,0)
Courtesy: NDVIA
Collections of warps are known as thread
block
Figure 3.2. An Example of CUDA Thread Organ
Memory hierarchy is in the form of
registers, constant memory, global memory,
and textures.
– registers: fastest level in the hierarchy,
a limited amount of space.
– constant memory: a subset of device
memory, cannot be modified at
run-time by a device.
– global memory: permits read and
write operation from all threads,
but is uncached and has long
latencies.
– textures memory: a subset of the
device memory, read-only on the
device, faster cached reads, allows
addressing through a specialized
texture unit.
Proposal PET/MRI Framework
A novel framework for MRI-based attenuation
correction for PET/MRI:
First – the registration of two different MR Images
(T1-weighted and T2-weighted).
Second – the co-registration of MR and CT
performed in conjunction of Mutual Information
(MI) based on gradient vector flow intensities
(GVFI).
Third – matching pair of MR and CT image
intensities and create corresponding look-up
table to map MRI intensities to pseudo-CT
values.
Last – derive a pseudo-CT image set by applying
the look-up table to the MR image set.
GVFI-MI
GVFI-MI can incorporate spatial information into mutual information to
increase the quality of the image with a higher success rate than traditional
MI.
Applying MI to gradient images would seem a logical solution to
incorporating spatial information.
Computed from the gradient of edge map, the GVF field not only keeps
large magnitude near the edge, but also extends into homogeneous regions
in the computational diffusion process.
Its property can help to improve the capture range of image gradient in the
registration procedure.
The GVF-intensity (GVFI) map of the image as
GVFI(x, y) = I(x, y) + GVF(x, y)
I(x, y) – given an image
GVF(x, y) – the magnitude of the vector from the GVF field
GVF(x, y) can be found by solving the following two Euler equations:
µ 2 u – ( f x2+ f y2 ) (u – f x ) = 0
µ  2 v – ( f x2 + f y2 ) (v – f y ) = 0
GVF(x, y) is derived from the vector field v which present the information
from GVF field as
GVF(x, y) = u 2  v 2
The algorithm to compute GVFI is as follows:
1. Compute edge map f and its gradients, fx and fy
2. Compute SMF (squared magnitude of the gradient field)
3. Initialize u as fx, and v as fy
4. While iteration < N do
5.
Compute Laplacians of u and v (noted as Lu and Lv)
6.
Update u: ui+1 = µ × 4 × Lu – SMF × (ui – fx)
7.
Update v: vi+1 = µ × 4 × Lv – SMF × (vi – fy)
8. End while
9. Compute GVF(x, y) = u 2  v 2 and GVFI(x, y) = I(x, y) + GVF(x, y)
The algorithm has three parameters to be determined: regularization
parameter µ, number of iterations N, and edge map method.
Determine the value of µ and N experimentally.
Edge map f(x, y) is derived from the image I(x, y).
Definitions of four edge map methods as
f1 = | Ix | + | Iy|
f2 = | Ix + Iy|
f3 = | (Gσ I)x | + | (Gσ I)y|
f4 = | (Gσ I)x |2 + | (Gσ I)y|2
Our MI-based multi-modality image
registration:
MI is computed on the GVFI
features instead of the intensities of
the original images.
Histogramming is used to compute
2D joint distribution of GVFI.
Bilinear interpolation is employed
throughout the registration procedure
Downhill simplex method is
selected as the optimization strategy.
Next steps
MR and CT image intensity transformation can be performed by nonproprietary histogram matching algorithm.
Based on image intensity match to create corresponding look-up table to
map MRI intensities to pseudo-CT values.
Applying the look-up table to the MR image to derive a pseudo-CT
image set
CUDA Algorithm for NMI registration
The registration procedure spent 90-95% of its run-time on mutual
information computation.
Medical image registration has a high level of data parallelism and image
data can be mapped onto the GPUs platform.
Based on our normalized mutual information method, the tasks are
classified into four CUDA kernels as follows:
Transformation – This group performs coordinate transform, affine
transform, and mapping matrix to establish spatial correspondence
between two images.
Interpolation – This group involves iteratively transforming image A with
respect to image B while optimizing the MI measure which is
calculated from corresponding voxel values.
Histogram – This group computes a joint histogram of the pairs of images
to evaluate the mutual information.
Optimization – This group detects optimization of estimate transformation
to evaluate its similarity
The CUDA implementation consists of four stages:
1) Allocate data memory on the device and transfer them from the host to
the device.
2) Set up the function kernel configuration.
3) Launch function kernel(s) and store the result in the device memory.
4) Transfer data from the device memory to the host memory.
GVF-based MI Experimental Results
We evaluated the performance of
GVF-MI compared to intensitybased MI, and test the robustness
of GVF-based MI with regard to
noise.
In the experiments:
– µ is set as 0.2.
– N is set as 80.
– The standard deviation σ for
Gaussian function Gσ in f3 and f4 is
set as 3.5.
(a) T1 image
(c) GVFI of (a)
(b) T2 image
(d) GVFI of (b)
T1 and T2 image pair and their GVFI
The robustness experiment: four pairs of
slices extracted from T1 and T2 images
were used, each at difference noise
levels, 3%, 5%, 7%, and 9%.
In the experiments:
– T1 slice was taken as reference image,
while T2 slice as floating image.
– We observed that the curves are
smooth at each noise level, within the
rotation rang of [-400, 400] about xy-axis,
and translation range of [-40, 40] on xaxis.
Success rates of registration methods based
on MI and GVFI-based MI
Dataset
Success Rate
MI
GVFI-MI
0% noise
85.3%
90%
3% noise
86.7%
90.7%
5% noise
85.7%
90%
7% noise
84.3%
93.7%
9% noise
85.3%
89.3%
CUDA Implementation Experimental
Results
The experiments involved the data sets of 7 patients, each consisting of
Computed Tomography (CT) and six Magnetic-Resonance (MR) volumes.
On a PC, having a 2.40 GHz Intel® Core™ 2 Quad CPUs, and 4 GB
DDR2 memory with NVIDIA’s GeForce 9600 GT graphic card.
All CT images were registered to the MR images using the MR image as
the reference image on PC
Run the registration procedure on both for the CPU-base platform (C
program), and the GPUs platform (the CUDA program).
Experimental results showed that the GPU implementation improves the
registration computational performance with a speedup factor of 23.4×
Comparison of GPU and CPU-based implementation for the
registration procedure runtimes
Run Time Range for 41 pairs data set on CPU and GPU
Run Time
Average
(mins.)
(mins.)
CPU-based
7.41 ~ 18.34
12.20
1
GPU-based
0.33 ~ 0.633
0.5
23.4
Architecture
Speedup
Conclusion
Based on the preliminary result on the GVF-MI, it shows promise for our
novel PET/MRI framework.
Create corresponding look-up table to map MRI intensities to pseudo-CT
values.
CUDA shows its potential as a medical image processing tool based on
the preliminary result.
We can use the CUDA on our novel PET/MRI framework
implementation to speed up the image processing procedure.
References
Kinahan PE, Townsend DW, Beyer T, Sashin D. Attenuation correction for a combined 3D PET/CT scanner. Med Phys.
1998; 25:2046–2054.
Chow PL, Rannou FR, Chatziioannou AF. Attenuation correction for small animal PET tomographs. Phys Med Biol. 2005;
50:1837–1850.
Diego Vivancos Gallego. Development of Pre-Processing Methods for Attenuation Correction Algorithms in Small Animal
PET/MR Imaging. Diploma Thesis, Techinsche Universität München, 2004.
Matthias Hofmann, Bernd Pichler, Bernhard Schölkopf, Thomas Beyer. Towards quantitative PET/MRI: a review of MRbased attenuation correction techniques. Europe Journal of Nuclear Medicine Mol. Imaging. 36(Suppl 1):S93 - S104, 2009.
Kops RE, Herzog H. Alternative methods for attenuation correction for PET images in MR-PET scanners. Nuclear Science
Symposium Conference Record, 2007. NSS’07. IEEE. Vol. 6. p. 4327–30. doi:10.1109/NSSMIC.2007.4437073.
Hofmann M, Steinke F, Scheel V, Charpiat G, Farquhar J, Aschoff P, et al. MRI-based attenuation correction for PET/MRI:
a novel approach combining pattern recognition and atlas registration. J Nucl Med 2008; 49:1875–83.
NVIDIA. Nvidia cuda compute unified device architecture. Programming Guide, Version 2.0. NVIDIA, 2008.
Yujun Guo. Medical image registration and application to atlas-based segmentation. PhD dissertation, 2007.
C. Xu and J. L. Prince. Snakes, shapes, and gradient vector flow. IEEE Trans. On Image Process, 7(3):359-369, Mar. 1998.
W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling. Numerical Recipes in C. Cambridge University Press,
1988.