1.1.

Light Transport in Tissue

Light propagation in tissue is dominated by multiple scattering and in the near-infrared (NIR) part of the spectrum, with an especially low absorption window between the wavelengths of 610 to $940\mathrm{nm}$ . Here, the absorption is orders of magnitude lower than that observed in the visible (blue-green), ultraviolet, or infrared regions. Light transport in all wavelength regions is thought to be accurately modeled by radiation transport methods (both deterministic and stochastic models), which can predict the optical fluence patterns based on knowledge of the microscopic interaction coefficients. These coefficients are the absorption coefficient ${\mu }_a$ , representing the probability of photon absorption per unit pathlength, the scattering coefficient ${\mu }_s$ , representing the probability per unit pathlength of a scattering event, and the microscopic phase function $P\left(\theta \right)$ , which describes the probability distribution for the angular direction of outgoing scattered photons. Radiation transport models can be simplified if the absorption is low and the fluence is being detected in the far field [i.e., distance $\gg 1∕({\mu }_a+{\mu }_s)$ ]. In this case, the fluence can be modeled well by an assumption of isotropic scattering using a modified scattering coefficient, called the transport or reduced scattering coefficient, defined as ${\mu }_s={\mu }_s(1-g)$ , where $g$ is the average cosine of the scattering angle:

In addition to analytic light transport modeling, stochastic methods to predict light fluence have been developed and formalized for many years, ^{6, 7, 8, 9} and have lead to an improved understanding of light transport in tissue, especially over smaller distances where the light is not a true diffusive field. These forward models provide the basis for understanding the physical limitations of imaging with light through tissue,^{10, 11} and are utilized in model-based imaging of tissue with light. Application of this forward modeling approach to clinical or preclinical imaging methods has not been successful yet. Thus, in this review, diffusion-model-based imaging is analyzed and characterized in terms of resolution, contrast, and detectability of objects.

1.2.

Imaging with Light: Planar Techniques

The mechanisms to form functional and physiologically relevant optical images of tissue with NIR light fall largely into two broad categories, namely planar imaging and tomographic image reconstruction. In planar imaging, similar to single photon emission corrugated tomography (SPECT) nuclear medicine imaging, the remitted light from tissue is imaged with a planar device, such as a charge-coupled device (CCD) when using optical photons, which inherently captures the image pixilated by the available detector resolution. This type of imaging is predominant in clinical procedures such as endoscopy, colposcopy, and ophthalmology,^{12, 13, 14} and experimental studies in reflectance, fluorescence, and bioluminescence imaging. The image resolution and performance of these procedures is highly dependent on the geometry and the specific optical system used to acquire data. In fluorescence imaging, the issue of background suppression is most often the dominating factor affecting image quality, since rejecting the excitation light from a diffuse body is a challenging problem. In bioluminescence imaging of tissue, though there is little background signal, the intensity of the emitted light is orders of magnitude lower than that seen in intrinsic optical tomography,^{15, 16, 17} and thus low signal to noise typically dominates imaging performance, but the noise is from the readout or dark noise of the camera itself. In cases with ample signal, the shot noise will dominate, but this is often not the case in practical imaging where the light signal is low. In all of these applications, imaging performance is still further affected by the external shape of the tissue and the system’s ability to compensate for irregular boundaries, as well as the light propagation through overlying tissues. Indeed, accurate knowledge of these parameters can improve the resulting images, and throughout NIR imaging of tissue, the problems associated with knowledge of the tissue boundaries and optical properties have been a key problem. As is discussed in Sec. 2, even though the hardware may have good resolution in the absence of scattering, ultimately as the depth increases the resolution of localized emitting tissues below the surface will be governed more by the light transport than by the specifics of the hardware.

In absorption-based and fluorescence imaging, there is a unique feature of having to introduce a light field into the tissue, and then having to remove this light field from the tissue prior to detection of the true remitted signal. The most significant problem of all is specular reflectance from the surface, which can be orders of magnitude larger than signals remitted from the interior of the tissue. To reduce specular reflection, using a crossed polarizer is standard, yet this does not completely remove the reflection, and it still does not compensate for differences in remittance due to the curvature of the tissue. Interpreting all of these issues in the presence of light that has traveled an unknown indirect path between the time it went into and was remitted out of the tissue is problematic. While radiation transport models can be used, they become most challenging for prediction at complex surfaces, and the entire interpretation of background, reflections, tissue coupling, and radiation transport can significantly degrade image quality if not carefully designed. Here, the focus is on using image analysis tools to interpret image quality, which may or may not correctly account for or compensate for all of these nontrivial issues. Planar imaging in particular as compared to tomographic imaging is prone to errors due to curvature of the tissue boundary. Corrective methods have been successfully proposed, such as using a phase signal to correct for tissue curvature,¹⁸ focusing on the second derivative of the overall remission image,¹⁹ or utilizing scanner systems to outline the surface volume.²⁰ Thus the problem is not insurmountable, yet in many commercial systems this issue is largely ignored. Again, as in biolumeinscence imaging, the deeper that the objects being imaged are located, the more they are affected by blurring from the light transport process, and inherently the lower the resolution will be. A more comprehensive analysis of this is provided in Secs. 2, 3.

1.3.

Imaging with Light: Tomographic Techniques

In tomographic imaging, the remitted signals are acquired, and an analytical or numerical algorithm is used to calculate images of the interaction coefficients or the distribution of tissue constituents that would give rise to the measured data. In-vivo absorption tomography, scatter tomography, ^{21, 22, 23, 24, 25} and fluorescence tomography ^{26, 27, 28, 29} are now all experimentally demonstrated, and tools for quantitative analysis of these systems are required. The issues involved in interpreting tomographic images are distinctly different from planar imaging, since the assumption of a nonscattered field for resolution assessment is rarely possible. Indeed, resolution assessment itself implies that the contrast is high and scattering is insignificant, such that the dominant factor is the physical limitations of resolution. Computed tomographic imaging is almost always done in a manner that sacrifices spatial resolution for improved contrast resolution; an unavoidable tradeoff due to the need to choose a finite number of projections to keep the computational burden within what is feasible. In comparison, projection tomography without computational recovery can often use an arbitrarily large number of projects to maximize detectability, as there is little computational burden. However, in both cases, detectability of the object is more important than resolution assessment. Thus, characterizing contrast resolution becomes more important than high resolution assessment of imaging performance. This assessment is readily achieved through contrast-detail analysis of the imaging system discussed in detail within Sec. 3.

In many cases, the recovery of images is done iteratively using a Newton method, requiring inversion of a highly ill-posed and ill-conditioned matrix. After acquisition of a set of measurements of the outgoing optical flux at the boundary ${\Phi }_m$ , these data are simulated in a calculation numerically or analytically with a light transport model ${\Phi }_c$ . This calculation is based on an initial estimate of the property distribution, and then the difference between the ${\Phi }_m$ and ${\Phi }_c$ values is minimized by iteratively updating the set of spatial absorption, scattering, or emission coefficients $\mu$ . This is done by solving a matrix equation by precalculating a linear sensitivity matrix, also known as the weight matrix or Jacobian matrix $\mathbf{J}$ , which is used to map a change in $\mu$ to a change in $\Phi$ , providing the update equation:

2.1.

Line Pair Imaging Analysis

Line pair analysis is the standard in most imaging systems, where the assumption of extremely high contrast makes sense, and where the ultimate spatial resolution of the system needs to be tested. Phantoms or test fields that do not cause significant background scatter, and have effectively infinite attenuation contrast, are standard tools to assess the limiting spatial resolution. The United States Air Force (USAF) resolution test chart, shown in Fig. 1a , is classically used as a universal test field containing a range of line pairs per millimeter. It is the ideal test field for mesoscopic or telescopic imaging fields, where the field size ranges from the order of millimeters to near $5\mathrm{\mu }\mathrm{m}$ . Below the micron range though, this test field does not provide sufficient resolution for useful imaging analysis, and precalibrated microscopic line pair test fields can be purchased. Resolution limits for a given imaging system are measured by visually assessing the number of line pairs per millimeter that are discernable in an image. More quantitative measures can also be applied, such as determining the line pairs per millimeter grouping that shows a discernable decrease between the two dark lines. In the limit of high resolution coherent light, this could be determined with the Rayleigh criterion, where the peak of one object appears to overlap the first Airy disk trough of the other.³⁸

Fig. 1

Images of the USAF test field, imaged (a) through air and (b) through $2\mathrm{mm}$ of 1% intralipid solution. The ability to use this type of resolution tool in a diffusing medium is very limited.

This method of line pair imaging analysis has been used in several biomedical papers mainly focused on thin tissue imaging. Examples include investigations into the use of polarization filtering,⁴⁰ ultra-fast time-resolved detection,⁴¹ or coherent gating methods,^{42, 43} all of which focus on rejecting multiple scattered light from transmittance or reflectance signals.^{44, 45} The use of line pair resolution testing at depths beyond the transport scattering length $({\ell }^{\prime }=1∕{\mu }_s^{\prime })$ is not very informative, since the anticipated resolution is larger than the maximum width between line pairs in the USAF test field. Line pair analysis is usually restricted to imaging in high resolution situations, since imaging multiple lines requires a large field of view that has high resolution and can thus image many lines within a field. For imaging systems with lower to moderate spatial resolution, or systems that image in a scatter dominated field, it is more useful to discuss the resolution as assessed by the point spread, line spread, or edge spread function, or as interpreted by the Fourier transform of these, the modulation transfer function. This is discussed in the next section.

2.2.

Point/Line/Edge Spread Functions

Assessment of the point spread function (PSF), line spread function (LSF), or edge spread function (ESF) can provide nearly equivalent information in imaging systems that have a linear and spatially homogeneous response. Unfortunately, this latter criteria is rarely true in most useful imaging systems, and so a comprehensive assessment of these systems would involve assessing the PSF at multiple locations in the imaging field⁴⁶ to assess the spatial response. For example, chromatic aberration effects lead to imaging field responses that vary in space, and so assessment of the PSF at different areas in the field is often required to characterize a system. Similarly, in diffuse imaging, the point spread function can easily be distorted by the boundaries of tissue, as the diffuse response function is altered by orders of magnitude when near a boundary.

Point spread function imaging in tomography has been a standard practice to assess imaging system and algorithm quality.⁴⁷ However, many studies have blurred the line between point imaging and circular region imaging when defining spatial resolution. When infinitely high contrast is used, the size of the object being imaged should not affect the size of the resulting image, but rather as the object size decreases, the resulting image response decreases in magnitude. PSF determination can be achieved using black spherical inclusions in the medium to simulate the small infinite absorbers present. These types of studies are tedious, but when done exhaustively can provide fundamentally important information about the imaging response and the field.^{10, 48, 49}

Measurements of the line spread function have been completed in many early studies of diffuse tomography, and illustrate the classic banana-shaped sensitivity function observed between the source and detector. As a line is translated across the path between source and detector and oriented perpendicular to the path of the light travel, the change in intensity obtained in this profile study is exactly proportional to the intensity of the photon path, or the sensitivity function of each source-detector pair. ^{10, 50, 51, 52, 53, 54, 55} These basic measurements can be obtained and compared to computational predictions to illustrate the accuracy of the forward model in predicting the Jacobian matrix.⁵⁴

One useful empirical observation applicable to imaging through diffusely transmitting slabs is that the lowest spatial resolution is most likely limited by the broadest point of the photon path, which has been estimated at 20% of the slab thickness.⁵⁶ While this estimation is empirical and only true in diffuse regimes, it provides a “rule of thumb” to estimate the resolution at the center of a slab for a point source and point detector. There are analytic predictions of the photon paths for both reflectance and transmittance imaging.^{54, 57} These predicted paths are effectively weights representing the ensemble of statistical paths taken by photons, and are analogous mathematically to the adjoint Jacobian matrix, as developed for image reconstruction.^{4, 33} These mathematical expressions and algorithms provide a useful way to estimate the lower resolution limit for given source-detector and boundary geometries.

The use of an edge spread function provides a strong contrast that is readily imaged, and analytic derivations of the transmission past an opaque edge in an otherwise scattering medium have been shown^{58, 59, 60} for infinite medium geometries. Measurements in bounded tomography regions show considerably less predictable response,⁶¹ but still provide strong insight into the nonlinear response across the imaging field. The disadvantage of this approach is that most diffuse reconstruction programs are inherently based on the concept of the Newton method, where the background medium is assumed to be homogeneous and the reconstruction algorithm function is to recover the embedded regions within this. This scenario is not true when a large and distributed field such as the edge of a large absorber is present, and so the response of most diffuse image reconstruction algorithms to a large flat object can be substantially different than the response to smaller round-shaped objects. Therefore, results of an ESF assessment in diffuse imaging should probably be interpreted alongside similar measures of the PSF or LSF.⁶¹

2.3.

Modulation Transfer Function

In most medical imaging settings, the data of PSF/LSF/ESF are interpreted in the frequency domain by Fourier transforming the dataset to provide a modulation transfer function (MTF). While the information content is the same, representing the spatial frequencies provides a direct linear method to spatially filter or modify the response function. In many systems, the standard approach has been to use a measurement of the LSF and transform this to the MTF.⁶² However, in certain systems, generating a line that is thin enough to test the system satisfactorily may be problematic due to constraints on the setup, and so it is often simpler to generate a sharp edge for measurement of the ESF. The MTF for each of the PSF/LSF/ESF curves in Fig. 2 are shown in Fig. 3 . It can be seen that similar information is provided in both the spatial and frequency domains for this simple case. In this case, the response at the center of the phantom is seen to have lower spatial frequencies in content, corresponding to wider values of the PSF/LSF/ESF functions. This difference is understood to be caused by the diffusive path between source and detector, which is widest when most distant from a source or detector. Narrowing the source detector distance is the only physical way to decrease PSF/LSF/ESF values, and points nearest either source or detector will always have the lowest values. Decreasing scatter or increasing absorption will reduce the PSF/LSF/ESF values and increase the spatial frequency bandwidth as well, throughout the entire imaging field.

Fig. 2

Graphs of the (a) point spread function (PSF), (b) line spread function (LSF), and (c) edge spread function (ESF) for “pencil-beam” transmission through a 60-mm slab, having diffuse interaction coefficients of ${\mu }_a=0.01{\mathrm{mm}}^{-1}$ ${\mu }_s^{\prime }=1.0{\mathrm{mm}}^{-1}$ . In all three cases, two locations were analyzed using a target near the edge ( $2\mathrm{mm}$ inside the surface) and then at the center of the slab.

Fig. 3

MTF profiles of the graphs shown in Fig. 2, illustrating the information contained as a function of spatial frequency, with LSF having highest resolution, PSF having next highest, and ESF having the lowest. Resolution is always considerably worse in the interior of the diffusing medium than at the edge near a source or detector.

2.4.

Application of Resolution Testing in the Field of Diffuse Tomography

In diffuse light imaging, it has long been recognized that light follows a statistical path in which the predominant path between source and detector is a line surrounded by a banana-shaped distribution.⁶³ This spreading of the photon paths is induced by the inherent multiple scattering present, and decreases in width in a medium with lower scattering or increased absorption. The effect of increased absorption is somewhat counterintuitive, but generally leads to a loss of photons that have traveled farther in tissue that subsequently narrows the average path of travel. These distributions have been studied by many investigators, and specifically quantified by papers in the early 1990’s. ^{30, 51, 52, 56, 64, 65} Imaging of edges has not proven all that useful, as the wide spread of photons really limits the ability to visualize the edge of objects clearly, and the spatial variation in the resolution ultimately complicates the analysis.

In diffuse tomography imaging, it is easier to resolve a smooth circular heterogeneity embedded in a field than step changes.⁶⁶ This is because of the fact that heterogeneities appear as symmetrically Gaussian filtered objects in the image. Almost all papers in the field of diffuse optical imaging have focused on assessing resolution by placing point objects or line objects in the field to assess spatial resolution. ^{4, 30, 64, 67, 68, 69, 70, 71, 72, 73} This focus has emerged from a fundamental limitation in the field of diffuse tomography stemming from the fact that all currently used reconstruction algorithms are derived at some level from perturbation theory. The Born, Rytov, and Newton methods for minimizing an objective function are all based on perturbing an initial field to find the solution. This approach is inherently optimized for imaging point objects, and the ill-posed nature of the problem, combined with significant regularization, leads to a solution that is significantly smoother than the original test field.

Once a system or algorithm is established in its ability to recover point objects, extension to multiple objects has been a common theme; however, this step is both important and problematic. The most significant problem is the nonlinear response of the measured field to multiple or extended inhomogeneities, requiring an infinite number of heterogeneity configurations to fully analyze system performance.

Perhaps the only reasonable approach to characterizing the imaging field response to multiple heterogeneities is to simulate the expected distributions of values possible in vivo and use this as the limited calibration of the system and corresponding algorithm. Even with these measures taken, it is critical to evaluate these distributions with the full range of object sizes and contrasts expected in vivo, as is discussed in Sec. 3 on contrast-detail analysis.

2.5.

Analysis of Luminescence and Fluorescence Diffuse Imaging Resolution

Imaging of the minimum spatial resolution is only reasonable when an effectively infinite contrast is expected. Fluorescence protein imaging or bioluminescence imaging are two of the few situations in optical imaging in vivo where it may be reasonable to expect nearly infinite contrast, when the background emission issues might be neglected or corrected. ^{74, 75, 76, 77} When cells are specifically transfected or modified to express an optical signal, such as a specific organ or a tumor, the background emission in the neighboring organs should be effectively zero. In green or red fluorescent protein (GFP or RFP) imaging, the background and leakage of excitation light through the filters does provide the most significant background signal; however, this can be significantly reduced when wavelength-dependent fitting or wavelength-based background subtraction is used. In bioluminescence, little real background is present in most cases, and background is often simply the dark noise in the camera or light leakage into the enclosure from the room. Thus, the spatial resolution of bioluminescence or fluorescence protein imaging in vivo can be assessed by point spread function or line spread function imaging, yet little study of this has been reported. One comprehensive paper on this issue by Troy ⁷⁷ showed effective point spread functions as measured in phantoms and in vivo, using small numbers of cells to assess the minimum detectable number of photons and cells. This analysis illustrated that bioluminescence is a more sensitive imaging technique in the remission geometry, by a considerable margin, due to the decrease of fluorescent protein imaging sensitivity caused by background autofluorescence. However, recent reports of fluorescence imaging in the transmission geometry will likely be more sensitive. In most applications of fluorescence or bioluminescence, the actual resolution was not the most important parameter in distinguishing the two systems, but actually the sensitivity. Resolution of bioluminescence and fluorescence appeared to be similar, because the photon spread within a spectral window was effectively equivalent. This study and other similar studies focus on minimum contrast or signal detectable, because the issue of resolution is not governed by system constraints, but rather by the physical constraints of the light transport in tissue. Resolution limits in this regime have less to do with system design than with the depth of the object to be resolved in the tissue. While the resolution of objects at the surface of a tissue can, in principle, be as high as the diffraction limit of light (i.e., near $250\mathrm{nm}$ ) given sufficient contrast, the presence of tissue motion and the quality of the imaging system typically contribute to the real imaging resolution being lowered to typically near 1 to $2\mathrm{\mu }\mathrm{m}$ when imaging at the surface of the tissue. Again, resolution clearly degrades by orders of magnitude in just a few millimeters of depth into the tissue, due to the overwhelming presence of scattering.

3.1.

Contrast-Detail Curves

Contrast-detail (C-D) analysis is commonly used to determine the performance of medical imaging systems and is an effective method for assessing the imaging capabilities of prototype systems. This technique was introduced to medical imaging in the 1970’s and has since been used extensively for CT ^{78, 79, 80, 81} magnetic resonance imaging,⁸² ultrasound,⁸³ mammography,^{84, 85} fluoroscopy,⁸⁶ whole body x-ray systems,⁸⁷ as well as imaging displays.^{88, 89} This technique is used to quantify the combined performance of the imaging system and the image reader in detecting objects representing a clinically relevant range of sizes and contrasts within a domain, focusing on assessing the lower limits of each possible range.^{90, 91} A C-D graph of minimum detectable contrast level for all sizes of objects provides limiting data on two major regimes of system operation, namely 1. the spatial resolution for high contrast objects (high contrast, small object size), and 2. the lower level of contrast detectable for larger-sized objects.⁸¹

Typical contrast-detail study test fields contain a series of objects representing a range of contrasts and diameters often in a regularly spaced pattern. The phantom can have either discretely sized objects or continuously varying sizes, and each object size is repeated with varying contrasts relative to the background. A theoretical representation of a C-D phantom is shown in Fig. 4a , with object size increasing from left to right and contrast increasing vertically. The phantom is imaged tomographically by simulating x-ray computed tomography transmission data with the addition of varying amounts of Gaussian distributed noise in each detector. The noise and geometry of the system contribute to the resulting image quality, as represented by the recovered images in Figs. 4b through 4f.

Fig. 4

An example of contrast-detail analysis is shown for a simulated computed tomography test phantom having six different sized objects at six different contrast values from the background. This type of phantom is used routinely to test imaging system performance, by imaging and determining the minimum contrast detectable for each object size. As seen in (b) and (c), the reconstructed images show a degradation of the image quality due to noise and sampling error, and in (c) the line of minimum detection is shown. In (d), (e), and (f), a 5% noise level was used, and the number of sampling (or projection) angles was systematically decreased to show the degradation of the image and the resulting increase in the position of the contrast-detail curve.

In assessing each image, the lower limit of detection for a given size is determined for each column, so the minimum “detectable” contrast level is determined for each size. This is repeated for all sizes, and a curve (as shown in the images of Fig. 4) is displayed representing the minimum detectable range of contrasts for all sizes. In this analysis, “contrast” represents the real object contrast prior to image reconstruction, relevant to the imaging modality, rather than the contrast within the generated image. In Fig. 4, it can be seen that an increase in noise results in an increased contrast required to detect an object of a given size. Furthermore, a decreased number of projection angle samples results in a higher C-D curve, and therefore lower contrast resolution. In tomography systems, the ability to resolve subtle contrast changes is typically most important, and a lower C-D threshold curve indicates the better system contrast resolution and therefore a higher level of imaging performance.⁹⁰

3.2.

Contrast-Detail Analysis in Medical Imaging

Detection thresholds usually are determined by human readers,⁹² but automated assays of the image can also be based on contrast-to-noise ratio (CNR) and have been used in prototype systems. ^{93, 94, 95, 96} In the clinical setting, human readers are the gold standard through which comparisons are made, since all radiological images are read by radiologists for diagnosis. There has been significant research to develop automated algorithms that systematically “detect” objects by acting as ideal observers or appropriately mimicking human observers.⁹²

Contrast-detail phantoms are commercially available for all conventional clinical imaging systems, and the American College of Radiology has established guidelines for phantoms, which should be used for calibration and validation of specific imaging modalities. Common applications of these studies include periodically assuring clinical system image quality, optimizing developing technologies, and comparing intersystem performance. Figure 5 is an example of contrast-detail curves illustrating the typical performance of a mammography system⁸⁵ and an x-ray CT system.^{97, 98, 99} Objects with sizes and contrasts in the range above and to the right of the curves are considered detectable, while those below and to the left are too small or have too little contrast to be detected in the image. Therefore, imaging systems with curves closer to the $x$ and $y$ axes indicate a lower contrast detection threshold for all objects and are considered to have better imaging performance for the phantom being imaged.

Fig. 5

Contrast-detail curves showing the values for a mammography system and for a CT system, illustrating the inherent strengths and weaknesses of each system. The mammography system is strongest for imaging smaller objects, below $2\mathrm{mm}$ diam, as lower contrasts can be detected. In comparison, CT cannot image objects smaller than $1\mathrm{mm}$ , but can image subtle changes in contrast better when the objects are larger than $2\mathrm{mm}$ .

As a performance measurement tool for standard radiography and mammography, contrast-detail analysis is used for scheduled quality assurance, optimizing system settings,^{98, 99} assessing and comparing digital and film-screen detection methods, ^{87, 100, 101, 102, 103, 104} and determining minimum requirements for image viewing and digital file storage protocols. ^{88, 105, 106, 107} Contrast-detail analysis is an efficient means to complete these studies, since a given C-D plot often requires only a single image and minimal time commitments by professional readers. Some C-D studies in radiography seek to optimize the tradeoff between image contrast and x-ray dose to the patient by extending the analysis to include total dose levels. This is particularly important for computed tomography systems, where patient dose is a major concern. A series of studies in the late 1970’s and early 1980’s by Cohen were the first to directly apply contrast-detail analysis to the assessment of CT scanners. ^{78, 79, 80, 90, 108} Faulkner published a fairly comprehensive set of contrast-detail results to explore the effect on CT imaging performance using different filtering techniques, reconstruction algorithms, and CT scanners.⁸¹ Contrast-detail analysis has also been used in assessing ultrasound, ^{83, 109, 110, 111, 112} fluoroscopy, ^{86, 113, 114, 115} and magnetic resonance (MR) systems.^{82, 116} The key in this process is to utilize a contrast-detail phantom, which is representative of the size and contrast scale of the tissues that the system will be used to image routinely. The shape of the C-D curve then allows comparison between systems with the same settings and analysis of which applications would be suitable for the system.

3.3.

Contrast-Detail Applications in Planar Versus Tomographic Imaging

Contrast-detail studies for planar imaging systems with relatively large and flat response fields can be completed using a phantom containing the full set of object contrasts and diameters. This allows the full C-D analysis to be performed with a single, or limited number of, image(s). Applying contrast detail to NIR tomography presents a more difficult problem, since the presence of multiple optical contrast objects in the test field has a substantial effect on image quality¹¹⁷ due to the inherent soft field nature of the problem. However, this has been addressed by using a series of images, each with a single object within the test field.⁹⁴ The object size and contrast are varied between images, and detection thresholds are extracted and compiled to produce the C-D curve. Additionally, simulated studies are easily completed to quantify expected contrast-detail system performance for best case and a variety of other simulated conditions.^{95, 96}

Further imaging performance assessments for emerging technologies should be reported following current medical imaging protocols. Contrast-detail analysis can provide a reasonably comprehensive method for quantitatively and systematically assessing system capabilities. Though other techniques, such as receiver operating characteristic (ROC) analysis, may be more rigorous, these tests are not always practical due to professional reader time constraints. ^{91, 109, 112, 118} The relative ease of producing contrast-detail plots makes this type of analysis attractive as a preliminary or primary method for ranking and optimizing medical imaging systems.

The key factor in the observation of higher spatial resolution in planar imaging is that tomographic image reconstruction spatial resolution is limited by the number of projects that can be obtained, due to the image reconstruction problem becoming overwhelming computationally at higher numbers of projections. This computational limit on the image reconstruction often limits the resolution that can be achieved; however, in projection imaging there is typically no major computation required, and so larger numbers of projections can be taken, and hence the resulting resolution can be maximized for objects near the surface of the tissue. This is certainly true in x-ray imaging throughout tissue volumes, but the scattering process in diffuse tomography complicates this a little further. It is generally believed that diffuse image reconstruction for interior objects (distant from any surfaces) could be recovered with higher resolution in tomography, than with projection imaging.

3.4.

Applications in Endogenous Versus Exogenous Contrast Imaging

With the recent advent of absorption, scattering, and fluorescence tomography methods, ^{26, 27, 119, 120} and potentially bioluminescence tomography,^{121, 122} it is important to fully understand the capabilities and limitations of each approach for characterizing tissues, in terms of the minimum detectable contrast, or cell number in vivo. As stated before, it should be anticipated that in most cases tomography is not the best way to improve spatial resolution, but rather is the optimal way to recover low contrast information (i.e., optimal contrast resolution) when a limited number of projects are able to be obtained due to computational limits. An illustration of this is shown in Fig. 6 , where the contrast-detail curve of absorption-based imaging is compared to tomographic imaging. Recent breakthroughs in enzymatically activatable fluorphores have created a significant interest in the possibilities for molecular imaging.²⁶ Systematic characterization of fluorescence tomography is important for the implementation of this modality to imaging tumors, where specific markers of cellular or vascular expression may be localized.

Fig. 6

Contrast-detail curves are shown for diffuse optical imaging simulations, comparing transmission imaging with projection data through a 6-cm-diam slab, to NIR tomography of an 8-cm cylinder. These geometries are chosen to mimic that of breast imaging, and the test object is placed at the center of the imaging field to simulate the most difficult lesion to detect. The minimum detectable contrast for each size object is shown as a square, illustrating that for reconstruction tomography, the minimum detectable contrast is lower than that for projection imaging. Objects smaller than $4\mathrm{mm}$ are not easily simulated in these calculations, but are presumably more readily detected with projection imaging than with tomographic imaging, as is observed in Fig. 5.

The contrast-detail response of fluorescence tomography is expected to be similar to that of near-infrared absorption tomography, as the path of light propagation is similarly diffuse in nature. Computational studies of this type have recently been completed⁹⁶ and the results, shown in Fig. 7 , illustrate a best-case C-D curve for imaging through $86\mathrm{mm}$ of tissue.

Fig. 7

Contrast-detail analysis of diffuse imaging demonstrating the difference in sensitivity for objects located at the center and near the edge of the imaging field. These simulations are based on fluorescence tomographic imaging, and are similar quantitatively to absorption-based imaging. The simulation field size was $86\mathrm{mm}$ diam, using scattering and absorption parameters typical of soft tissue imaging in the near-infrared ( ${\mu }_s^{\prime }=1.0{\mathrm{mm}}^{-1}$ , ${\mu }_a=0.01{\mathrm{mm}}^{-1}$ ). The points present the contrast at which a minimum contrast-to-noise ratio of 3 is recovered in the resulting tomographic images.

Future studies with experimental systems can demonstrate how the experimental apparatus and extension to real tissue affects the contrast-detail performance reported in Fig. 7. Ultimately, this will help determine the potential role of the system in a research or clinical setting.

Though most current studies in fluorescence optical tomography focus on small animal studies, the high tissue absorption and small size of these test subjects results in a less diffuse light field, which departs from the diffusion approximation of photon propagation and thus can complicate the image reconstruction process. This problem has been well studied, yet no clear solution exists other than attempting to model the light propagation with radiation transport theory or Monte Carlo models.^{10, 123} An emerging approach is to incorporate high absorption coefficients into a diffusion-based model. This technique has been reported with visible wavelengths and has been shown to work considerably well in small animal imaging.¹²⁴ Further work in this area is ongoing, and a clearer analysis of contrast-detail characterization of these systems would be a significant benefit.

3.5.

Contrast-Detail Application in Assessment of Hybrid Imaging Systems

In recent years, there has been a significant interest in developing near-infrared imaging systems that are coupled to standard clinical systems, or systems that can contribute structural prior information to the image reconstruction process. This has been shown with MRI ^{21, 125, 126, 127, 128} and ultrasound.^{129, 130, 131} The hypothesis driving the development of these hybrid techniques is that the accuracy of the image or the image property values will be somehow improved. While this is implied, it has not clearly been proven, and must be demonstrated for each system. Indeed, there is thought to be a delicate balance between forcing the solution to converge to a solution imposed by the a priori constraints, versus allowing the constraints to minimally coerce the solution toward the most accurate image. In recent studies, Brooksby ^{21, 125, 126, 127, 128} have demonstrated in-vivo imaging with MRI-coupled NIR tomography, and used the MRI information to segment the adipose from the glandular tissue, providing input spatial information about the two tissue regions. This approach provides important a priori information, which then improves the image reconstruction. Analysis of this improvement can be accomplished through contrast-detail analysis, thereby quantifying the improvement in contrast that can be expected for each size object, given a decision criterion. This is demonstrated in Fig. 8 with contrast-detail curves for images generated with and without structural a priori information. These curves were determined using $\mathrm{CNR}=3.0$ as the threshold for assessing object detection; however, this study can also be completed using human observers or a computer-generated “ideal observer.”

Fig. 8

Contrast-detail analysis of diffuse imaging with and without a priori information about the fat and fibroglandular layer of the breast, estimating the detectable levels of contrast for given object sizes. With the inclusion of a priori fat/glandular tissue layers, the minimum contrast required for each size is decreased to achieve the same CNR value.

The shift of the curve to the left in Fig. 8 with the inclusion of a priori information from MRI is a quantitative indication that the CNR of the imaging system is superior to the imaging capability without a priori information. The experimental extension of this work is ongoing, and use of contrast-detail analysis in this application will demonstrate improved contrast resolution with increasing levels of spatial constraints that can be implemented. As hybrid imaging systems become more established, as are PET/CT systems, this type of analysis will become increasingly important to quantitatively evaluate the system capabilities.

4.1.

Sensitivity, Specificity, and Receiver Operating Characteristic Analysis

The receiver operating characteristic (ROC) methodology has been widely used to address the clinical efficacy of medical imaging systems. ^{132, 133, 134, 135} In an ROC study, the readers view a cohort of normal and abnormal radiology images and assign numeric ratings (typically four to six) to each image as an indication of their confidence level that the image shows a clinical abnormality. The resulting rating data are then analyzed, summarized, and plotted on an ROC curve. This graph reveals the relationship between the true-positive fraction (TPF) and the false-positive fraction (FPF) as the reader’s confidence level varies. Summary measures of the curve are typically used as an objective measure of the ability of the reader to detect objects in the images, representing the quality of the medical imaging modality when applied in a human diagnosis task. These summary measures include the area under the entire curve and the partial area under the curve in a particular region of interest.

ROC analysis is widely employed to evaluate observer diagnostic performance in a possible situation in which two alternatives exist. In this case, the classification of the stimulus (patients’ real condition) and response (radiologists’ diagnosis) have only one of two possible choices, normal (nondiseased) or abnormal (diseased). For a diagnostic test of $N$ patients, of which $n_1$ patients are abnormal and $n_0$ patients are normal (here $N=n_0+n_1$ ), there is a $2\times 2$ diagnostic table that completely describes the observer’s diagnostic performance, as shown in Fig. 9a . The table is a listing of all the possible combinations for a pair of binary variables, and the data represent the number of occurrences of combinations of the two variables. In Fig. 9a, the two variables are the patient’s real condition and the observer’s diagnosis. The observer gives $m_1$ positive and $m_0$ negative diagnosis readings (where also $N=m_0+m_1$ ) for patients of which $n_1$ are diseased and $n_0$ are nondiseased. Within this, the true positives (TP) and false positives (FP) are the number of diseased and nondiseased patients who are diagnosed as diseased (i.e., positive), respectively, while the true negatives (TN) and false negatives (FN) are the number of nondiseased and diseased patients who are diagnosed as not diseased (i.e., negative), respectively.

Fig. 9

Statistical analysis of data given two distributions. (a) The $2\times 2$ diagnosis table. The observer gives $m_1$ positive and $m_0$ negative diagnosis to $N$ patient images of which $n_1$ are diseased and $n_0$ are nondiseased. (b) Example of probability density distributions of an observer’s confidence in a diagnostic test, which are analyzed by translating the threshold criteria through all possible values, and then plotting the sensitivity versus specificity for the test. (c) The resulting graph from (b) represents the ROC curve, and imaging tests that maximize the area under this curve are considered more beneficial for accurate clinical classification of lesions. The square point on the curve corresponds to the observer’s confidence threshold line shown in (b).

From the $2\times 2$ diagnostic table, the sensitivity and specificity can be calculated and are commonly used as indications of discriminatory accuracy of the diagnostic study. The sensitivity is the conditional probability of a positive diagnosis $(T+)$ , given that the patient is in fact diseased or abnormal $(D+)$ , i.e.:

In practice, the diagnostic discrimination capacity of an observer in a specific diagnostic test is usually not perfect, because diagnoses are made from various states of symptom or evidence. In other words, the diagnosis depends on the confidence level of diagnostic evidence, i.e., the confidence threshold. Thus, it is more informative and meaningful to design the diagnostic tests on a confidence rating scale, either on a fixed number of discrete response categories or a continuous test variable, and then calculate different sensitivity and specificity pairs, which are used to generate a receiver operating characteristic (ROC) curve.

The probability density distribution function of a radiologist’s confidence in a positive diagnosis for a particular diagnostic task is shown schematically in Fig. 9b. The degree of overlap of the diseased and nondiseased distribution functions completely determines the ability of the test to distinguish diseased patients from nondiseased. As shown in Fig. 9b, for a specific decision or confidence threshold value $x_c$ , the sensitivity and specificity values can be calculated. As $x_c$ increases, the specificity increases at the expense of sensitivity. To graphically present the relationship of sensitivity and specificity, the ROC curve is generated [Fig. 9c], which plots sensitivity versus false positive rate (FPR), defined as (1-specificity). This graph provides the sensitivity and specificity for a given imaging system/reader combination for a range of reader confidence thresholds.

4.2.

Interpretation of the Receiver Operating Characteristic Curve

In practice, human-generated ROC curves for a full range of confidence thresholds require a substantial time commitment from the professional image readers involved. This is commonly addressed by reducing the number of confidence thresholds required to be assessed by the reader and completing the curve by using various techniques to model the probability density distributions shown in Fig. 9b. This provides a more continuous representation of the sparse study data. There are several algorithms that have been published to construct ROC curves based on discrete or continuous test data. These algorithms can be divided into two basic categories, nonparametric or parametric, depending on whether the implementation of the algorithm assumes a parametric model. The empirical, nonparametric approach is used to calculate the ROC curve using empirically determined histogram distributions, in which there is no need for structural assumptions and parameters for modeling or fitting. Though the empirical, nonparametric method is easy to implement and robust in general cases, it does not provide a smooth fitted curve and there are no standard statistical measurements, such as confidence levels, available for evaluation. A nonparametric kernel smoothing technique can significantly improve the empirical nonparametric method,^{134, 136} but must be done with a careful interpretation of how the smoothing kernel affects the data. In this method, a kernel function estimating the densities of the distribution functions in the diseased and nondiseased populations and a bandwidth is applied and optimized to numerically represent the distribution functions. This analysis results in the generation of a smooth and optimal ROC curve. Having stated this, parametric modeling of the data is most often chosen and provides a strong proven approach to generating smooth ROC curves with statistically useful estimates of the confidence intervals. The standard normal distribution is most commonly used, and provides a mean and standard deviation value, which are used in generation of the confidence interval lines.

One of the great practical challenges in ROC analysis in a typical observer performance study is how to deal with the variation of skill levels among different observers. Although the basic concept of ROC analysis has been understood since the early 1980’s,¹³⁷ the analytical techniques or research tools available at the time had limited practical applicability, until the introduction of the so-called multiple-reader multiple-case (MRMC) ROC paradigm¹³⁸ in the early 1990’s. The MRMC ROC paradigm uses a theoretical model and proposes applying procedures such as the jackknife method (jackknife readers or jackknife cases) or the bootstrap method, to the area under the ROC curve obtained for each reader. These approaches have allowed inclusion of multiple readers into ROC analysis.

4.3.

Location Receiver Operating Characteristic Analysis

In standard ROC methods, the major focus is assessing the diagnostic utility of the medical images, where the complexity of the target object location is often eliminated by clearly specifying the region of interest (ROI) in the images. Recent developments in localization-response ROC (LROC) analysis offers more understanding of medical imaging methodology, in terms of measuring the ability to detect and correctly localize the actual target within the reconstructed image. These developments include simultaneous ROC/LROC fitting¹³⁹ and alternative free-response ROC (AFROC) analysis.¹⁴⁰ The LROC measures the probability of successfully detecting and locating objects within images, versus the probability of falsely detecting objects in normal images, as a function of the detection criteria.

Despite the essential simplicity of the fundamental concept of ROC/LROC analysis, the reproducibility and repeatability of human observers depends on many nonimage-related conditions such as experience, physiological or ambient conditions, number of tests, etc.¹⁴¹ ROC/LROC analysis therefore must be carried out with careful attention to these possibly confounding conditions.

4.4.

Receiver Operating Characteristic/Location Receiver Operating Characteristic Analysis in Diffuse Optical Imaging

In diffuse optical imaging, little attention has been paid to assessing the detectability of objects, as assessed by human observers. Yet ROC analysis can be readily applied here, and has been shown useful in assessing diffuse images of arthritic finger joints.¹⁴² It is also likely that assessing how humans perceive nonlinearly reconstructed images may help interpret how image reconstruction algorithms should be tailored. As part of this assessment, the ROC analysis of simulated NIR tomography images were completed and are shown in Song ¹⁴³ The purpose of the study was to evaluate diffuse optical imaging by determining the limitation of human readers’ ability to detect objects within the image by ROC and LROC analysis. Given a fixed system noise level, four key parameters determine the quality of a reconstructed image: 1. the size of the heterogeneities, 2. the absorption and scattering coefficients’ contrast between the heterogeneities and the homogeneous tissue, 3. the number of reconstruction iterations used in the image formation, and 4. the location of the region of interest (ROI) in heterogeneous images. Human observer performance is reported in terms of the related area under the curve (AUC) value with error correction of ROC and LROC curves.

Typical reconstructed images are shown in Fig. 10 . Each image contains an inclusion with a different size and contrast level located to the right of the image field. These images are shown to illustrate the visible quality of the tomographic images, and how decreasing size or contrast diminishes detectability of the object within the reconstructed image. The human observer analysis was completed on a large number of images, similar to the ones in Fig. 10. The resulting ROC data of area under the curve (AUC) is shown in Fig. 11 , using 600 images for each object size, and four observers, as discussed in Song ¹⁴³ The data are separated into studies, one in which the contrast was held constant while the size was varied [Fig. 11a], and the other in which the effect of varying contrast was considered for a constant anomaly size [Fig. 11b]. For each of these studies, the location of the object was randomly moved around within the imaged domain, and control images with no objects were also used. Not surprisingly, as the object size or contrast decreases, the ability of humans to detect the objects decreases. One interesting observation is that the LROC AUC values drop to almost zero for the smallest anomaly sizes and lowest anomaly contrast values, as seen in Figs. 11a and 11b, while the ROC AUC values remain high in these regimes. These numbers indicate that humans are still able to detect the presence of local abnormalities in the images, even when their ability to determine the location of the object is diminished or near zero.

Fig. 10

Diffuse optical tomography images reconstructed from simulated data showing the effect of ROI size and contrast on image quality. The ROI location was held constant in each test field, as was the reconstruction iteration number, while the size and absorption contrast was varied. In the top row of images, the size was fixed at $12\mathrm{mm}$ diam, and the contrast values were varied in the original data with (a) 1.1, (b) 1.4, and (c) 2.0. In the bottom row of images, the contrast was fixed at 2.0 and the size was varied with (d) $4\mathrm{mm}$ , (e) $10\mathrm{mm}$ , and (f) $16\mathrm{mm}$ diam. The background optical properties were ${\mu }_a=0.004{\mathrm{mm}}^{-1}$ and ${\mu }_s^{\prime }=1.0{\mathrm{mm}}^{-1}$ .

Fig. 11

ROC analysis of diffuse optical tomography images with area under the curve (AUC) values of ROC and LROC and their standard errors. In (a), the heterogeneity size study is shown. The objects all have the same contrast $C=2.0$ , and six iterations were used. In (b), the heterogeneity contrast study is shown. The diameter of ROI was equal to $10\mathrm{mm}$ for all, and six iterations were used. In (c), the reconstruction process iteration number was studied. In (d), the results of a heterogeneity location study are shown, where the object location is varied from the center to the edge of the field. In (c) and (d), the object was a fixed diameter of $10\mathrm{mm}$ and the absorption contrast was varied randomly from 1.1 to 2.0.

In Fig. 11c, ROC and LROC AUC analysis was used to examine the effect of the number of iterations in the nonlinear reconstruction algorithm by presenting the observers with image sets at different levels of iterations. The object size was maintained at $10\mathrm{mm}$ , and the contrast and location of the object were varied randomly. The AUC values for both ROC and LROC are nearly constant, but with a subtle but significant rise in LROC at lower iteration numbers. This indicates that observers are better able to localize the region at fewer numbers of iterations. This may be somewhat counterintuitive; however, it is explained by understanding that even though the modeled and measured data more closely match after more iterations, these reconstructed images may also contain a higher degree of spatial noise.

In Fig. 11d, the effect of the location of the object was assessed by similar LROC and ROC analyses. It is well known that the imaging field response is highly spatially variant in diffuse tomography, and this is demonstrated directly from the images using LROC curves. The ability to localize objects increases (i.e., the AUC of the LROC increases) as the object moves toward the edge of the imaging field. This reflects the fact that the diffuse imaging provides a better reconstruction of objects closer to the edge of the field.

Further use of ROC analysis in detection of objects will allow improved assessment of the factors that limit the use of diffuse tomography methods in detection tasks.

4.5.

Receiver Operating Characteristic Analysis of Mammography with Adjuvant Imaging

Near-infrared tomography is now in clinical trials at several centers for breast cancer tumor imaging. The use of NIR as a screening tool has been the subject of research for several years. A clinical study of multispectral NIR tomography in breast cancer imaging was recently summarized by Poplack ¹⁴⁴ In this study, a six wavelength, multispectral imaging system was used to recover images of breast hemoglobin, oxygen saturation, water, and scattering properties in more than 100 patients. ROC analysis was not completed on this study; however, quantification of the contrast for each tumor region relative to the background breast tissue was completed. The tumor location information from mammographic images, such as size, distance from the chest wall, and depth within the breast, were provided prior to NIR imaging to allow quantification of the ROI values. Typical images of total hemoglobin (HbT), oxygen saturation (StO2), water, scattering amplitude, and scattering power values are shown in Fig. 12a .

Fig. 12

(a) In this figure, a representative image is shown for a cancer case, with a focal increase in hemoglobin, decrease in oxygen saturation, increase in water, and increase in scattering at the site of the tumor. These values were quantified for a series of tumors used in the study, and an ROC curve was developed for classification of cancer tumors relative to benign tumors, based on the hemoglobin concentration value. The curve (b) is shown for both small $(\mathrm{diameter}<=6\mathrm{mm})$ and large $(\mathrm{diameter}>6\mathrm{mm})$ breast cancer tumors, to illustrate that detectability of smaller tumors is not feasible with this type of tomography, but that detection of larger tumors may be clinically beneficial.

Using the compiled values of contrast in the cancer and benign tumors, the curve of sensitivity and specificity was generated using a parametric analysis of the relative contrast data. The decision criterion was continuously varied to draw the complete ROC curve. Figure 12b shows the ROC curves based on the normalized hemoglobin values of both cancer tumors relative to benign tumors. The data in the graph were processed for two tumor sizes, those above $6\mathrm{mm}$ and those below $6\mathrm{mm}$ , to illustrate the fact that the size of the tumor plays a major factor in detectability of the cancer. If the tumor size was greater than $6\mathrm{mm}$ , the NIR tomography images have reasonably high diagnostic accuracy, with the AUC equal to 0.88. However, if the tumors were smaller than $6\mathrm{mm}$ , the system would fail to distinguish normal and abnormal cases.

The analysis presented here cannot strictly be termed ROC analysis, due to the lack of a human observer detecting the presence or estimating the likelihood of an object being present; however, the quantification of the relative contrast values allowed analysis of the sensitivity-specificity tradeoff, similar to what is determined in a formal ROC curve analysis. The inability to achieve high AUC performance for all tumor sizes indicates that the imaging modality may be inadequate as a screening tool to detect small cancers. Alternatively, relatively high AUC values for larger tumors indicate that the modality is efficient at differentiating malignant versus benign tumors above $6\mathrm{mm}$ in diameter. Further analysis of diffuse tomography used for an imaging medium to large sized cancers is likely to be a promising avenue for this modality. A similar analysis for fluorescence tomography would also be beneficial.