1.

Introduction

Near-infrared (NIR) diffuse optical tomography (DOT) relies on functional processes, and provides unique measurable parameters with potential to enhance breast tumor detection sensitivity and specificity. For example, several groups have demonstrated the feasibility of breast tumor characterization based on total hemoglobin concentration, blood oxygen saturation, water and lipid concentration and scattering. ^{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17}

The functional information derived with DOT is complementary to structural and functional information available to conventional imaging modalities such as magnetic resonance imaging (MRI), x-ray mammography, and ultrasound. Thus the combination of functional data from DOT with structural/anatomical data from other imaging modalities holds potential for enhancing tumor detection sensitivity and specificity. To achieve this goal of data fusion, two general approaches can be employed. The first, concurrent imaging, physically integrates the DOT system into the conventional imaging instrument. This approach derives images in the same geometry and at the same time. The second approach, nonconcurrent imaging, employs optimized stand-alone DOT devices to produce 3-D images that must then be combined with those of the conventional imaging modalities via software techniques. In this case, the images are obtained at different times and often in different geometries.

Thus far a few DOT systems have been physically integrated into conventional imaging modalities such as MRI,^{18, 19, 20, 21, 22} x-ray mammography,⁵ and ultrasound⁹ for concurrent measurements. By doing so, however, these DOT systems are often limited by the requirements of the “other” imaging modality, for example, restrictions on metallic instrumentation for MRI, hard breast compression for x-ray mammography, limited optode combinations for ultrasound (and MRI, x-ray), and time constraints. On the other hand, among the stand-alone DOT systems available today, only a few attempts have been made to quantitatively compare DOT images of the same breast cancer patient to those of other imaging modalities^{9, 22} obtained at different times, because the nonconcurrent coregistration problem presents many challenges. It is therefore desirable to develop quantitative and systematic methods for data fusion that utilize the high-quality data and versatility of the stand-alone imaging systems.

In this work, we introduce and demonstrate a novel software platform for combining nonconcurrent MRI and DOT: the optical and multimodal imaging platform for research assessment and diagnosis (OMIRAD). To our knowledge, OMIRAD is the first software prototype platform capable of fusing and jointly analyzing multimodal optical imaging data with x-ray tomosynthesis and MR images of the breast.^{23, 24, 25, 26, 27, 28} The OMIRAD platform enables multimodal 3-D image visualization and manipulation of datasets based on a variety of 3-D rendering techniques. Through its ability to simultaneously control multiple fields of view, OMIRAD can streamline quantitative analyses of structural and functional data. OMIRAD is the result of four years of work to develop and test a prototype platform specifically designed for multimodal optical data visualization, fusion, and analysis, including the ability to share data and analysis results across several institutions. It is possible that other postprocessing software exists, but to our knowledge, none of it has been integrated into a full-fledged software prototype, nor tested within several leading institutions in optical imaging research such as within the Network for Translational Research in Optical Imaging (NTROI). Our preliminary study takes an important step toward improved diagnosis and treatment of breast cancer patients with DOT and MRI. Coregistration, facilitated by the software platform, combines structural and functional data from multiple modalities. Segmentation and fusion will also enable a-priori structural information derived from MRI to be incorporated into the DOT reconstruction algorithms. The combined MRI/DOT dataset provides information in a more useful format than the sum of the individual datasets, and we expect the platform to have substantial impact on the standardization of diffuse optical imaging systems, and therefore on the translation of optical imaging research prototypes into viable clinical systems.

The work is organized as follows. We first summarize the structure of the DOT system. We then present the multimodal software platform developed to enable nonconcurrent multimodal data analysis. We describe procedures used for image coregistration and segmentation, and we demonstrate the approach using a simulated phantom breast model. Last, we apply the platform to coregister DOT and MRI images of three breast cancer patients.

2.

Diffuse Optical Tomography System Overview

A schematic of the DOT instrument is shown in Fig. 1 .²⁹ This hybrid continuous-wave (cw) and frequency-domain (FD) parallel-plane DOT system has been extensively characterized for breast cancer imaging using tissue phantoms and normal breast.³⁰ The breast is softly compressed between the source plate and a viewing window, to a thickness of $5.5\text{to}7.5\mathrm{cm}$ . The breast box is filled with a matching fluid (i.e., Intralipid and indian ink) that has optical properties similar to human tissue. Four laser diodes (690-, 750-, 786-, and $830\text{-}\mathrm{nm}$ wavelength), amplitude modulated at $70\mathrm{MHz}$ , were used as light sources. We use a grid of $9\times 5=45$ source positions with a spacing of $1.6\mathrm{cm}$ . For (cw) transmission detection, we sample a $24\times 41=984$ grid of pixels from the charge-coupled device (CCD), which corresponds to a detector separation of about $3\mathrm{mm}$ on the detection window. For remission detection (FD), we use a $3\times 3=9$ grid ( $1.6\text{-}\mathrm{cm}$ spacing) of detector fibers located on the source plate. Remission detection is used to determine the average optical properties of the breast. These values are used as an initial guess for the nonlinear image reconstruction. The CCD data are used for the image reconstruction. For each source position and wavelength, FD measurements were obtained via nine detector fibers on the source plate, and cw measurements were obtained simultaneously via CCD camera in transmission. The amplitude and phase information obtained from the FD measurements are used to quantify bulk optical properties, and the cw transmission data are used to reconstruct a 3-D tomography of optical properties within the breast.

Fig. 1

Schematic of the parallel plate DOT instrument with patient subject lying in the prone position (refer to Ref. 29 for the technical details).

To reconstruct the absorption and scattering image, an inverse problem associated with the photon diffusion equation is solved by iterative gradient-based optimization.³¹ The algorithm reconstructs chromophore concentrations $(C_{\mathrm{Hb}},C_{\mathrm{Hb}{\mathrm{O}}_2})$ and scattering coefficients directly using data from all wavelengths simultaneously. We use a variation of the open-source software package Time-resolved optical absorption and Scattering tomography (TOAST) for these reconstructions. TOAST determines the optical properties inside a medium by adjusting these parameters, such that the difference between the modeled and experimental light measurements at the sample surface is minimized.³² Images of physiologically relevant variables, such as total hemoglobin concentration (THC), blood oxygenation saturation $\left(\mathrm{St}{\mathrm{O}}_2\right)$ , and scattering are thus obtained. Figure 2 shows an example of a 3-D distribution of THC $\left(\mu \mathrm{M}\right)$ in a patient breast with an invasive ductal carcinoma.

Fig. 2

3-D distribution of THC ( $\mu \mathrm{M}$ concentration) in a patient breast containing an invasive ductal carcinoma. Consecutive 2-D patient slices are adjacent.

The resulting DOT dataset is a finite element (FE) model containing on average 50,000 nodes and 200,000 tetrahedral elements. Each node is associated with the reconstructed physiological values such as THC and $\mathrm{St}{\mathrm{O}}_2$ . To facilitate registration of DOT and MR images, the FE model is automatically resampled into a 3-D voxelized volume. The smallest bounding box surrounding the FE model is identified; this volume is divided into voxels ( ${128}^3$ by default). Every voxel is associated to the tetrahedral element to which it belongs, and finally, using the element’s shape functions, the correct physiological value is interpolated at the location of the voxel.

3.

Software Platform Description

The OMIRAD platform enables multimodal integration and visualization of data from DOT and MRI. Figure 3 describes a typical workflow that a user can follow:

The visualization platform showing same patient MRI and DOT (blood volume) datasets before registration is exhibited in Fig. 4 . After the appropriate color transfer functions are applied, one can clearly observe the location of the invasive ductal carcinoma diagnosed in this patient breast.

Fig. 4

Visualization platform showing MRI and DOT (blood volume) datasets from the same patient before registration. After the appropriate transfer functions are applied, one can clearly observe the location of the invasive ductal carcinoma diagnosed in this patient breast. The following components are shown from left to right: 1. orientation cube, 2. transfer function editors, 3. data attribute windows, 4. volume rendering window, 5. MPR windows, and 6. command tabs.

4.

Three-Dimensional/Three Dimensional Diffuse Optical Tomography to Magnetic Resonance Imaging Image Registration Algorithm

3-D DOT/3-D MRI image registration presents several new challenges. Because registration of DOT to MR acquired nonconcurrently has not been extensively studied, to our knowledge no standard approach has been established for this problem. In this work, we present a new approach, specifically devised for addressing the challenges of DOT/MR image registration. DOT images have much lower anatomical resolution and contrast than MRI, and the optical reconstruction process typically uses a geometric model of the breast. In our case, the constraining geometric model of the breast is a semiellipsoid. At the University of Pennsylvania, the patient breast is compressed axially in the DOT imaging device and sagitally in the MRI machine, and, of course, the breast is a highly deformable organ.

For this task, we require that registration be automatic with little prior user interaction and be robust enough to handle the majority of patient cases. In addition, the process should be computationally efficient for applicability in practice, and yield results useful for combined MRI/DOT analysis.

At this time, automatic image registration is an essential component in medical imaging systems. The basic goal of intensity-based image registration techniques is to align anatomical structures in different modalities. This is done through an optimization process, which assesses image similarity and iteratively changes the transformation of one image with respect to the other, until an optimal alignment is found.³⁴ Computation speed is a critical issue and dictates applicability of the technology in practice. Although feature-based methods are computationally more efficient, they are notoriously dependent on the quality of the extracted features from the images.³⁵

In intensity-based registration, volumes are directly aligned by iteratively computing a volumetric similarity measure based on the voxel intensities. Since the amount of computations per iteration is high, the overall registration process is very slow. In the cases where mutual information (MI) is used, sparse sampling of volume intensity could reduce the computational complexity while compromising the accuracy.^{36, 37} In Ref. 38, a projection-based method for 2-D/2-D image registration is proposed. In this method, the projections along the two axes of the image are computed. Horizontal and vertical components of the shift are then computed using a 1-D cross-correlation-based estimator. They show the method is robust in the presence of temporal and spatial noise, and computationally efficient compared to the 2-D correlation-based shift estimator. In Ref. 39, the authors formulate 3-D/3-D registration cost function as the summation of three 2-D/3-D optimization cost functions. The optimization is then done concurrently on the sum of the cost functions, which are identically parameterized. Furthermore, images are preprocessed to extract a binary segmentation. Projection images from the binary segmentation are used for computing similarity measures. The key is to choose a well-behaved similarity measure that can robustly characterize a metric for the volumes.⁴⁰ To make such an algorithm practical, the computational time must also be reduced. In Ref. 40, researchers suggest random sampling of the volume datasets and computation performance based only on these random samples to decrease the computational load. In Ref. 41, authors propose a hybrid technique, which selects a set of high-interest points (i.e., landmarks) within the volume and tries to do registration based on those points only.

Let us consider two datasets to be registered to each other. One dataset is considered the reference and is commonly referred to as the “fixed” dataset. The other dataset is the one onto which the registration transformation is applied. This dataset is commonly referred to as the “moving” dataset. Registration of volumetric datasets (i.e., fixed and moving) involves three steps: first, computation of the similarity measure quantifying a metric for comparing volumes; second, an optimization scheme that searches through the parameter space (e.g., 6-D rigid body motion) to maximize the similarity measure; and third, a volume warping method that applies the latest computed set of parameters to the original moving volume to bring it a step closer to the fixed volume.

Our proposed approach is a novel combination of the methods described in Refs. 42, 43: we compute 2-D projection images from the two volumes for various projection geometries, and set up a similarity measure with an optimization scheme that searches through the parameter space. These images are registered within a 2-D space, which is a subset of the 3-D space of the original registration transformations. Finally, we perform these registrations successively and iteratively to estimate all the registration parameters of the original problem.

We further optimize the performance of projection and 2-D/2-D registration similarity computation through the use of graphics processing units (GPU). Details and general validation of this novel approach have been recently presented.⁴⁴ Multiple 2-D signatures (or projections) can represent the volume robustly depending on the way the signatures are generated. An easy way to understand the idea is to derive the motion of an object by looking at three perpendicular shadows of the object (see Fig. 5 ).

Fig. 5

Generation of 2-D signatures from 3-D volumes, in the three mutually orthogonal directions.

Figure 6 provides an illustration of different transformation models used in medical image registration: rigid, affine, and free-form transformations. Nonrigid registration, depending on complexity, can be classified in two ways: 1. affine transformations [see Fig. 6c], which include nonhomogeneous scaling and/or shearing; and 2. free-form transformations [see Fig. 6d], which include arbitrary deformations at the voxel level. These transformations can be based on intensity, shape, or material properties. The dominant transformation observed across the MR and DOT datasets is due to the difference in compression axis (lateral compression for MR versus axial compression for DOT); this transformation can be modeled using affine parameters. DOT images do not possess enough local structure information for computation of a free-form deformation mapping to register a DOT to an MR dataset.

Fig. 6

Illustration of different transformation models applied to a 2-D image: (a) original image, (b) rigid transformation, (c) affine transformation, and (d) free-form transformation.

Given the previous challenges, we used the following parameters in the nonrigid registration algorithm.

Mathematically, the estimate of the nine degrees-of-freedom (DOF) homogeneous transformation matrix ${\mathbf{T}}^9$ is initially given by

Since here the similarity metric is mutual information, i.e., $S^2\equiv h\left(A\right)+h\left(B\right)-h(A,B)$ , Eq. 2 can be rewritten as:

Entropy is a measure of variability and is defined as: $h\left(x\right)\equiv -\int p\left(x\right)\ln p\left(x\right)dx$ , and $h(x,y)\equiv -\int p(x,y)\ln p(x,y)dxdy$ ,⁴⁰ where $p\left(x\right)$ is the probability density function (PDF) of variable $x$ , and $p(x,y)$ is the joint PDF of variables $x$ and $y$ . The entropy $h$ is discretely computed as:

The algorithm flowchart is shown in Fig. 7 . Figure 7a shows the global registration flowchart. For a number of iterations $n$ (typically $n=3$ ), the three mutually orthogonal 2-D signatures are generated (sagittal, coronal, and axial) for both the fixed and moving volumes. After each 2-D signature generation, the moving 2-D signature is registered to the fixed 2-D signature. This process is shown schematically in Fig. 7b, and explained in detail next.

Fig. 7

Flowcharts for the registration process: (a) global registration flowchart, and (b) 2-D signatures registration flowchart.

First, the $\Delta$ variables are initialized:

Then, step $k=1$ to $m$ is as follows.

The variables are initialized at the beginning of the registration process and are typically set to the following: $n=3$ , $m=40$ , $\Delta \text{scale}̱\text{initial}=4\mathrm{mm}$ , $\Delta \mathrm{trans}̱\text{initial}=4\mathrm{mm}$ , $\Delta \mathrm{rot}̱\text{initial}=4\mathrm{deg}$ , $\text{divider}̱\text{threshold}=40$ .

5.

Breast Magnetic Resonance Imaging Image Segmentation

Our proposed segmentation approach is based on the random walker algorithm. In this case, the segmentation technique requires little user interaction and is computationally efficient for practical applications.

This algorithm, originally developed in Ref. 45 and extended Ref. 46 to incorporate intensity priors, can perform multilabel, semiautomated image segmentation. Given a small number of pixels with user-defined labels, one can analytically (and quickly) determine the probability that a random walker starting at each unlabeled pixel will first reach one of the prelabeled pixels. By assigning each pixel to the label for which the greatest probability is calculated, high-quality image segmentation can be obtained (Fig. 8 ). This algorithm is formulated in discrete space (i.e., on a graph) using combinatorial analogs of standard operators and principles from continuous potential theory, allowing it to be applied in arbitrary dimensions.

Fig. 8

Breast MRI 3-D image segmentation based on “random walkers:” (a) segmenting fatty from nonfatty tissue, and (b) segmenting tumor from nontumor tissue.

Usually ${\mathrm{T}}_1$ -weighted MR imaging is performed. Images show lipid as bright and parenchyma as dark. Tumor also tends to be dark. Minimal user initialization is required. We have developed a workflow customized for this type of breast MRI segmentation.

The algorithm takes two minutes on average for a MRI volume of size $256\times 256\times 50$ . This algorithm can be used to distinguish fatty from nonfatty tissue and tumor from nontumor tissue, as shown in Fig. 8. Since the objective of this initial study is to show how nonconcurrent data may be registered and jointly analyzed, we should note that we did not use the results of MRI segmentation to improve the DOT reconstructions. Rather, we used the MRI segmentation to isolate the tumor tissue in the image.

One significant advantage of spatially registering DOT to MRI data is the ability to treat anatomical information from MRI data as prior information in the DOT chromophore concentration and scattering variables reconstruction process. By segmenting fatty from nonfatty tissue in a MR dataset for example, we can provide a-priori data about the tissue, which interacts with light in a DOT imaging device. This information can further be incorporated in solving the inverse problem associated with the photon diffusion equation, and lead to a more precise reconstruction of physiological quantities (such as hemoglobin concentration). Note that in this work, the MR segmentation is not incorporated in the DOT reconstruction process. This approach is the subject of an ongoing study to be published in the near future.

6.

Initial Validation of Registration Using a Simulated Phantom Model

6.1.

Building a Simulated Phantom Model

To obtain reference results, we tested our methodology using a virtual model of the breast. This model consists of a hemispherical form representing the breast and containing a second sphere of twice the background intensity, representing the tumor. The diameter of the tumor is $25.6\mathrm{mm}$ (20% of the spherical form diameter), and the diameter of the spherical form is $128\mathrm{mm}$ (see Fig. 9 ).

Fig. 9

Compressed models: (a) 3-D sagittal perspective view of superimposed MRI (dark volume) and DOT (light volume) models. (b) Sagittal cross section of MRI model going through the center of the tumor. (c) Spatially corresponding sagittal cross section of DOT model.

The semispherical model is first compressed in the axial direction to simulate the DOT image. The initial model is again compressed in the sagittal direction to simulate the MR image. The amount of compression used is 25% for both optical and MR images, respectively, in the axial direction (along the $z$ axis) and the sagittal direction (along the $x$ axis).

For the axial compression, the $z$ component of the voxel size was decreased by 25% and the $x$ and $y$ components are proportionally increased to keep the same volume size as the uncompressed model. The sagittal compression is simulated in a similar way by decreasing the $x$ component of the voxel size by 25%, and the $z$ and $y$ components are proportionally increased to keep the same volume size. The new tumor center position after compression is determined by multiplying the tumor center position in pixels, by the new voxel size (see Fig. 9).

6.2.

Simulation Description

The experiments conducted test the registration algorithm’s sensitivity to changes in breast compression, translation, and rotation between the MRI and DOT datasets, mainly due to patient positioning differences.

6.2.1.

First set of simulations: incremental translations along the $x$ , $y$ , and $z$ axis

Initial translations along the $x$ axis are applied incrementally to the DOT image, i.e., the moving volume. The registration algorithm was tested after each translation. These translations simulate the difference in patient placement between the two image acquisition processes (translations of $+∕-50\mathrm{mm}$ ). The simulation is repeated with translations applied along the $y$ and $z$ axes, and the registration is tested for each translation. Figure 10 shows the visual results of translations along the $z$ and $x$ axes. The MR model (top row) is the fixed volume in the simulation, and therefore remains unchanged. The DOT model is the moving volume shown in the center row in Fig. 10. This center row shows different initial starting points for the DOT model. In Figs. 10c and 10d, the tumor appears very small because the cross sections shown are spatially corresponding to those of the MR model, and show the edge of the tumor. Note that the bottom row (the DOT model after registration) should look as much as possible like the top row (the MR model).

Fig. 10

Simulation of compression and translation. Spatially corresponding cross sections of MRI model (top row), DOT model before registration (center row), and DOT model after registration (bottom row). (a) and (b) coronal cross sections: DOT model is translated $\pm 50\mathrm{mm}$ along $z$ direction. (c) and (d) Axial cross sections: DOT model translated $\pm 50\mathrm{mm}$ along $x$ direction.

6.2.2.

Second set of simulations: incremental rotations about the $x$ axis

Several incremental rotations about the $x$ axis (clockwise and counterclockwise directions) are applied to the DOT volume, and the registration algorithm is tested after each rotation. This is repeated for rotations around the $y$ and $z$ axes, and the algorithm is tested for each rotation step. Figure 11 shows examples of rotations applied and the resulting alignments (rotations of $+∕-18\mathrm{deg}$ are shown). Here again, the MR model (top row) is the fixed volume in the simulation, and therefore remains unchanged. The DOT model is the moving volume shown in the center row in Fig. 11. This center row shows different initial rotations for the DOT model. In Fig. 11c, the tumor appears very small because the cross sections shown are spatially corresponding to those of the MR model, and show the edge of the tumor. Note that the bottom row (the DOT model after registration) should look as much as possible like the top row (the MR model).

Fig. 11

Simulation of compression and rotation. Spatially corresponding cross sections of MRI model (top row), DOT model before registration (center row), and DOT model after registration (bottom row). (a) Sagittal cross sections: DOT model is $\text{rotated}\pm 18\mathrm{deg}$ about the $x$ axis. (b) Coronal cross sections: DOT model is $\text{rotated}\pm 18\mathrm{deg}$ about the $y$ axis. (c) Axial cross sections: DOT model is $\text{rotated}\pm 18\mathrm{deg}$ about the $z$ axis.

6.2.3.

Third set of simulations: incremental axial compression of the simulated diffuse optical tomography dataset

Different incremental amounts of compression are applied to the DOT images in the axial direction (along the $z$ axis). To simulate the axial compression, the $z$ component of the voxel size was decreased by 10% for each test, and the $x$ and $y$ components are proportionally increased to keep the same volume size as the uncompressed model. The range of compression used is from 0% compression to 40% compression with a step of 10% for each simulation. Note, no figure is shown in this section.

7.

Application to Nonconcurrent Magnetic Resonance Imaging and Diffuse Optical Tomography Data of Human Subjects

A study involving three patients was performed. This study provides an initial answer to a vital question regarding MRI/DOT data analysis: how can functional information on a tumor obtained from DOT data be combined with the anatomical information about the tumor derived from MRI data?

Three MRI and three DOT (displaying THC) datasets are used in this experiment.

Fig. 3

Typical user workflow for the OMIRAD platform.

All DOT image acquisitions are similar and show the patient total hemoglobin concentration (THC). The procedure described in the typical workflow (Fig. 3) was used for visualizing, editing, and registering the MRI and DOT datasets. However, since the objective of this initial study is to show how nonconcurrent data may be registered and jointly analyzed, we did not use the results of MRI segmentation to improve the DOT reconstructions.

A quantitative analysis of the resulting data is not trivial. We propose a simple analysis method that provides valuable functional information about the carcinoma. Using the MRI/DOT registered data, we calculate the differences in total hemoglobin concentration (THC) between the volumes inside and outside the segmented tumor, as follows.

The computed quantities are described in Fig. 12 . The results are shown in next section.

Fig. 12

THC distribution in a DOT dataset showing the resulting computed quantities after DOT-MRI image registration.

8.

Results and Discussion

8.1.

Phantom Model Validation

For most registration tasks, the most significant error measure is the target registration error (TRE), i.e., the distance after registration between corresponding points not used in calculating the registration transform. The term “target” is used to suggest that the points are typically points within, or on, the boundary of lesions. The registration algorithm gives us the absolute transformation $T_{\text{result}}$ that should be applied to the DOT volume to be aligned to the MRI volume. This transformation is applied to the tumor center and 26 neighboring points. The points are arranged on a cube in which the tumor is inscribed. The cube shares the same center as the tumor (as shown in Fig. 13 ). It has a side length of $25.6\mathrm{mm}$ (equal to the diameter of the tumor). The point positions resulting from the application of the absolute transformation are then compared to the corresponding point positions resulting from the application of the ground truth transformation $T_{\mathrm{GT}}$ , which gives us the expected point positions. This allows us to determine the average TRE for each simulation. The TRE is computed as the average Euclidian distance between the 27 pairs of points $(P_{\mathrm{GT}}^i,P_{\text{result}}^i)$ :

Eq. 5

$$\mathrm{TRE}=\frac{1}{27}\sum _{i=1}^{27}d(P_{\mathrm{GT}}^i,P_{\text{result}}^i).$$

Fig. 13

26 points arranged on the cube and used to compute the TRE. The tumor is inscribed in the cube, and shares the same center as that of the cube, noted with an ×.

The volume of the tumor after registration is also compared to the initial one, and the percentage error is computed. It is important to note that the range of translations chosen during simulations is $40\mathrm{mm}$ (from $-20\text{to}20\mathrm{mm}$ ) to keep a realistic aspect of the simulations. Indeed, the translations represent the patient displacements during the image acquisition, so a range of $40\mathrm{mm}$ is reasonable. Also, the range of rotations chosen is $36\mathrm{deg}$ (from $-18\text{to}18\mathrm{deg}$ ) for the same reasons as before.

Tables 1, 2, 3 show the percent volume errors and the resulting target registration errors. As can be observed, the algorithm is more sensitive to rotations than translations, as the error exceeds 5% in some instances. This is explained by the fact that the algorithm uses 2-D signatures of the 3-D volume. By applying a rotation to the volume, the shape of the 2-D signature changes, whereas by applying a translation, the signature is moved compared to the volume of reference while keeping the same form. The change in form due to rotation makes the convergence of the algorithm more difficult. However, the higher level of rotations (more than $\pm 10\mathrm{deg}$ ) will seldom be encountered in reality, where patients usually lie prone in a reproducible manner, and will not cause such high levels of initial rotation. We conducted tests at these higher rotations to explore limitations of the registration technique. For certain points, the error rate increases considerably. This is also explained by the use of the 2-D signatures. Indeed, when the displacement of the image exceeds the limit of the projector that captures the signature, part of the information on volume is lost, leading to a potential divergence of the algorithm. Even though the algorithm is not strictly volume preserving, because of the scaling transformation, the volume percent error shows that within the practical range of deformations, the tumor volume is preserved within an average of about 3% of its original size, which is a reasonable error. Finally, the error due to compression is always under 5%.

Table 1

Percent volume errors with respect to the original moving volume, and resulting target registration error due to incremental translations applied along the x , y , and z axes.

Translation amount(mm)	Translation along x axis		Translation along y axis		Translation along z axis
	Volume%error	AverageTRE(mm)	Volume% error	AverageTRE(mm)	Volume% error	AverageTRE(mm)
$-20$	2.97	3.77	2.59	0.60	1.05	0.89
$-10$	2.27	1.76	1.47	0.87	$-1.51$	2.34
0	1.79	2.62	1.79	2.62	1.79	2.62
10	4.51	3.02	2.80	0.71	1.24	1.00
20	3.95	3.03	$-1.43$	3.03	5.55	4.21

Table 2

Percent volume errors with respect to the original moving volume, and resulting target registration errors due to incremental rotations applied about the x , y , and z axes.

Rotation amount(Degrees)	Rotation about x axis		Rotation about y axis		Rotation about z axis
	Volume% error	AverageTRE(mm)	Volume%error	AverageTRE(mm)	Volume% error	AverageTRE(mm)
$-18$	7.52	7.45	2.42	2.66	3.29	10.59
$-9$	10.08	11.31	0.71	0.90	7.00	4.29
0	1.79	2.62	1.79	2.62	1.79	2.62
9	2.88	2.58	4.77	2.98	5.67	1.77
18	0.70	0.70	$-0.34$	$-0.34$	4.24	4.24

8.2.

Patient Study

Figure 14 shows the computed statistical quantities as well as the difference measures. As expected, all DOT datasets show average tumor THC values, one to three standard deviations higher than the average breast THC values. The results also show large variability in average breast THC values from one patient to another (varying from $21\text{to}31\mu \mathrm{M}$ ). This justifies the use of our difference measure $\mu$ , which defines a normalized quantity allowing interpatient comparisons.

Fig. 14

Statistical values computed in the registered DOT datasets: the segment middle points are the average THC values ( $\alpha$ inside the breast, $\beta$ inside the tumor). The segment endpoints represent one standard deviation spread $\sigma$ and $\mu$ is the difference measure (distance from $\alpha$ to $\beta$ in terms of $\sigma )$ .

These results confirm that the tumor areas in the patient breasts exhibit significantly higher THC than their surroundings. Figures 15, 16, 17 show the visual results of the registration algorithm when applied to real patient datasets. We show superimposed MRI and DOT images (3-D renderings and 2-D fused images) before and after registration. As can be qualitatively ascertained from the figures, registration has greatly improved the alignment of the DOT and MRI datasets. The images also show a significant overlap between the location of the tumors in the MRI and DOT datasets. Patient 3 shows particularly good correlation between the two modalities.

Fig. 15

Superimposed MRI and DOT images (3-D renderings and 2-D fused images) of patient 1 before and after registration. The 2-D fused images show the cross sections going through the center of the tumor.

Fig. 16

Superimposed MRI and DOT images (3-D renderings and 2-D fused images) of patient 2 before and after registration. The 2-D fused images show the cross sections going through the center of the tumor.

Fig. 17

Superimposed MRI and DOT images (3-D renderings and 2-D fused images) of patient 3 before and after registration. The 2-D fused images show the cross sections going through the center of the tumor.

The combination of DOT and MR image resolution, the registration technique, and the segmentation accuracy in MR all affect the final outcome. Can registration errors significantly affect the quantification of the results? Certainly this is possible, but a larger-scale study is required to better characterize the effect of the registration error. Indeed, this issue is significant when automatic segmentation is employed and when quantitative values are derived, especially in the case of small pathologies. Variations in the target registration error (TRE) cause variations in the overlap of the MR segmentation to the THC in the DOT dataset, which in turn cause variations in the quantification of the computed difference measure $\mu$ . However, because the THC is a slowly varying quantity in the DOT dataset, we expect small variations in $\mu$ .

To test this hypothesis, we simulated variations in the TRE by translating incrementally the MR segmentation area in the direction of maximum THC gradient in the DOT dataset. This enabled us to assess the upper bound of the quantification error due to TRE variations. The MR segmentation area was translated 1, 2, 3, 4, and $5\mathrm{mm}$ . Then the different statistics were computed again, and variations in the difference measure $\mu$ are shown in Fig. 18 .

Fig. 18

Variations in $\mu$ for each patient, due to translations of the MR segmentation area inside the THC DOT dataset.

Table 3

Percent volume errors with respect to the original moving volume, and resulting target registration errors due to incremental axial compression of the moving volume.

% Amount ofaxialcompression	Volume % error	Average TRE (mm)
0	4.91	1.21
10	3.66	0.87
20	1.80	1.06
30	$-0.36$	2.37
40	0.11	2.76

As Fig. 18 shows, in all cases the difference measure decreases in amplitude as the translation distance is increased. This shows that the MR segmentation area is translated away from the THC “hotspot” in the DOT datasets. The variations of $\mu$ from the baseline $(\text{translation}=0\mathrm{mm})$ in all cases are less than 15%, and $\mu$ remains equal to or larger than 1, i.e., the average THC inside the segmentation area remains more than one standard deviation away from the overall dataset average THC. Even though these results are limited to only three patients, they exhibit a relative robustness of the registration-segmentation-quantification approach to errors in automatic registration and segmentation. It is also worth noting that these results may apply more generally to patients with breast cancer tumors of sizes within the size range tested, between 2 and $5\mathrm{cm}$ , which is typical.

The most efficient way to improve the coregistration technique in the near future will be to provide additional structural information on the DOT dataset. One way to achieve this goal is to provide a more accurate surface map of the patient’s breast as it is scanned in the DOT device, using stereo cameras for example. A more precise surface map may enable us in the near future to use more complex voxel-based and/or physics-based nonrigid registration techniques, and achieve better MR-DOT dataset alignment, especially for patients whose cancer tumors are smaller than $1\mathrm{cm}$ in size.

9.

Conclusion

A software platform is developed and tested for multimodal integration and visualization of DOT and MRI datasets. The platform, OMIRAD, enables qualitative and quantitative analysis of DOT and MRI breast datasets, and combines the functional and structural data extracted from both types of images. To our knowledge, no other software systems combine nonconcurrent DOT and MRI data, as this developed platform does. The present work introduces this platform and its applications. Multiple simulation results are obtained and demonstrate the robustness of the algorithm.

The method integrates advanced multimodal registration and segmentation algorithms along with a straightforward and well-defined workflow. The method requires little prior user interaction, and is robust enough to handle a majority of patient cases, computationally efficient for practical applications, and yields results useful for combined MRI/DOT analysis. This method presents additional advantages:

We also present an initial patient study that explores an important question: how can functional information on a tumor obtained from DOT data be combined with the anatomy of that tumor derived from MRI data? The study confirms that the tumor areas in the patient breasts exhibit significantly higher THC than their surroundings. The results show significance in intrapatient THC variations, and justify the use of our normalized difference measure.

The completion of this study and the use of the OMIRAD platform contribute to the important long-term goal of enabling a standardized direct comparison of the two modalities (MRI and DOT). We anticipate it will have a positive impact on standardization of optical imaging technology through establishing common data formats and processes for sharing data and software, which in turn will allow direct comparison of different modalities, validation of new versus established methods in clinical studies, development of commonly accepted standards in postprocessing methods, creation of a standardized MR-DOT technology platform, and eventually, translation of research prototypes into clinical imaging systems.