Translator Disclaimer
1 March 2008 Highly resolved diffuse optical tomography: a systematic approach using high-pass filtering for value-preserved images
Author Affiliations +
We attempt to develop a systematic scheme through adopting high-pass filtering (HPF) to well resolve value-preserved images such as medical images. Our approach is derived from the Poisson maximum a posteriori superresolution algorithm employing the HP filters, where four filters are considered such as two low-pass-filter-combination based filters, wavelet filter, and negative-oriented Laplacian HP filter. The proposed approach is incorporated into the procedure of finite-element-method (FEM)-based image reconstruction for diffuse optical tomography in the direct current domain, posterior to each iteration without altering the original FEM modeling. This approach is justified with various HPF for different cases that breast-like phantoms embedded with two or three inclusions that imitate tumors are employed to examine the resolution performances under certain extreme conditions. The proposed approach to enhancing image resolution is evaluated for all tested cases. A qualitative investigation of reconstruction performance for each case is presented. Following this, we define a set of measures on the quantitative evaluation for a range of resolutions including separation, size, contrast, and location, thereby providing a comparable evaluation to the visual quality. The most satisfactory result is obtained by using the wavelet HP filter, and it successfully justifies our proposed scheme.



Over the last several decades, there has been great enthusiasm in developing medical imaging techniques to assist physicians in detecting and diagnosing tumors and diseases. Today, the efforts drive toward developing imaging systems employing noninvasive, nonradioactive, and relatively low cost instrumentations. Near-infrared (NIR) diffuse optical tomography (DOT) imaging is such an imaging modality that NIR light is used to probe biological tissues and it is promising to continuously monitor the status of tissues using NIR imaging. Therefore, the realization of NIR DOT as a viable clinical imaging modality would be a beneficial advancement in medical diagnosis. Basically, both the absorption and scattering tomographic images are evaluated in an NIR imaging system, thereby relating absorption properties to the oxygen saturation of hemoglobin content and water content, as well as scattering properties to the scatter size and density or the mitochondrial compartment and blood glucose concentration.1, 2, 3, 4, 5

An NIR spectral window exists from about 650to1000nm wherein the absorption is relatively small, which enables transillumination of NIR radiance through biological tissues. With a difficulty arising from strongly scattering effects in human tissues, the contrast and resolution of optical images are severely reduced. Compared with conventional x-ray mammography, magnetic resonance imaging (MRI), and ultrasound imaging all with acceptable resolutions (100mm) , but low intrinsic contrast (101) , NIR imaging possesses exceptionally high intrinsic contrast (1012) , but exhibits inferior spatial resolutions (101mm) as a result of highly scattering nature of biological tissues.4, 6 Many efforts have been made to improve NIR optical tomographic image resolution through different ways. 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34 Hebden and Deply7 proposed a method using the least-squares fit between the temporal-distribution measures of transmitted light and a model of the diffusion equation to enhance time-resolved imaging. Moon and Reintjes8 applied the Markov-chain technique to enhance optical image resolution. Jiang and Paulsen9 and Jiang 10, 11, 12 improved diffuse optical images in the direct current (dc) domain using the scheme with total variation minimization, dual mesh, and low-pass spatial filtering to achieve a satisfactory result. Recently researchers have adopted hybrid modalities to attain high-resolution NIR tomographic images by the use of a priori structural information available from MRI (Refs. 13, 14, 15, 16, 17, 18) or ultrasonic imaging.19, 20, 21 Especially, the structure information acquired from MRI was incorporated with a Laplacian-type regularization integrated in the inversion-computation process.22, 23 Additionally, the spectral priors acquired from various source wavelengths were combined with the reconstruction process, validating improvement over spatial priors.24

For the enhancement of image reconstruction, Kanmani and Vasus25 used a nonlinear approximation of the perturbation equation through adding the second term involving the Hessian in the Taylor expansion instead of a linear perturbation model that adopts only the first order derivatives (the Jacobian), which is solved by using conjugate gradient search. Furthermore, Jiang26 reconstructed optical images using the third-order diffusion equation, providing more stable inverse solutions. Pogue 27 improved diffuse optical images with spatially variant regularization in the radial orientation, thereby minimizing high-frequency noise and producing constant image resolution and contrast. Brooks 28 and Zhang 29 obtained accurate reconstruction images by the joint use of measurement-model agreement, amplitude, and total variation type constraints. Guven 30, 31 proposed an adaptive multigrid algorithm for the enhancement of image resolutions where two-level meshes were generated to provide high resolution of the region of interest. Stott 32 presented a technique to improve optical images through using simultaneous calibration of optode positions that were sensitive to image quality. Furthermore, Ntziachristos 33 and Intes 34 employed a fluorescent diagnostic agent, indocyanine green (ICG) to enhance heterogeneity contrast for obtaining better resolutions prior to optical image reconstruction. As to the background information about image processing techniques for the enhancement of reconstructed optical-property images, especially applied in this study, some monographs35, 36, 37 and related papers38, 39, 40 are valuable. Three referred to books are rather appropriate for the beginner, especially the third one, and three reference papers are actually the origin of the idea resulting in the proposed algorithm presented in this paper. Further, more references were reviewed and introduced.

In this paper, the design of a high-pass filtering (HPF) method to enhance optical images is studied. Based on the viewpoint of image processing, generally, visual quality perception is preferred to actual image values. Moreover, as is known, the effect of an HP filter applied on an image to be processed yields a different image, which cannot be used when true optical property values are required. In this paper, we attempt to develop a systematic scheme through adopting HP filters for value-preserved images such as medical images. Therefore, as simply implementing a specific HP filter to resolve NIR DOT images, it is not suitable that the procedure of the conventional approach takes routine steps like HPF the original image to be weighted and then histogram equalization. As can be understood, this conventional image-processing procedure is performed on optical property images between a reconstruction and a true distribution. For instance, one may be more interested in estimating the true values than obtaining the visual effect. As a result, an approach to systematically implementing HPF is demanding. Additionally, reconstructed highly resolved images that preserve true values are extremely expected. In this paper, first we investigate the properties of HP filters that can be classified into two types, i.e., low-pass filter (LPF) combined form and a wavelet-like filter, respectively. To preserve the true value distribution of optical images, the approach proposed and realized is derived from the Poisson maximum a priori (Poisson MAP) superresolution algorithm for the application of HPF on absorption- and diffusion-coefficient DOT images. Following this, the proposed approach is incorporated into the finite-element-method (FEM)-based image reconstruction in the continuous wave (cw) domain. Simulation results and their corresponding evaluations are demonstrated and investigated by comparing reconstruction with and without filtering.

This paper aims to (1) develop an approach derived from the Poisson MAP superresolution algorithm to systematically implement HP filters on optical property images; (2) justify the proposed approach with various HP filters performed on breast-like optical heterogeneity; (3) demonstrate the resolution performance of this approach under certain extreme conditions, significantly improving the reconstruction performance even in the absence of a priori information or modified reconstruction algorithms; and (4) define a set of measures for the evaluation of computation resolutions on the separation, size, and location of inclusions, and the contrast of inclusions to background. Additionally, further discussions on these measures are also provided. The paper is organized as follows. Section 2 briefly describes processing with HP filters that are used to enhance an image. Following this, a novel approach that starts from the Poisson MAP superresolution algorithm is derived. Section 3 implements four HP filters on several DOT images, and presents both qualitative and quantitative discussions of the reconstructed images. Finally, in Sec. 4 we draw conclusions and discuss future works.


Theoretical Analysis of the Proposed Approach

Following from the introduction this section concerns image processing. As is known, a linear image enhancement technique can improve image visual quality but cannot preserve its true values, whereas nonlinear image restoration can obtain an improved and value-preserved image, but is time-consuming. Here, we propose an approach derived from the Poisson MAP superresolution algorithm. This approach is incorporated with the procedure of FEM-based image reconstruction to obtain resolution-enhanced images. In this section, conventional image processing is first addressed, then a novel approach to implementing HP filters is derived, and finally our proposed approach is integrated with DOT image reconstruction.


Conventional Image Processing

In image processing, image enhancement is always used to improve image visual quality. The techniques of contrast enhancement, histogram equalization, and HPF are usually adopted. Contrast enhancement conducts an operation to expand the contrast of features of interest. The procedure of histogram equalization, basically, transforms the histogram distribution of an image into an output image with an equal number of pixels at each gray level. This causes a ragged histogram to become flat. HPF is exactly a transfer function with a unit at dc frequency and higher gains associated with larger frequencies. Usually, edge enhancement can be regarded as an alternative to HPF, sharpening the edge but with overshoot.

As described in some image-processing monographs,35, 36, 37 the HPF applied to improve image quality follows some routine steps such as the so-called high-frequency emphasis filtering,

Eq. 1

and histogram equalization, where Hhef denotes high-frequency emphasis filtering in the frequency domain (FD), Hhp is an HP filter, a is an offset, and b is a weighting number (usually, b> a ). Note that the capital to represent a filter is the corresponding frequency function. Therefore, 1Hlp is adopted, a complementary filter to low-pass filtering in the FD for an HP filter.

To achieve both improved image visual quality and preserve true value distribution in biological applications, we argue in this paper that an in-depth investigation is required to cope with the challenges of emergent biomedical imaging modalities. Before proposing our scheme, we first define HP, filters and classify them for the convenience of following discussions. There are two types of hhp to be performed, respectively. One is a differential filter through the combination of two LPFs, and the other is an intrinsic (wavelet-like) HP filter. It is sensible that an HP filter can be described as

Eq. 2

where hlp1 and hlp2 denote LPFs. More precisely, this represents a narrow full width at half maximum (FWHM) LPF with a larger amplitude (A1) subtracted by a broad FWHM LPF with a smaller amplitude (A2) . Both HP filters must comply with two rules of thumb37 as follows:

Eq. 3

The difference between these filters is that the hlp1 of the former rolls off faster than that of the latter. In this study, Gaussian functions with various standard deviations (σ) are employed for hlp . Thus, the HP filter shown in Fig. 1 can be formulated as

Eq. 4

where g1(r)=[A1(2πσ12)12]exp(r22σ12) and g2(r)=[A2(2πσ22)12]exp(r22σ22) , respectively, and σ1<σ2 .

Fig. 1

Two LPF-combined HP filters (a) δg2 and (b) g1g2 , and (c) a wavelet-like filter, where δ is the delta function.



LPF-combined HP filter (σ1<σ2)

This filter is usually determined with a smaller σ1(1) , as shown in Fig. 1a; this can also be determined with a larger σ1 , as shown in Fig. 1b. If we let σ1 approach zero, hhp1 narrows further to an impulse, then Eq. 4 can be expressed in the frequency domain as 1Hlp .


Wavelet-like HP filter

A dilated wavelet-like function expressed as

Eq. 5

where a is a dilated factor, as depicted in Fig. 1c, and can be used as an HP filter.


Negative-oriented Laplacian HP filter

Alternatively, a 3×3 negative-orientated Laplacian edge operator

Eq. 6

in a form similar to a wavelet is also considered and employed as an HP filter is this study.

As can be seen in Fig. 1, a wavelet-like HP filter has sharp sidelobes rather than a LPF-combined HP filter.


Novel Approach Implementing HPF for Optical Tomography

To find the fundamental theoretical basis to explain our proposed novel approach, an attempt to derive the value-preserving HPF technique begins with the Poisson MAP superresolution algorithm. Mathematically, the algorithm38, 39, 40 is given as

Eq. 7


Eq. 8

is regarded as the correction term during the iterative restoration progress; ⊗ represents a convolution; * represents correlation; h denotes the point spread function (PSF); g is the observed image; and the subscript n is the number of iteration. Additionally, f̂0 defined as g is the initial guess of iteration, and f̂N is the final superresolved image. In terms of the operation of Poisson MAP, it is an iterative algorithm, where successive estimate of the restored image is obtained through the multiplication of current estimate by such a quantity close to one that is a function of the interpolated image divided by a convolution of the current estimate with the PSF. Using Taylor series expansion, Eq. 7 can be expanded to be approximate as

Eq. 9

where the intermediate obtained image f̂n can be expressed as adding the previous one f̂n1 with a correction increment Δf̂n1 . It is nontrivial to further explore the correction increment

Eq. 10

where gf̂n1h can be reduced to an increment Δf̂n1 in a condition of decreasing correction rate. Assuming f̂n1h approaches a constant as f̂n1 has a simple distribution. Thus, Eq. 10 is approximated to

Eq. 11

Through point-by-point multiplying both sides of Eq. 11 with Δf̂n1 , we obtain

Eq. 12

and then reorganize Eq. 12 to yield

Eq. 13

As defined previously, h is the PSF like an LPF, and a Gaussian function is employed in the study. Additionally, the operation of correlation is equivalent to take a convolution due to the symmetry of function h . Therefore, Eq. 13 can be derived to the following equation with the HPF definition

Eq. 14

We here consider the quantity obtained through the convolution of an image and HPF as a new correction increment, i.e.,

Eq. 15

Thus, Eq. 15 is equivalent to

Eq. 16

To further simplify Eq. 16 for numerical evaluation, we assume that the denominator is a positive number relative to Δf̂n1 , and finally get an approximate solution of the correction increment for using HPF, as follows

Eq. 17

where the denominator is simplified with the norm of Δf̂n1 multiplied by w , a weight number 10, used in computation as well as herein the symbols ⟨∣ and ∣⟩ stand for the state and ⟨ zx and xy ⟩ represent the operations of a point-by-point product ( z and x ) and a convolution ( x and y ), respectively, resulting in the other state. In considering Eq. 17 used in the computation of NIR DOT imaging, heterogeneities are treated as a perturbation to homogeneous background for a phantom, and incremental values of both absorption and scattering coefficients are estimated from a projection of a high-frequency enhancement to the original increment.


NIR DOT Image Reconstruction Incorporating with Novel Approach

Compared with other medical imaging modalities, NIR imaging requires the solution of an inverse problem. In NIR DOT imaging, the fundamental equation governing the propagation of light in biological tissues is the Boltzmann transport equation (BTE) to model the optical characteristics of the scattering and absorption.

The BTE is an integrodifferential equation, so it is rather difficult to obtain solutions to the BTE under general conditions. With the use of approximation techniques by assuming the experimental material or tissues have highly scattering properties and that the input radiance is isotropic and modulated under a 1-GHz frequency, the BTE can be reduced to an easily solvable form of the diffusion approximation. In NIR imaging, mappings of the absorption and/or scattering coefficients can be evaluated by using an FEM to invert the diffusion approximation. The FEM-based image reconstruction in the dc domain is concluded with the following equations. More derivation details can be found in Ref. 41.

As described previously, the physical process can be deduced from a diffusion equation:

Eq. 18

where S(r,ω) and Φ denote the source and the radiance, respectively; and μa , c , and D are the absorption coefficient, the wave speed in the medium, and the diffusion coefficient, respectively. To solve Eq. 18, the boundary condition DΦn̂=αΦ (flux in fact) and the FEM are applied. Since only dc data are considered, ω is set as a null; i.e., the imaginary part is to vanish from the subsequent equations. Thus, the following discrete equations in a matrix form,

Eq. 19

can be obtained. Obviously, the forward solution, Φ , can be evaluated through Eq. 19. In terms of the physical process, the radiance matrix is quantitatively and qualitatively dependent on the source matrix and the optical-property matrix, respectively, where the optical-property matrix is the inertia of the material in spite of relating to the wavelength. Furthermore, the following two equations can be derived for the computation of image reconstruction, i.e.,

Eq. 20


Eq. 21

where the superscripts I and b denote interior and boundary nodes, and Dk for k=1,2,,K and μl for l=1,2,,L are the reconstruction parameters for the optical property profile. For the inverse problem to update absorption/scattering coefficients, the partial differentiation of boundary radiance to the parameters of interest, Φbμl or ΦbDk , must be obtained from Eqs. 20, 21. The Newton-Raphson technique regularized by a Levenberg-Marquardt algorithm and with the Tikhonov regularization parameter is adopted to iteratively update the diffusion and absorption coefficients, i.e.,

Eq. 22

where Jacobian matrix J denotes J(ΦbDk,Φbμl) , Δχ means Δχ(ΔDk,Δμl) , and λ is a Tikhonov regularization parameter of the Jacobian matrix. As described, this inversion generally requires the construction of the Jacobian matrix; actually, the Jacobian represents a highly underdetermined system of equations. Although it is possible to obtain a least-squares solution to underdetermined systems of equations, the resulting images are usually inaccurate relating to inferior resolution. The procedure of FEM-based image reconstruction in the dc domain is illustrated in Fig. 2 . As indicated, the proposed approach is merely implemented once, subsequently posterior to each iteration without altering the original FEM modeling. The following section illustrates the comparison and effectiveness of the incorporated resolution-enhanced schemes.

Fig. 2

Flowchart of NIR DOT image reconstruction incorporated with the HPF approach.



Results and Discussion

The phantoms employed to justify our proposed technique incorporate two or three inclusions with various sizes, locations, and separations, illustrated in Fig. 3 , where R denotes the radius in millimeters. In this paper, four HPFs and four phantom cases were performed. The numerical simulations of multiinclusion phantoms provide further information concerning the spatial resolution (separation, size, and location) and the contrast resolution beyond that of the single-inclusion case. Of the phantom, the background absorption (μa) and reduced scattering (μs) values are about 0.0025 and 0.25mm1 , respectively, while the maximum absorption and reduced scattering for the inclusion are 0.025 and 2.5mm1 , if we assume the contrast ratio of the inclusion to background is 10:1, because high contrast results in much more overlapping effects than low contrast, although a contrast of 2 to 10 was used throughout other published works.

Fig. 3

Schematic diagram for the dimensions of four different test cases in simulation. (a) to (d) are cases 1 to 4, respectively, where R is radius in millimeters.


As depicted in Fig. 3, cases 1 and 2 and cases 3 and 4, respectively, have three and two inclusions separated by a similar distance but of different sizes. As the separation resolution of inclusions is examined, several (two or three) embedded inclusions are necessary, and different inclusion sizes are considered as well. To test the limitation of each HPF employed here, the phantoms of cases 1 and 4 with larger inclusions and closer to the phantom center were designed, compared with case 2 and case 3 designs. For the convenience of discussion, we denote M0 to M4 as the reconstructions with the schemes using nonfiltering, δg2(σ2=1.5) , g1g2(σ1=0.75,σ2=1.5) , wavelet (a=0.5) , and Laplacian HPF in their 2-D form, respectively. Currently, absorption-coefficient images are presented for our cw image reconstruction algorithm.

In FEM-based image reconstruction, the homogeneous background ( μa=0.0025mm1 , μs=0.25mm1 ) was adopted as an initial guess. For both the forward and inverse processes, 256 elements and 257 nodes were used, associated with a desktop PC with a 3.6-GHz CPU and 4Gbytes of RAM, respectively. Thirty iteration assignments were employed for each case as the normalized increasing rate, i.e., mean value of (Φn+1Φn)Φn2 , reaches less than 102 , where each iteration takes about 2min . Meanwhile, the absorption- and diffusion-coefficient images were updated concurrently in spite of the fact that reconstruction began from a homogeneous condition and only the acquired dc data were employed.

First, a qualitative investigation of the reconstruction performance of each case is presented in Sec. 3.1. Following this in Sec. 3.2, we describe quantitative performance measures for various HPFs for a range of resolutions including separation, size, contrast, and location. Finally, we further investigate and discuss the significance of the proposed measures in Sec. 3.3.


Examples Illustration


Case 1

Figure 4 shows a set of reconstructed absorption-coefficient images [Figs. 4a, 4b, 4c, 4d, 4e] and quantitative information [Figs. 4f, 4g, 4h, 4i, 4j] for the images along with their corresponding circular transaction profiles. Comparing the reconstructed absorption images, it is obvious that all of the reconstructions show roughly correct images of the inclusions and all of the reconstruction techniques can highly resolve images and separate inclusions, except the result using M0. However, the M1 to M4 schemes underestimate the computed absorption coefficients of the inclusions. For further inspection, M4 generated highly ringing artifacts between inclusions. At this phase, it is not easy to speculate about the causes of such artifacts that might be referred to as ‘false’ inclusions and concluded as a wrong judgment.

Fig. 4

Case 1, 2-D reconstructed absorption images (a) without HPF (M0) and (b) to (e) with M1, M2, M3, M4 filtering, respectively; (f) to (j) 1-D sectional profiles corresponding to (a) to (e), where the solid lines are the designed and the dotted lines represent the reconstructed schemes.



Case 2

It can be found that compared with only the image reconstruction employed, considerable improvement is observed in the reconstructed images, as illustrated in Figs. 5a, 5b, 5c, 5d, 5e when the HPF approach is invoked. Evidently, M0 retains highly blurred inclusions, while the other reconstruction schemes can better differentiate inclusions, and the M4 scheme overestimated the absorption coefficients. Again, ringing artifacts are produced surrounding inclusions and their optical properties are lower than background levels, as depicted in Figs. 5f, 5g, 5h, 5i, 5j.

Fig. 5

Reconstructed case 2 images, with (a) to (j) as described for Fig. 4.


From the results for cases 1 and 2 note that schemes with filtering can discriminate even small size inclusions, whereas scheme M0 cannot meet even the basic requirements of image reconstruction, especially for small inclusions (case 2).


Case 3

Compared with previous two cases, this case was designed as a phantom with three smaller inclusions. Several improved images were obtained by using appropriate filtering, as shown in Figs. 6b, 6c, 6d, 6e . Likewise, M2 resulted in a worse-resolved image than the others with HP filtering. Negative artifacts occurred in each reconstructed image, as depicted in Figs. 6g, 6h, 6i, 6j. It is well noted that M4 overestimated the inclusion amplitudes, which yields a higher inclusion-to-background contrast.

Fig. 6

Reconstructed case 3 images, with (a) to (j) as described for Fig. 4.



Case 4

In this highly challenging case, a phantom with two closest-separation inclusions was designed. As shown in Figs. 7a, 7b, 7c, 7d, 7e , all reconstructed images underestimated inclusions, and offered relatively poor resolution for two separate inclusions. This is rather competitive for these employed filters. Based on a quantitative comparison, as depicted in Figs. 7i and 7j, the M3 and M4 schemes demonstrate better resolution discrimination to separate longer and closer inclusions in comparison with case 3.

Fig. 7

Reconstructed case 4 images, with (a) to (j) as described for Fig. 4.


From the results of cases 3 and 4 for a phantom with inclusions of both small size and close separation, it can be concluded that the wavelet-like HP filtering (M3) demonstrates the best spatial-resolution capability to the inclusions.

This evidently shows that the enhancement of reconstruction through the incorporation of our proposed HPF approach can effectively improve computed images. As already illustrated, the wavelet-like HP filtering schemes (M3 and M4) further yield better results than the LPF-combined HP filtering schemes (M1 and M2). In the aspects of sensitivity and stability of evaluation, the M3 scheme yielded results closest to the true absorption property compared to the other schemes. However, scheme M4 visually characterizes the inclusion-to-background contrast best.


Performance Investigation

In terms of the optical properties within the inclusion and the background, note that the image reconstruction not only pursues qualitative correctness but also obtains favorably quantitative information about the optical properties of either the inclusions or the background. The parameters of interest, such as size, contrast, and location variations associated with image quantification measures are most frequently investigated and discussed.

Several measures42, 43, 44, 45 have been used to evaluate the performance of the NIR imaging algorithms or systems. Song 42 used the contrast-to-noise ratio (CNR), which is defined as the difference between the region of interest (ROI) and the background region values of the optical properties divided by the average variation in the background, where one inclusion was considered. Furthermore, informative works (Pogue 44) provided an overview of the three major methods utilized for image analysis in the imaging science and medical physics communities, which lie in the areas of the spatial resolution, the contrast detail (CD) analysis, and human perception of images. Briefly, the first one relates to the modulation transfer function (MTF) profile; the second one, to the CD curve (contrast versus size) obtained by human observation; and the last concerns the receiver operating characteristic (ROC) curve and location receiver operating characteristic curve (LROC) obtained by the human observer detection of abnormalities. In our cases, however, several inclusions in the background of a phantom were considered and further investigations of the contrast, size, separation, and location were conducted so that it is essential that these four terms are respectively defined and discussed.

To provide a quantitative assessment for these reconstructed images through using various HPF approaches, we designate and formulate four measures over the ROI for the evaluation of these filtering schemes, and these measures are normalized to be in-between a null and unit with the ratio of the reconstructed to the original images. To interpret these measures in detail, we describe them using Fig. 8 , where 1ḎROI is chosen as the line segment between the two outmost nodes of the inclusions, and 2ḎROI is the possibly smallest region around and/or covering inclusions. Moreover, although it is usually a difficult task to define a separation resolution, here we regard an inclusion and a separation as a bump and a cave, respectively.

Fig. 8

Diagram for the explanation of defined measures.



Contrast resolution (Rcont.1D,2D)

The measure Rcont.1D,2D is defined to evaluate the resolution on the contrast of optical property values with the inclusions to the background region especially between inclusions.

Eq. 23


Eq. 24

where max¯ and min¯ denote the average of maxima and minima over all the selected regions as the superscripts ( Δincl. and Δincl¯ ). Also incl. and incl¯ correspondingly represent inclusions and complementary inclusions as well as Δincl. and Δincl¯ are chosen with several nodes around central area of incl. and incl.


Separation resolution (Rsep1D,2D)

The measure Rsep1D,2D is designed to evaluate the resolution on the separation between inclusions.

Eq. 25

where MSE is the mean square error over all the selected region as the superscript (incl¯) and Baseline is used with 0.025mm1 ; Ori.2Baseline is, here, 0.0025to0.025mm1 , and Recon.2Ori. is the reconstructed value in the region incl¯ to 0.0025mm1 .


Size resolution (Rsize1D,2D)

The measure Rsize1D,2D is designed to evaluate the resolution on the size over all inclusions.

Eq. 26

where MSE is over the selected region incl. and baseline is used with 0.0025mm1 . Note that the FWHM usually operated manually and subjectively is not adopted for the evaluation of inclusion size. Here, attempt to automatically estimate this resolution with the idea of the capacity rate for the term of interest.


Location resolution (Rlocat1D,2D)

The measure Rlocat1D,2D is defined to evaluate the resolution on the location over all inclusions.

Eq. 27


Eq. 28

where CM¯ is the average of the center of mass over all the selected region as the superscript (incl.).

Note that the resolutions, Rsep1D,2D , Rsize1D,2D , and Rlocal1D,2D include a multiplication operation by Rcont.1D,2D to avoid the low-contrast reconstruction with high-consistent inclusions or complements. Furthermore, it is expected that the resolution is higher because the R value is approaching more closely to a unit.

Based on the preceding definitions, the evaluated resolutions of 1-D profiles and 2-D images with multiinclusions are listed in Tables 1, 2, 3, 4 for each phantom case, respectively. The quantities for Rsep1D,2D , Rsize1D,2D , and Rlocal1D,2D are small because the defined measures are quite strict. For overall cases, it is found that location resolution is above 0.95 prior to Rolocat1D,2D multiplied by Rcont.1D,2D , and less difference exists between these two; whereas the contrast, separation, and size resolution have comparable differences. Figures 9, 10, 11, 12 illustrate comparisons of the separation, size, and contrast resolutions among various HP filters to clarify our observation. Overall, the resolutions obtained in cases 1 to 3 are better than those in case 4, as expected. Basically, our approach demonstrates the effectiveness of separation and size resolution rather than contrast resolution. A discussion of each individual case follows.

Fig. 9

Case 1: (a) 1-D measures and (b) 2-D measures, where solid, +, ⋆, △, and ○ lines represent using schemes of M0 to M4, respectively.


Fig. 10

Case 2 with (a), (b) and line key as described for Fig. 9.


Fig. 11

Case 3 with (a), (b) and line key as described for Fig. 9.


Fig. 12

Case 4 with (a), (b) and line key as described for Fig. 9.


Table 1

Case 1 separation, size, contrast, and location resolutions for various filtering on 1-D and 2-D conditions.

1-D 2-D
Sep. 0Sep.Size 0SizeContrastLoc. 0Loc.Sep. 0Sep.Size 0SizeContrast Locx . 0 Locx . Locy . 0 Locy .

Table 2

Case 2 separation, size, contrast, and location resolutions for various filtering on 1-D and 2-D conditions.

1-D 2-D
Sep. 0Sep.Size 0SizeContrastLoc. 0Loc.Sep. 0Sep.Size 0SizeContrast Locx . 0 Locx . Locy . 0 Locy .
M40.870.870.401.000.74 0.07 0.501.000.98
0.590.590.630.61 0.19 0.710.70

Table 3

Case 3 separation, size, contrast, and location resolutions for various filtering on 1-D and 2-D conditions.

1-D 2-D
Sep. 0Sep.Size 0SizeContrastLoc. 0Loc.Sep. 0Sep.Size 0SizeContrast Locx . 0 Locx . Locy . 0 Locy .

Table 4

Case 4 separation, size, contrast, and location resolutions for various filtering on 1-D and 2-D conditions.

1-D 2-D
Sep. 0Sep.Size 0SizeContrastLoc. 0Loc.Sep. 0Sep.Size 0SizeContrast Locx . 0 Locx . Locy . 0 Locy .


Case 1

Figure 9a shows the resolution performance of schemes M3, M1, M2, M4, and M0, respectively. The results show that the M4 scheme yielded false inclusions. For the 2-D condition, the revealed performance is similar to that in the 1-D condition except for the M0 and M4 schemes. Obviously, this evaluation is consistent with that for Fig. 4 based on the visual perception.


Case 2

Figure 10a shows the performance of schemes M4, M1, M3, M2, and M0, respectively. The performance in the 2-D condition is similar to that in the 1-D condition except for the M4 scheme. Unfortunately, a negative value occurs in the 2-D condition [Fig. 10b], which means that M4 highly overestimated the inclusion size.


Case 3

Obviously, the performance of case 3 is similar to that of case 2, shown as Fig. 11. For the same reasons as in case 2, the highly overestimated effect of M4 attenuates the measure values. Basically, other filtering schemes obtain the measure distribution as expected.

Generally speaking, the M1 and M3 schemes perform better than the M2 and M4 schemes on either of the defined measures or the visual perception for cases 2 and 3.


Case 4

In this case, Fig. 7 shows that only wavelet-like HP filtering is able to resolve images well. As expected, Fig. 12 shows a better performance of schemes M4 and M3 than that of schemes M1, M2, and M0 on both the 1-D and the 2-D measures.

In summary, it can be seen that the evaluations depicted in Figs. 9, 10, 11, 12 using our defined measures are quite consistent with those evaluations based on visual perception on Figs. 4, 5, 6, 7.

Case 1 is the only example that is resolved to some extent without having to use HPF (M0). However, scheme M0 made some measure evaluation better than others since the corresponding M0 reconstructions have a nearly uniform distribution. In spite of this, the measures we defined remain effective for most of the 1-D and 2-D cases.


Evaluation on Defined Measures

In an aspect of quantitative discussions on resolution, we employ these four measures to explain the effectiveness of each proposed filtering. Particularly, accurate demonstrations for the 1-D condition are almost fully matched with the evaluation in quality. To an extent, most are also promising for the 2-D condition. In other words, the evaluation implies that our defined measures are quite acceptable. For further inspection, these measures can be seen based on individual inclusion or separation as well. In this subsection, more discussion of the defined measures is given. First, the contrast resolution can be also defined as

Eq. 29

where mean is to find the average value of the selected regions as the superscripts (incl. and incl¯ ). A definition in this manner, however, is not suitable for our cases 1 to 3. For further investigation, Eq. 23 can be concluded as

Eq. 30

Rcont.1D,2D{> cnormalsituation=cnocontrast<cabnormalsituation},
where c is equal to 1(μincl.μincl.) (0.1 is used here) and the abnormal situation, here, means the optical value of the inclusion is smaller than that of the separation region. Likewise, the separation and size resolution can be defined as

Eq. 31

It is found that Eq. 31 eventually regards a reconstructed “inclusion” as a reverse cave with values ranging between 0 and 1, and Eq. 31 always gives positive values. When considering Eqs. 25, 26, it can be proven that both Rosep.;size1D,2D are always smaller than a unit and, moreover, are negative values to denote a high underestimation or overestimation. Note that our 2-D ROI can be determined automatically using a computer program but not with manual selection, whereas the FWHM is not adopted in this paper because it is selected manually. Finally, the location resolution is regulated by Eq. 28. Prior to adjustment, the positive or negative errors to the unit can be explained as denoting the multiinclusion position in a reverse direction.


Concluding Remarks

We proposed and implemented a resolution-enhancing technique using HPF incorporated with the FEM-based inverse computation to obtain highly resolved NIR diffuse optical images in a systematical manner. As mentioned previously, our approach, derived from the Poisson MAP, was justified by various HPFs for different designated phantoms. Qualitative visual perception and quantitative evaluations of the reconstructions also validate the proposed approaches.

Obviously, the wavelet-like HP filtering is superior to the LPF-combined HPF, as shown in Figs. 4, 5, 6, 7. In summary, the approach to use the wavelet-like HP filtering, M3, is recommended in terms of its resolving ability and computational stability. It is observed that the M4 scheme demonstrates a high resolution result as well, but reveals worse stability than the M3 scheme. Additionally, a small inclusion-to-background diameter ratio, 2:20, is detectable and distinguished.

Due to the variation in the choice of σ1 and σ2 associated with each filter, various filters result in different reconstruction results. In this paper, we did not attempt to conduct a wide comparison and an extensive study over a range of HP filters and phantom cases, but rather chose to begin with two categories of filters and a set of more-or-less extreme cases. Although the resolutions of absorption images enhanced with our proposed techniques were presented, this approach remains effective to improve the scattering images for the frequency-domain DOT imaging system as well. In future work, a thorough investigation of HP filters used in the proposed approach will be conducted to find one or several appropriate filters. Moreover, the resolution limit should be specified with a set of designed cases. A further study is also required to identify exact causes of the negative-value artifacts shown in Figs. 5f, 5g, 5h, 5i, 5j, 6f, 6g, 6h, 6i, 6j which are dependent either on the filters or on the cases themselves. Owing to the lack of a sound method for quantitative evaluation, it is believed that even to objectively define a measure corresponding to visual perception is quite complicated. Alternatively, four reasonable measure definitions were considered and defined to provide an initial basis for quantitative evaluations, from which further explorations of an individual inclusion or separation ROI can begin. Briefly, our proposed measures mainly provide approximate information, and areas remain for further investigation and improvement.


The authors would like to acknowledge the funding support from the grants by the Veteran General Hospital/University System of Taiwan Joint Research Program (VGHUST95-P4-13, VGHUST96-P4-17) and the National Science Council (NSC 95-2221-E-236-002) in Taiwan.



G. A. Millikan, “The oximeter, an instrument for measuring continuously the oxygen saturation of arterial blood in man,” Rev. Sci. Instrum., 13 434 –444 (1942). 0034-6748 Google Scholar


J. S. Maier, S. A. Walker, S. Fantini, M. A. Franceschini, and E. Gratton, “Possible correlation between blood glucose concentration and the reduced scattering coefficient of tissues in the near infrared,” Opt. Lett., 19 2062 –2064 (1994). 0146-9592 Google Scholar


M. Kohl, M. Cope, M. Essenpreis, and D. Böcker, “Influence of glucose concentration on light scattering in tissue-simulating phantoms,” Opt. Lett., 19 2170 –2172 (1994). 0146-9592 Google Scholar


B. W. Pogue, S. D. Jiang, H. Dehghani, C. Kogel, S. Soho, S. Srinivasan, X. M. Song, T. D. Tosteson, S. P. Poplack, and K. D. Paulsen, “Characterization of hemoglobin, water, and NIR scattering in breast tissue: analysis of intersubject variability and menstrual cycle changes,” J. Biomed. Opt., 9 (3), 541 –552 (2004). 1083-3668 Google Scholar


B. Brooksby, S. D. Jiang, H. Dehghani, B. W. Pogue, K. D. Paulsen, J. Weaver, C. Kogel, and S. P. Poplack, “Combining near-infrared tomography and magnetic resonance imaging to study in vivo breast tissue: implementation of a Laplacian-type regularization to incorporate magnetic resonance structure,” J. Biomed. Opt., 10 051504 (2005). 1083-3668 Google Scholar


P. J. Cassidy and G. K. Radda, “Molecular imaging perspectives,” J. R. Soc., Interface, 2 133 –144 (2005). 1742-5689 Google Scholar


J. C. Hebden and D. T. Deply, “Enhanced time-resolved imaging with a diffusion model of photon transport,” Opt. Lett., 19 (5), 311 –313 (1994). 0146-9592 Google Scholar


J. A. Moon and J. Reintjes, “Image resolution by use of multiply scattered light,” Opt. Lett., 19 (8), 521 –523 (1994). 0146-9592 Google Scholar


H. Jiang and K. D. Paulsen, “A finite element based higher-order diffusion approximation of light propagation in tissues,” Proc. SPIE, 2389 608 –614 (1995). 0277-786X Google Scholar


H. Jiang, K. D. Paulsen, and U. L. Österberg, “Optical image reconstruction using DC data: simulations and experiments,” Phys. Med. Biol., 41 1483 –1498 (1996). 0031-9155 Google Scholar


H. Jiang, K. D. Paulsen, U. L. Österberg, and M. S. Patterson, “Frequency-domain optical image reconstruction in turbid media: an experimental study of single-target detectability,” Appl. Opt., 36 52 –63 (1997). 0003-6935 Google Scholar


H. Jiang, K. D. Paulsen, U. L. Österberg, and M. S. Patterson, “Improved continuous light diffusion imaging in single- and multi-target tissue-like phantoms,” Phys. Med. Biol., 43 675 –693 (1998). 0031-9155 Google Scholar


B. W. Pogue and K. D. Paulsen, “High-resolution near-infrared tomographic imaging simulations of the rat cranium by use of a priori magnetic resonance imaging structural information,” Opt. Lett., 23 1716 –1718 (1998). 0146-9592 Google Scholar


B. W. Pogue, T. McBride, C. Nwaigwe, U. L. Österberg, J. F. Dunn, and K. D. Paulsen, “Near-infrared diffuse tomography with a priori MRI structural information: testing a hybrid image reconstruction methodology with functional imaging of the rat cranium,” Proc. SPIE, 3597 484 –492 (1999). 0277-786X Google Scholar


V. Ntziachristos, “Concurrent diffuse optical tomography, spectroscopy and magnetic resonance imaging of breast cancer,” University of Pennsylvania, (2000). Google Scholar


H. Dehghani, B. W. Pogue, and K. D. Paulsen, “Development of hybrid NIR/MRI imaging system algorithm: use of a-priori information for tumor detection in the female breast,” 657 –660 (2002). Google Scholar


H. Xu, H. Dehghani, B. W. Pogue, K. D. Paulsen, and J. F. Dunn, “Hybrid MR/near infrared imaging of the murine brain: optimization of optical fiber arrangement and use of a-priori knowledge,” 74 –77 (2002). Google Scholar


B. A. Brooksby, H. Dehghani, B. W. Pogue, and K. D. Paulsen, “Near-infrared (NIR) tomography breast image reconstruction with a priori structural information from MRI: algorithm development for reconstructing heterogeneities,” IEEE J. Sel. Top. Quantum Electron., 9 199 –209 (2003). 1077-260X Google Scholar


Q. Zhu, T. Durduran, V. Ntziachristos, M. Holboke, and A. G. Yodh, “Imager that combines near-infrared diffusive light and ultrasound,” Opt. Lett., 24 1050 –1052 (1999). 0146-9592 Google Scholar


P. Guo, D. Piao, Q. Zhu, and J. Fikiet, “A combined 2-D ultrasound and NIR imaging system,” 77 –78 (2000). Google Scholar


M. J. Holboke, B. J. Tromberg, X. Li, N. Shah, J. Fishkin, D. Kidney, J. Butler, B. Chance, and A. G. Yodh, “Three-dimensional diffuse optical mammography with ultrasound localization in a human subject,” J. Biomed. Opt., 5 237 –247 (2000). 1083-3668 Google Scholar


B. Brooksby, S. Jiang, H. Dehghani, B. W. Pogue, K. D. Paulsen, J. Weaver, C. Kogel, and S. P. Poplack, “Combining near-infrared tomography and magnetic resonance imaging to study in vivo breast tissue: implementation of a Laplacian-type regularization to incorporate magnetic resonance structure,” J. Biomed. Opt., 10 051504 (2005). 1083-3668 Google Scholar


B. Brooksby, B. W. Pogue, S. Jiang, H. Dehghani, S. Srinivasan, C. Kogel, T. D. Tosteson, J. Weaver, S. P. Poplack, and K. D. Paulsen, “Imaging breast adipose and fibroglandular tissue molecular signatures by using hybrid MRI-guided near-infrared spectral tomography,” Proc. Natl. Acad. Sci. U.S.A., 103 (23), 8828 –8833 (2006). 0027-8424 Google Scholar


B. Brooksby, S. Srinivasan, S. Jiang, H. Dehghani, B. W. Pogue, and K. D. Paulsen, “Spectral priors improve near-infrared diffuse tomography more than spatial priors,” Opt. Lett., 30 (15), 1968 –1970 (2005). 0146-9592 Google Scholar


B. Kanmani and R. M. Vasu, “Diffuse optical tomography through solving a system of quadratic equations without re-estimating the derivatives: the Frozen-Newton method,” Proc. IEEE Int. Workshop on Biomedical Circuits and Systems, S2.2-17 –20 (2004) Google Scholar


H. Jiang, “Optical image reconstruction based on the third-order diffusion equations,” Opt. Express, 4 241 –246 (1999). 1094-4087 Google Scholar


B. W. Pogue, T. O. McBride, J. Prewitt, U. L. Österberg, and K. D. Paulsen, “Spatially variant regularization improves diffuse optical tomography,” Appl. Opt., 38 2950 –2961 (1999). 0003-6935 Google Scholar


D. H. Brooks, R. J. Gaudette, E. L. Miller, C. A. DiMarzio, D. Boas, and M. Kilmer, “An admissible solution approach for diffuse optical tomography,” 333 –337 (2000). Google Scholar


Y. H. Zhang, D. H. Brooks, and D. Boas, “A multi-resolution admissible solution approach for diffuse optical tomography,” 1005 –1008 (2002). Google Scholar


M. Guven, B. Yazici, X. Intes, and B. Chance, “An adaptive multigrid algorithm for region of interest diffuse optical tomography,” 823 –826 (2003). Google Scholar


M. Guven, B. Yzazici, X. Intes, and B. Chance, “An adaptive v-grid algorithm for diffuse optical tomography,” 95 –96 (2003). Google Scholar


J. J. Stott, J. P. Culver, S. R. Arridge, and D. A. Boas, “Optode positional calibration in diffuse optical tomography,” Appl. Opt., 42 3154 –3162 (2003). 0003-6935 Google Scholar


V. Ntziachristos, A. G. Yodh, M. Schnall, and B. Chance, “Concurrent MRI and diffuse optical tomography of breast after indocyanine green enhancement,” Proc. Natl. Acad. Sci. U.S.A., 97 2767 –2772 (2000). 0027-8424 Google Scholar


X. Intes, J. Ripoll, Y. Chen, S. Nioka, A. G. Yodh, and B. Chance, “In vivo continuous-wave optical breast imaging enhanced with indocyanine green,” Med. Phys., 30 1039 –1047 (2003). 0094-2405 Google Scholar


W. K. Pratt, Digital Image Processing, Wiley, New York (1991). Google Scholar


R. J. Schalkoff, Digital Image Processing and Computer Vision, Wiley, New York (1989). Google Scholar


R. C. Gonzalez, R. E. Woods, and S. L. Eddins, Digital Image Processing, Prentice Hall, Upper Saddle River, NJ (2004). Google Scholar


B. R. Hunt and P. J. Sementilli, “Description of a Poisson imagery super-resolution algorithm,” 196 –199 (1992). Google Scholar


M.-C. Pan, “Improving a single down-sampled image using probability-filtering-based interpolation and improved Poisson maximum a posteriori super-resolution,” EURASIP J. Appl. Signal Process., 2006 97492 (2006). 1110-8657 Google Scholar


M.-C. Pan, “A novel blind super-resolution algorithm for restoring Gaussian blurred images,” Int. J. Imaging Syst. Technol., 12 (6), 239 –246 (2002). 0899-9457 Google Scholar


K. D. Paulsen and H. Jiang, “Spatially varying optical property reconstruction using a finite element diffusion equation approximation,” Med. Phys., 22 691 –701 (1995). 0094-2405 Google Scholar


X. Song, B. W. Pogue, S. Jiang, M. M. Doyley, H. Dehghani, T. D. Tosteson, and K. D. Paulsen, “Automated region detection based on the contrast-to-noise ratio in near-infrared tomography,” Appl. Opt., 43 (5), 1053 –1062 (2004). 0003-6935 Google Scholar


S. C. Davis, B. W. Pogue, H. Dehghani, and K. D. Paulsen, “Contrast-detail analysis characterizes diffuse optical fluorescence tomography image reconstruction,” J. Biomed. Opt., 10 (5), 050501-1 –3 (2005). 1083-3668 Google Scholar


B. W. Pogue, S. C. Davis, X. Song, B. A. Brooksby, H. Dehghani, and K. D. Paulsen, “Image analysis methods for diffuse optical tomography,” J. Biomed. Opt., 11 (3), 033001-1 –16 (2006). 1083-3668 Google Scholar


C. A. Kelsey, R. D. Moseley Jr., J. F. Garcia, F. A. Mettler Jr., T. W. Parker, and J. H. Juhl, “ROC and contrast detail image evaluation tests compared,” Radiology, 154 (3), 629 –631 (1985). 0033-8419 Google Scholar
©(2008) Society of Photo-Optical Instrumentation Engineers (SPIE)
Min-Cheng Pan, Chien-Hung Chen, Liang-Yu Chen, Min-Chun Pan D.V.M., and Yi-Ming Shyr "Highly resolved diffuse optical tomography: a systematic approach using high-pass filtering for value-preserved images," Journal of Biomedical Optics 13(2), 024022 (1 March 2008).
Published: 1 March 2008

Back to Top