Next: 3. Starting and Running Up: 1 A. Overview Previous: 1. Introduction Contents

Subsections

2. Background

2.1 Applications of Medical Imaging Analysis: A Brief Overview

BioImage Suite can be applied to help analyze images in a large number of domains. We briefly review some of the key application areas below.

Structural Neuroimaging: Measurement of brain structure is important for characterizing differences associated with normal variation (sex, handedness), development and aging, and pathology. Anatomic magnetic resonance images (MRI) are used to measure cortical and subcortical structure. Key problems in neuroimaging include the measurement of gray and white matter volume in the brain as well as segmentation and quantification of individual structures such as the hippocampus or caudate nucleus. Structural differences between groups of images (e.g. autism patients vs controls) can be quantified using non-rigid deformation based techniques (e.g. deformation based morphometry).

Diffusion Tensor Magnetic Resonance Imaging: Diffusion tensor MR imaging (DTI), allows the quantification of oriented tissue in terms of its diffusion properties as well as with fiber tracking. While typically applied to white matter in the central nervous system, DTI is also used for characterizing muscle (cardiac, skeletal) as well as other structures. For example, the structural changes associated with cancer can be measured with DTI. The key processing steps are diffusion tensor estimation at each location, diffusion property calculation and fiber tracking. In addition, multisubject registration and segmentation may also be needed for quantification.

Functional Neuroimaging: Functional magnetic resonance imaging (fMRI) measures blood flow changes associated with brain function. It has been widely used since its development in 1990 for localizing functional regions of the brain, characterizing normal variation and understanding pathological differences. fMRI is also used clinically as an aid for neurosurgical planning. fMRI activation detection and comparison involves a variety of processing steps including registration for motion and distortion correction and the formation of group composite functional maps for the analysis and comparison of subject groups. Cortical Blood flow measurements when combined with more conventional blood oxygenation dependent contrast provide a more quantitative measure of brain activity through calculation of cerebral metabolic rate of oxygen consumption (CMRO). Positron emission tomography (PET) requires many of these same processing steps, in addition to methods for the quantification of blood perfusion and receptor kinetics [15]. The latter requires extensive processing.

Cardiac Deformation Analysis: The quantification of cardiac deformation patterns from sequences of 3D MR and 3D echocardiographic images is a key step toward both understanding normal heart function and for evaluating the effects of acute coronary occlusions and post cardiac insult remodeling. The key image analysis methods here are segmentation and non-rigid deformation estimation.

Angiography: The detection and quantification of different parts of the vasculature has important applications in the understanding of many diseases such as cancer and diabetes.

Image Guided Neurosurgery: Image guided neurosurgery can leverage virtually the entire spectrum of medical imaging and image analysis techniques including segmentation, registration, fMRI and MRS analysis etc. In epilepsy neurosurgery, a particular focus application, localization of intra-cranial electrodes used to detect seizure sites is an important capability.

Image Guided Prostate Radiotherapy: The key image analysis problem in prostate radiotherapy is the localization of the prostate and nearby soft tissue structures such as the bladder and rectum via the registration of pre-therapy CT images and intra-therapy CT acquisitions. Such registrations may also require the use of articulated transformation models to handle the movement of different bony structures in the hip and pelvic regions. The registration can then be used to update the radiation plan (which is based on pre-therapy images).

Abdominal Fat Quantification: The quantification of intra-abdominal visceral fat is becoming a key step in the understanding of the relationship between obesity, insulin resistance and the development of Type II diabetes. These methods currently require image segmentation (in terms of voxel classification) as well as the use of magnetic resonance spectroscopic methods (MRS) for the quantification of intra-organ and intra-muscle fat which is dispersed in the tissue and not easily detected using magnetic resonance imaging methods.

Metabolic Imaging: Magnetic resonance spectroscopy (MRS) is used to quantify biochemicals (metabolites) in vivo by measuring the magnetic resonance chemical shift spectrum and quantifying the characteristic peaks corresponding to chemicals of interest. (a.) Proton MRS can be used to measure the concentrations of cerebral metabolites, such as NAA (N-acetyl aspartate) and creatine, and with improved MR methodology also glutamate, glutamine and GABA ( $\gamma$ -aminobutyric acid). Proton MRS is also applicable in the rest of the body for water/fat imaging which is important for the study of diabetes [122] (b.) POCE (Proton Observed Carbon Edited) MRS can be used to look at metabolic fluxes, like the tricarboxylic acid (TCA) cycle flux. (c.) $^{31}$ P CSI is a key technique in understanding energy metabolism in a wide range of disorders, including epilepsy. MRS is typically performed in a small number of voxels. Our goal is to implement efficient reconstruction methods that will make true spectroscopic chemical shift imaging (CSI) a practical modality. In addition, accurate segmentation (e.g. gray/white matter in the brain [25]) is essential for quantifying metabolite concentrations and registration methods are required to map the metabolic images to the underlying anatomy.

2.2 Medical Image Processing & Analysis

A large number of image processing and analysis methods has been developed over the last 20-30 years (a good but slightly dated critical review of the field can be found in Duncan & Ayache [29].) The following is a categorization of most of the key tasks in Medical Image Analysis. BioImage Suite aims to be a comprehensive image analysis software platform and does/will include representative algorithms in all of these categories.

Pre-processing: These initial image processing algorithms are often used in preparation for subsequent image analysis algorithms. They may also be used in their own right for the purpose of visualization. A good example is skull stripping [102] in neuroimaging which allows for the direct display of the cortex using volume rendering techniques. This visualization can be used for neurosurgical planning applications. Another example of this type of processing is MR bias inhomogeneity correction [100,8,115,111,43] for both brain and abdominal images which enables accurate tissue classification. We have developed innovative methods for bias field correction involving direct measurement of B1-fields [120].

Voxel Classification: These techniques label each voxel in the image according to tissue classes. Typical applications of these methods are the classification of brain tissue into gray matter, white matter and cerebrospinal fluid (CSF) [124,129,95,85,128,30] in neuroimaging and fat/lean muscle classification in abdominal imaging [80].

Deformable Surface Segmentation: Deformable surfaces either using explicit parameterizations [55,103] or levelset formulations [88,94,58,16,127] are applicable to all parts of the body as a method for measuring the shape and volume of different organs and sub-cortical brain structures [127]. In addition to automated methods, there are also applications of manual tracing and semi-automated methods in many practical applications both for complete segmentation or for the review/correction by an expert user of the results of automated segmentation such as in the BioImage Suite surface editor tool [79,73,12].

Rigid and Articulated Registration: Rigid registration allows for the alignment of different images of the same patient in a common coordinate space [10,89,109,126], that subsequently permits combined analysis and visualization of the data. Such methods can be used both for intra-modality registration (e.g. serial MR, serial CT) or for multimodality registration (e.g. MR/CT, MR/PET etc.) An extension of these methods is our recent articulated piecewise rigid work for mapping serial lower limb mouse CT images [75] while accounting for differences in joint angles.

Non-rigid Registration: Non-rigid registration methods generate transformations that allow the mapping of images of different patients to a common coordinate space [17,35,113,91,114,20,24,22]. In addition, such methods are also used to map post-surgical images to pre-surgical images of the same patient [110]. Once images from different patients are mapped into a common coordinate space, they can be used to generate voxel-wise group statistics for the study of within-group variability (both structure and function) as well as the quantification of statistical between-group differences (e.g. normal controls vs patients).

Diffusion Tensor MR Imaging (DTI): This novel imaging technique allows for the quantification and characterization of local brain tissue diffusion properties using a tensor model [9] and the construction of potential white matter axonal pathways connecting different brain regions (e.g. [54,67,60,50,51]).

Vessel Extraction and Tracking: The localization and quantification of vascular structure has important applications in the evaluation of peripheral and coronary artery disease, as well, as in the study of vascular growth and remodeling [31].

Cardiac Strain Analysis: The quantification of cardiac motion and strain is crucial for the understanding of cardiac function, injury and remodeling [23,63]. Image analysis in conjunction with biomechanical modeling [42,45] allows for the measurement of regional strain from cardiac cine MRI [5,97] and echocardiography[82].

fMRI Processing: Functional magnetic resonance imaging analysis involves a series of processing steps including: (i) motion correction of the timeseries raw T2-weighted images (e.g. [40,39,36,66], (ii) smoothing and denoising (iii) detection of activated regions [37,84,38,65,121,99,69,21] and (iv) distortion correction for mapping the functional images to anatomical MRI images [53,106].

Magnetic Resonance Spectroscopy and Spectroscopic Image Analysis: Typical MRS spectra are quantified by first approximating the spectrum with a linear combination of in vitro spectra (after shifting and broadening) measured with the exact same sequence: this is known as the LC-model [87]. Typical implementations of the LC-model, such as the one by Dr.<22968>>Robin de Graaf of the Yale MRRC, require approximately 10 minutes for each voxel. While the computational time is feasible if each study consists of only a handful of voxels, the problem becomes intractable for larger studies such as a technically feasible CSI acquisition consisting of 32 $\times$ 32 $\times$ 4 = 4096 voxels. Current acquisitions are commonly restricted to low resolution measurements [18]. While simpler techniques can be used the full power of CSI will not be realized until the same proven procedures that have been applied at the voxel level become available at the image level.

Compartmental Modeling: PET and MR Perfusion The dynamic images that PET provides can be used to model the PET tracer kinetics to produce quantitative estimates of physiological parameters. Estimation of these parameters, including blood flow, metabolic rates and receptor concentrations, requires weighted nonlinear least squares fits of the time-activity curve of each voxel to the appropriate compartmental model equations. These same techniques can be adapted for dynamic contrast enhanced magnetic resonance images (DCE-MRI) and MR perfusion imaging.

2.3 Software Development Related to Medical Image Analysis

Please Note: This section represents work in progress and is by no means complete.

In this section we review related work in software grouped into two categories (i) development libraries for numerical computation and visualization, as well as, libraries specific to image analysis and (ii) complete software packages (both commercial and free) that are currently available to support biomedical research.

a. General Development Libraries: Medical image analysis in general, and medical image analysis software in particular, rely on many core algorithmic techniques from a number of fields such as numerical analysis, graphical user interface design and computer graphics and visualization. Key examples include:

Numerical Processing: Much of the available code for numerical processing was originally developed in Fortran 77 often converted to C using the freely available f2c converter. A good example is the LAPACK linear algebra library [32] which is also available for C as CLAPACK. Other freely available libraries include the GNU Scientific Library [33], FITPACK (a spline fitting toolbox [26]), and the vnl:numerics library.

Computer Graphics and Visualization: The cornerstone of most scientific 3D Graphics and Visualization is the Open GL library [96,70] originally developed by Silicon Graphics. Open GL is a low level library and is often used within a higher level toolkit such as the Open Inventor toolkit [125], or the Visualization Toolkit (VTK) [92]. VTK has gained substantial popularity for medical image analysis as it also provides common image/surface processing algorithms, and functionality for some numerical processing in addition to a complete set of graphics routines. BioImage Suite relies on VTK for most of its visualization tasks.

Graphical User Interfaces: There is a large number of possibilities for graphical user interface design, although the list shrinks dramatically in the case of cross-platform (i.e. Windows/Unix/Mac OS) support. Possible choices include WxWindows, FLTK, and the TK toolkit which comes together with the Tcl scripting language. BioImage Suite uses the Tcl/Tk combination together with a custom C++ wrapper around Tcl to enable the creation of user interfaces in both a scripting language (Tcl) and C++ (as needed for complex components).

Medical Image Analysis Libraries The key development in this area over the last 3-4 years has been the development of the NIH/NLM Insight Toolkit (ITK). ITK focuses exclusively on medical image segmentation and registration, and provides implementation of many commonly used algorithms. BioImage Suite uses some aspects of ITK (e.g. the levelset libraries) as a basis for implementation of some methods and will make accessible from the user interface some of the algorithms implemented in ITK that complement what is currently available.

An alternative is the Matlab programming environment (MathWorks, Inc. [62]) which provides a programming language, a large number of numerical methods as well as basic graphical interface capabilities. Some complex software packages are developed entirely in Matlab (e.g. SPM). However, while Matlab is an excellent prototyping tool, algorithms that do not naturally fall into a matrix manipulation paradigm can be highly inefficient.

2.3.0.0.1 b. Complete Software Packages

There are a number of software packages available for medical image analysis that provide some overlapping capabilities; some are commercial but most are research-based. While a number of these packages offer strong features in specific areas, there is currently no software package that provides state-of-the-art methods in the array of techniques that are included in BioImage Suite.

The neuroimaging packages available are probably the most well developed. The Slicer [118] package has a similar design to BioImage Suite (built upon VTK/Tcl) and provides many valuable methods. Slicer focuses on registration and segmentation techniques geared for surgical planning, particularly in the brain. Brainsuite [2] focuses primarily on segmentation and visualization of anatomical brain MR images. SPM [104] is a Matlab package for fMRI and PET activation analysis. FSL [41] also focuses on image analysis and statistical tools for brain imaging, including fMRI and DTI. AFNI [3] was originally specific to fMRI image analysis, although has been recently expanded for other tasks such as DTI analysis. Freesurfer [112] is a set of command-line neuroimaging utilities for the semi-automated reconstruction of the cortical surface and overlay of functional data onto the reconstructed surface. Map3D and BioPSE [11] are packages primarily for bioelectric field modeling program but they have many useful visualization and analysis methods including FEM analysis, DTI analysis and basic image processing capabilities. AIR [4] is limited to image registration. NIHImage/Scion [47] provides primarily basic image processing functions.

Some commercial packages are available, both general purpose and application focused. Commercial packages can be expensive, may be tied to a specific niche market and their development cycles are often too slow for researchers to influence. They also typically do not have the benefits of open source. MEdX [64] is software for visualization, processing, and analysis. It primarily serves as an interface to functionality integrated from AIR, SPM and FSL. VIDA [117] is software primarily for 3D visualization with an emphasis on cardiac applications including image segmentation and cardiac mechanics. The general purpose commercial software package Analyze [7] has strengths in visualization, segmentation, registration and measurement. BrainVoyager [119] is, like SPM, specific to brain activation analysis.

2.4 3D Graphics and Volume Rendering

2.4.1 3D Graphics

Most researchers' first experience with image display consists of simply displaying a two-dimensional image slice (most likely in Matlab). While there are some issues, such as the relative orientation (row/column vs column/row) and the position of the origin, it is for the most part a fairly straightforward and intuitive procedure. The same applies to rendering curves and point landmarks on the image. The key behind this simplicity is the fact that, unsurprisingly, two-dimensional computer monitors are well suited to displaying two-dimensional content!

The move to 3D graphics and visualization requires, in some way, getting around the fact that for the most part our display units are two-dimensional. The most common ``illusion'' used in 3D graphics - and the one that is employed by the Visualization Toolkit (VTK) and hence by BioImage Suite - is that the program first generates a 3D world/scene consisting of three-dimensional entities/models - actors/props in VTK parlance - which have various appearance properties. The world is illuminated by a set of lights and the ``user'' looks at the world through the eyes of a virtual camera - See Figure 2.1 for an example. The image that gets displayed on the computer monitor is precisely the output of such a virtual camera. To recap, a world consists of (i) Actors, (ii) Lights and (iii) a camera. The exact output naturally depends critically, in addition to the Actors themselves, on the position and orientation of the camera. If the camera is looking in the wrong direction, even if all the actors have been generated by your program, it will not show anything on your computer monitor. Just as in the real world, when you have to take a picture of a scene, you point the camera towards that scene, in this 3D-world the camera has to be position correctly located (far enough from the object) and directed towards the scene.

Figure 2.1 shows a schematic of a camera and its properties. The frustum defined by the near and far planes forms the ``view volume. '' Everything inside the view volume is shown on the screen whereas everything outside the view volume is ``clipped'' out from the scene and is not shown on the screen. The projection of the scene/actors on the view plane is shown on the screen. The view plane can be the same as the near plane.

**Figure 2.1:** The Camera Model for 3D Graphics. A camera looks at the synthetic world (e.g. the bunny) and takes a picture of it - this is what is displayed on the screen/view plane. An additional advantage of this virtual world over the real world is that we can restrict the viewing volume to exclude objects that are ``too close'' or ``too far'' - hence visually cropping the scene.

The up-vector defines the ``up'' direction for the camera. Based on the up-vector, the image shown on the screen will be shown. Figure 2.2 shows a schematic of the result of having a up vector and an inverted up-vector. The inverted up-vector causes the projected image to be inverted. As can be seen, the camera in 3D graphics is very similar to a camera in the real world.

**Figure 2.2:** The Up-vector for the Camera. The left image shows the result of the camera vector pointing upwards. The right image shows the result of the up-vector pointing downwards. The projected image on the view plane is inverted.

All display operations in BioImage Suite are done using 3D Rendering, even when a single image slice is involved. This allows, for example, for fast zooming in and out as this operation simply involves changing the zoom on the ``virtual'' camera as opposed to redrawing the image bitmap scaled. This type of renderer also allows the seamless mixing of images and polygonal data (e.g. surfaces, curves) in the display windows.

2.4.2 Displaying 3D Images

**Figure 2.3:** An example of an image (top) and how a computer looks at it (bottom).

Most people begin to think of images as essentially a matrix of numbers - see for example Figure 2.3. The value of these numbers may be arbitrary (e.g. anatomical MR acquisitions) or have physical meaning (e.g. a fractional anisotropy map). To display these image we somehow must convert each number to a color. The only complication, at this stage, is the relationship between the row/column indices of the matrix and the i- and j-axes of the image, respectively.^2.1.

Notation: We will use i,j,k to refer to the internal axis of the image in voxels and x,y,z to be the coordinates in mm.

This type of radical simplicity disappears, unfortunately, once we have issues with mapping images. Examples include, image registration applications when we are looking to estimate a transformation between two images, statistical shape model building etc. In general, in medical imaging, we need to keep track not only the image intensities but additional attributes such as voxel dimensions (voxels need not be isotropic), image orientation (the relationship of the i-,j- and k-axis orientation to the human body - see Figure 2.4), the position of the voxels etc.

**Figure 2.4:** Standard Image Acquisition Orientations, Axial or Transverse, Coronal and Sagittal. The arrows indicate the z-axis direction (which may be inverted depending on the acquisition protocol), the x-axis and the y-axis are perpendicular to this. **Right:** Axial, coronal and sagittal slices. Oblique acquisitions are also sometimes used, in which the z-axis is rotated, or *obliqued* away from one of the standard acquisitions, e.g. coronal oblique.

2.4.3 Volume Rendering

Volume rendering is defined as the process of generating a 2D image directly from three-dimensional data. The basic idea is to simulate the absorption and emission of light as it goes through matter.

To simulate the passage of light through data, rays are cast from the image plane into the space containing the volume, as shown in Figure 2.5. For every pixel in the image plane, a ray is cast into the volume space that traverses the volume. At each point along the ray, the data is sampled to identify the density/intensity at that point. This value is different for different materials such as bone, tissue, fat and so on. Based on the intensity at that point along the ray, the color and opacity is identified using a lookup table called a colormap or transfer function.

2.4.4 Types of compositing functions

Color and opacity are accumulated along the ray and various compositing operators are used to get different effects. An X-ray image can be simulated by averaging the intensity values sampled along the ray. A MIP (Maximum intensity projection) image is obtained by preserving the maximum intensity value along the ray. MIP images are frequently used to visualize vascular structures as they can be clearly seen using this technique. Figure 2.6 shows an example of the two techniques.

**Figure 2.5:** This image depicts a schematic of the volume rendering process. For every pixel in the image plane, a ray is cast into the volume space. Along the ray, the volume rendering integral is evaluated to simulate the passage of light through matter. At each position along the ray, the volume/data is sampled to identify the value at that location. The location identifies the material properties such as absorption and emission. A color and opacity is obtained for every density by performing a lookup into a table, called the *colormap or transfer function*. This allows users to color different regions of the volume differently.

The compositing technique used for direct volume rendering is based on accumulating the color based on the color and opacity of the current voxel under consideration. For example, during ray traversal if in a CT scan, bone is encountered, it absorbs more energy than tissue or a thin vessel.

**Figure 2.6:** The left image shows an X-ray style image that can be generated by average the intensity values along the ray. The right image shows a MIP image that can be obtained by picking the maximum intensity observed along that ray. Such MIP images are widely used to visualize vessels.

2.4.5 Types of Volume Rendering

There are four types of volume rendering techniques: raycasting, texture mapping-based, shear warp and splatting. Raycasting and texture mapping-based techniques are the most widely used of the four. Raycasting is the technique in which an image is generated by casting rays into the volume space.

As graphics hardware got better with time, researchers devised a way to perform volume rendering using graphics hardware. This technique loads the data into graphics memory and utilizes 2D (and more recently 3D) texture maps to perform interpolation and blending. Since graphics hardware is extremely fast, this technique is much faster than the raycasting technique described earlier.

2D texture mapping based volume rendering

2D texture mapping based volume rendering uses the ability of graphics hardware for fast rendering of polygons and textures. Figure 2.7 shows a schematic of the process. The leftmost image shows polygons that are aligned to an axis. These polygons, which are collectively known as the proxy geometry, are texture-mapped based on the values in the data. A blending step is performed that blends the result of each texture-mapped polygon in a back-to-front manner to get a volume rendered image, as shown in the rightmost image in Figure 2.7.

**Figure 2.7:** This figure shows a schematic of the axis-aligned polygons on the left, that are texture-mapped and blended to get a volume rendered image of the data.

There are a few disadvantages with using such 2D-textures though. For this type of volume rendering, three copies of the data have to be stored in axis-aligned slices (one for each X, Y, Z-axes), as shown in Figure 2.8. Based on the view vector of the camera, the sets of data that is most perpendicular to the view vector is selected. The biggest drawback of using three sets of axis-aligned 2D textures is that as the viewer rotates the camera around the data, there is a visual popping artifact, when the algorithm switches from one set of axis-aligned slices to another.

**Figure 2.8:** Three axis-aligned copies of the data have to be stored in memory to allow 2D texture mapped volume rendering. The XY-aligned slices are used when the view vector is aligned to the Z-axis, the YZ-aligned slices are used when the view vector is aligned to the X-axis and similarly, the XZ-aligned slices are used for the Y-axis aligned view vector.

3D texture mapping based volume rendering

The advent of 3D textures helped solve the problems of maintaining 3 copies of the same data in memory as well the visual popping artifacts. The data can now be stored in a single 3D texture and instead of using axis-aligned slices, a new algorithm was developed to draw ``view-aligned'' slices, as shown in Figure 2.9. View-aligned slices ensured that slices were drawn orthogonal to the view vector at all times and ensured that there would be no visual artifacts.

**Figure 2.9:** View aligned slices eliminate the need for storing three sets of data. The view aligned slices are interactively drawn perpendicular to the view vector, as the user rotates the camera around the data.

Figure 2.10 shows a schematic of the 3D texture mapping-based volume rendering approach. The render geometry consists of polygons that are drawn perpendicular to the view vector as the viewer rotates around the data. The rendering of proxy geometry, texture mapping the data from the 3D texture onto it and the blending, all happen at run time and provide a volume rendered image of the data to the viewer.

**Figure 2.10:** A schematic of the 3D textures based volume rendering process. The view aligned slices are textured using the 3D textures and blended to provide a volume rendered representation of the 3D data.

Next: 3. Starting and Running Up: 1 A. Overview Previous: 1. Introduction Contents