1.

Introduction

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 $650\text{to}1000\mathrm{nm}$ 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 $(\sim {10}^0\mathrm{mm})$ , but low intrinsic contrast $(\sim {10}^{-1})$ , NIR imaging possesses exceptionally high intrinsic contrast $(\sim {10}^{1–2})$ , but exhibits inferior spatial resolutions $(\sim {10}^1\mathrm{mm})$ 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 Deply⁷ 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 Reintjes⁸ applied the Markov-chain technique to enhance optical image resolution. Jiang and Paulsen⁹ 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.²⁴

For the enhancement of image reconstruction, Kanmani and Vasus²⁵ 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, Jiang²⁶ reconstructed optical images using the third-order diffusion equation, providing more stable inverse solutions. Pogue ²⁷ 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 ²⁸ and Zhang ²⁹ 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 ³² presented a technique to improve optical images through using simultaneous calibration of optode positions that were sensitive to image quality. Furthermore, Ntziachristos ³³ and Intes ³⁴ 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 monographs^{35, 36, 37} and related papers^{38, 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.

2.

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.

2.1.

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

$$H_{\mathrm{hef}}=a+b{H}_{\mathrm{hp}},$$

and histogram equalization, where $H_{\mathrm{hef}}$ denotes high-frequency emphasis filtering in the frequency domain (FD), $H_{\mathrm{hp}}$ 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, $1-H_{\mathrm{lp}}$ 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 $h_{\mathrm{hp}}$ 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

$$h_{\mathrm{hp}}=h_{\mathrm{lp}1}-h_{\mathrm{lp}2},$$

where $h_{\mathrm{lp}1}$ and $h_{\mathrm{lp}2}$ denote LPFs. More precisely, this represents a narrow full width at half maximum (FWHM) LPF with a larger amplitude $\left(A_1\right)$ subtracted by a broad FWHM LPF with a smaller amplitude $\left(A_2\right)$ . Both HP filters must comply with two rules of thumb³⁷ as follows:

Eq. 3

$$H_{\mathrm{hp}}\left(0\right)=A_1-A_2\text{and}{\left(H_{\mathrm{hp}}\right)}_{\max }⩽A_1.$$

The difference between these filters is that the $h_{\mathrm{lp}1}$ of the former rolls off faster than that of the latter. In this study, Gaussian functions with various standard deviations $\left(\sigma \right)$ are employed for $h_{\mathrm{lp}}$ . Thus, the HP filter shown in Fig. 1 can be formulated as

Eq. 4

$$h_{\mathrm{hp}}\left(r\right)=g_1\left(r\right)-g_2\left(r\right),$$

where $g_1\left(r\right)=[A_1∕{\left(2\pi {\sigma }_1^2\right)}^{1∕2}]\exp (-r^2∕2{\sigma }_1^2)$ and $g_2\left(r\right)=[A_2∕{\left(2\pi {\sigma }_2^2\right)}^{1∕2}]\exp (-r^2∕2{\sigma }_2^2)$ , respectively, and ${\sigma }_1<{\sigma }_2$ .

Fig. 1

Two LPF-combined HP filters (a) $\delta -g_2$ and (b) $g_1-g_2$ , and (c) a wavelet-like filter, where $\delta$ is the delta function.

2.1.1.

LPF-combined HP filter $({\sigma }_1<{\sigma }_2)$

This filter is usually determined with a smaller ${\sigma }_1(⪡1)$ , as shown in Fig. 1a; this can also be determined with a larger ${\sigma }_1$ , as shown in Fig. 1b. If we let ${\sigma }_1$ approach zero, $h_{\mathrm{hp}1}$ narrows further to an impulse, then Eq. 4 can be expressed in the frequency domain as $1-H_{\mathrm{lp}}$ .

2.1.2.

Wavelet-like HP filter

A dilated wavelet-like function expressed as

Eq. 5

$${\psi }_a\left(r\right)=\frac{1}{\sqrt{3a}\sqrt[4]{\pi }}(1-\frac{r^2}{a^2})\exp (-\frac{r^2}{2a^2}),$$

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

2.1.3.

Negative-oriented Laplacian HP filter

Alternatively, a $3\times 3$ negative-orientated Laplacian edge operator

Eq. 6

$$\left[\begin{array}{ccc}0& -1& 0\\ -1& 4& -1\\ 0& -1& 0\end{array}\right]$$

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.

2.3.

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\text{-}\mathrm{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

$$\nabla \cdot D\nabla \Phi (r,\omega )-({\mu }_a-\frac{i\omega }{c})\Phi (r,\omega )=-S(r,\omega ),$$

where $S(r,\omega )$ and $\Phi$ denote the source and the radiance, respectively; and ${\mu }_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\nabla \Phi \cdot \widehat{n}=\alpha \Phi$ (flux in fact) and the FEM are applied. Since only dc data are considered, $\omega$ 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, $\Phi$ , 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

$$\left[\begin{array}{cc}{A}_{bb}-\alpha B_{bb}& A_{bI}\\ A_{Ib}& A_{II}\end{array}\right]\left\{\begin{array}{c}\frac{\partial {\Phi }_b}{\partial D_k}\\ \frac{\partial {\Phi }_I}{\partial D_k}\end{array}\right\}=\left[\begin{array}{cc}-\frac{\partial A_{bb}}{\partial D_k}& -\frac{\partial A_{bI}}{\partial D_k}\\ -\frac{\partial A_{Ib}}{\partial D_k}& -\frac{\partial A_{II}}{\partial D_k}\end{array}\right]\left\{\begin{array}{c}{\Phi }_b\\ {\Phi }_I\end{array}\right\}+\left\{\begin{array}{c}\frac{\partial C_b}{\partial D_k}\\ \frac{\partial C_I}{\partial D_k}\end{array}\right\},$$

and

Eq. 21

$$\left[\begin{array}{cc}{A}_{bb}-\alpha B_{bb}& A_{bI}\\ A_{Ib}& A_{II}\end{array}\right]\left\{\begin{array}{c}\frac{\partial {\Phi }_b}{\partial {\mu }_l}\\ \frac{\partial {\Phi }_I}{\partial {\mu }_l}\end{array}\right\}=\left[\begin{array}{cc}-\frac{\partial A_{bb}}{\partial {\mu }_l}& -\frac{\partial A_{bI}}{\partial {\mu }_l}\\ -\frac{\partial A_{Ib}}{\partial {\mu }_l}& -\frac{\partial A_{II}}{\partial {\mu }_l}\end{array}\right]\left\{\begin{array}{c}{\Phi }_b\\ {\Phi }_I\end{array}\right\}+\left\{\begin{array}{c}\frac{\partial C_b}{\partial {\mu }_l}\\ \frac{\partial C_I}{\partial {\mu }_l}\end{array}\right\},$$

where the superscripts $I$ and $b$ denote interior and boundary nodes, and $D_k$ for $k=1,2,\dots ,K$ and ${\mu }_l$ for $l=1,2,\dots ,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, $\partial {\Phi }_b∕\partial {\mu }_l$ or $\partial {\Phi }_b∕\partial D_k$ , 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

$$({\mathsf{J}}^T\mathsf{J}+\lambda I)\Delta \chi ={\mathsf{J}}^T({\Phi }^o-{\Phi }^c)={\mathsf{J}}^T\Delta \Phi ,$$

where Jacobian matrix $\mathsf{J}$ denotes $\mathsf{J}(\partial {\Phi }_b∕\partial D_k,\partial {\Phi }_b∕\partial {\mu }_l)$ , $\Delta \chi$ means $\Delta \chi (\Delta D_k,\Delta {\mu }_l)$ , and $\lambda$ 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.

3.

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 $\left({\mu }_a\right)$ and reduced scattering $\left({\mu }_s^{\prime }\right)$ values are about 0.0025 and $0.25{\mathrm{mm}}^{-1}$ , respectively, while the maximum absorption and reduced scattering for the inclusion are 0.025 and $2.5{\mathrm{mm}}^{-1}$ , 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, $\delta -g2({\sigma }_2=1.5)$ , $g1-g2({\sigma }_1=0.75,{\sigma }_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 ( ${\mu }_a=0.0025{\mathrm{mm}}^{-1}$ , ${\mu }_s^{\prime }=0.25{\mathrm{mm}}^{-1}$ ) 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\text{-}\mathrm{GHz}$ CPU and $4\mathrm{Gbytes}$ of RAM, respectively. Thirty iteration assignments were employed for each case as the normalized increasing rate, i.e., mean value of ${\mid ({\Phi }_{n+1}-{\Phi }_n)∕{\Phi }_n\mid }^2$ , reaches less than ${10}^{-2}$ , where each iteration takes about $2\min$ . 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.

3.1.

Examples Illustration

3.1.1.

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.

3.1.2.

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).

3.1.3.

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.

3.1.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.

3.2.

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 measures^{42, 43, 44, 45} have been used to evaluate the performance of the NIR imaging algorithms or systems. Song ⁴² 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 ⁴⁴) 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.

3.2.1.

Contrast resolution $\left(R_{\mathrm{cont.}}^{1\mathrm{D},2\mathrm{D}}\right)$

The measure $R_{\mathrm{cont.}}^{1\mathrm{D},2\mathrm{D}}$ 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

$$R_{\mathrm{cont.}}^{1\mathrm{D},2\mathrm{D}}=\frac{{({\overline{\max }}^{\Delta \mathrm{incl.}}∕{\overline{\min }}^{\Delta \overline{\mathrm{incl.}}})}_{\text{reconstruction}}}{{({\overline{\max }}^{\Delta \mathrm{incl.}}∕{\overline{\min }}^{\Delta \overline{\mathrm{incl.}}})}_{\text{original}}},$$

and

Eq. 24

$$R_{\mathrm{cont.}}^{1\mathrm{D},2\mathrm{D}}=2-R_{\mathrm{cont.}}^{1\mathrm{D},2\mathrm{D}},\text{if}1<R_{\mathrm{cont.}}^{1\mathrm{D},2\mathrm{D}}<2,$$

where $\overline{\max }$ and $\overline{\min }$ denote the average of maxima and minima over all the selected regions as the superscripts ( $\Delta \mathrm{incl.}$ and ${\Delta }_{\overline{\mathrm{incl}}}$ ). Also incl. and $\overline{\mathrm{incl}}$ correspondingly represent inclusions and complementary inclusions as well as $\Delta \mathrm{incl.}$ and ${\Delta }_{\overline{\mathrm{incl}}}$ are chosen with several nodes around central area of incl. and incl.

3.2.2.

Separation resolution $\left(R_{\mathrm{sep}}^{1\mathrm{D},2\mathrm{D}}\right)$

The measure $R_{\mathrm{sep}}^{1\mathrm{D},2\mathrm{D}}$ is designed to evaluate the resolution on the separation between inclusions.

Eq. 25

$$R_{\mathrm{sep}}^{1\mathrm{D},2\mathrm{D}}={\left\{[1-\frac{{\left({\mathrm{MSE}}^{\overline{\mathrm{incl.}}}\right)}_{\mathrm{Recon}.2\mathrm{Ori.}}}{{\left({\mathrm{MSE}}^{\overline{\mathrm{incl.}}}\right)}_{\mathrm{Ori.}2\text{Baseline}}}]R_{\mathrm{cont.}}^{1\mathrm{D},2\mathrm{D}}\right\}}^{1∕2}\equiv {({\mathit{Ro}}_{\mathrm{sep}}^{1\mathrm{D},2\mathrm{D}}\times R_{\mathrm{cont.}}^{1\mathrm{D},2\mathrm{D}})}^{1∕2},$$

where MSE is the mean square error over all the selected region as the superscript $\left(\overline{\mathrm{incl}}\right)$ and Baseline is used with $0.025{\mathrm{mm}}^{-1}$ ; Ori.2Baseline is, here, $0.0025\text{to}0.025{\mathrm{mm}}^{-1}$ , and Recon.2Ori. is the reconstructed value in the region $\overline{\mathrm{incl}}$ to $0.0025{\mathrm{mm}}^{-1}$ .

3.2.3.

Size resolution $\left(R_{\text{size}}^{1\mathrm{D},2\mathrm{D}}\right)$

The measure $R_{\text{size}}^{1\mathrm{D},2\mathrm{D}}$ is designed to evaluate the resolution on the size over all inclusions.

Eq. 26

$$R_{\text{size}}^{1\mathrm{D},2\mathrm{D}}={\left\{[1-\frac{{\left({\mathrm{MSE}}^{\mathrm{incl.}}\right)}_{\mathrm{Recon.}2\mathrm{Ori.}}}{{\left({\mathrm{MSE}}^{\mathrm{incl.}}\right)}_{\mathrm{Ori.}2\text{baseline}}}]R_{\mathrm{cont.}}^{1\mathrm{D},2\mathrm{D}}\right\}}^{1∕2}\equiv {({\mathit{Ro}}_{\text{size}}^{1\mathrm{D},2\mathrm{D}}\times R_{\mathrm{cont.}}^{1\mathrm{D},2\mathrm{D}})}^{1∕2},$$

where MSE is over the selected region incl. and baseline is used with $0.0025{\mathrm{mm}}^{-1}$ . 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.

3.2.4.

Location resolution $\left(R_{\mathrm{locat}}^{1\mathrm{D},2\mathrm{D}}\right)$

The measure $R_{\mathrm{locat}}^{1\mathrm{D},2\mathrm{D}}$ is defined to evaluate the resolution on the location over all inclusions.

Eq. 27

$$R_{\mathrm{locat}}^{1\mathrm{D},2\mathrm{D}}={[1-\frac{{\left({\overline{\mathrm{CM}}}^{\overline{\mathrm{incl.}}}\right)}_{\text{Reconstruction}}}{{\left({\overline{\mathrm{CM}}}^{\mathrm{incl.}}\right)}_{\text{Original}}}{R}_{\mathrm{cont.}}^{1\mathrm{D},2\mathrm{D}}]}^{1∕2}\equiv {({\mathit{Ro}}_{\mathrm{locat}}^{1\mathrm{D},2\mathrm{D}}\times R_{\mathrm{cont.}}^{1\mathrm{D},2\mathrm{D}})}^{1∕2},$$

and

Eq. 28

$${\mathit{Ro}}_{\mathrm{locat}}^{1\mathrm{D},2\mathrm{D}}=2-{\mathit{Ro}}_{\mathrm{locat}}^{1\mathrm{D},2\mathrm{D}},\text{if}1<{\mathit{Ro}}_{\mathrm{locat}}^{1\mathrm{D},2\mathrm{D}}<2,$$

where $\overline{\mathrm{CM}}$ is the average of the center of mass over all the selected region as the superscript (incl.).

Note that the resolutions, $R_{\mathrm{sep}}^{1\mathrm{D},2\mathrm{D}}$ , $R_{\text{size}}^{1\mathrm{D},2\mathrm{D}}$ , and $R_{\text{local}}^{1\mathrm{D},2\mathrm{D}}$ include a multiplication operation by $R_{\mathrm{cont.}}^{1\mathrm{D},2\mathrm{D}}$ 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 $R_{\mathrm{sep}}^{1\mathrm{D},2\mathrm{D}}$ , $R_{\text{size}}^{1\mathrm{D},2\mathrm{D}}$ , and $R_{\text{local}}^{1\mathrm{D},2\mathrm{D}}$ are small because the defined measures are quite strict. For overall cases, it is found that location resolution is above 0.95 prior to ${\mathit{Ro}}_{\mathrm{locat}}^{1\mathrm{D},2\mathrm{D}}$ multiplied by $R_{\mathrm{cont.}}^{1\mathrm{D},2\mathrm{D}}$ , 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 0Size	Contrast	Loc. 0Loc.	Sep. 0Sep.	Size 0Size	Contrast	Locx . 0 Locx .	Locy . 0 Locy .
M0	0.46	0.93	0.12	1.00	0.92	0.94	0.58	1.00	0.98
	0.23	0.33		0.35	0.73	0.74		0.76	0.76
M1	0.69	0.66	0.31	0.95	0.99	0.66	0.68	0.96	0.99
	0.46	0.45		0.54	0.82	0.67		0.81	0.82
M2	0.70	0.63	0.30	0.98	0.99	0.65	0.54	1.00	1.00
	0.46	0.44		0.55	0.73	0.59		0.73	0.73
M3	0.84	0.85	0.36	1.00	0.99	0.78	0.87	1.00	0.99
	0.55	0.55		0.60	0.93	0.82		0.93	0.93
M4	0.74	0.71	0.15	1.00	0.97	0.65	0.80	1.00	0.98
	0.33	0.32		0.38	0.88	0.72		0.89	0.88

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 0Size	Contrast	Loc. 0Loc.	Sep. 0Sep.	Size 0Size	Contrast	Locx . 0 Locx .	Locy . 0 Locy .
M0	0.84	0.85	0.10	1.00	0.95	0.87	0.24	1.00	0.98
	0.29	0.29		0.32	0.48	0.46		0.49	0.49
M1	0.94	0.95	0.33	0.98	0.95	0.89	0.95	0.96	0.99
	0.56	0.56		0.57	0.95	0.92		0.95	0.97
M2	0.86	0.91	0.17	0.99	0.97	0.90	0.83	0.99	0.96
	0.38	0.39		0.41	0.90	0.86		0.91	0.89
M3	0.88	0.94	0.17	1.00	0.97	0.95	0.70	1.00	0.99
	0.39	0.41		0.42	0.83	0.82		0.84	0.84
M4	0.87	0.87	0.40	1.00	0.74	$-0.07$	0.50	1.00	0.98
	0.59	0.59		0.63	0.61	$-0.19$		0.71	0.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 0Size	Contrast	Loc. 0Loc.	Sep. 0Sep.	Size 0Size	Contrast	Locx . 0 Locx .	Locy . 0 Locy .
M0	0.85	0.74	0.10	1.00	0.98	0.75	0.18	1.00	1.00
	0.29	0.27		0.32	0.42	0.36		0.42	0.42
M1	0.84	0.95	0.20	0.99	0.92	0.97	0.40	0.98	0.98
	0.41	0.43		0.44	0.61	0.62		0.62	0.62
M2	0.57	0.86	0.12	0.98	0.91	0.94	0.33	0.96	0.99
	0.26	0.32		0.34	0.55	0.56		0.57	0.57
M3	0.91	0.94	0.16	1.00	0.95	0.92	0.38	1.00	0.99
	0.39	0.39		0.41	0.60	0.59		0.62	0.61
M4	0.74	0.80	0.42	1.00	0.82	0.69	0.42	0.99	0.98
	0.56	0.58		0.65	0.59	0.54		0.65	0.65

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 0Size	Contrast	Loc. 0Loc.	Sep. 0Sep.	Size 0Size	Contrast	Locx . 0 Locx .	Locy . 0 Locy .
M0	0.85	0.81	0.09	1.00	0.97	0.82	0.19	0.99	1.00
	0.28	0.27		0.30	0.43	0.40		0.44	0.44
M1	0.79	0.83	0.09	0.98	0.97	0.84	0.40	0.97	0.96
	0.27	0.27		0.30	0.62	0.58		0.62	0.62
M2	0.63	0.81	0.08	0.98	0.94	0.85	0.31	0.96	0.99
	0.22	0.25		0.28	0.54	0.51		0.54	0.55
M3	0.83	0.89	0.11	1.00	0.97	0.87	0.41	1.00	1.00
	0.31	0.32		0.34	0.63	0.60		0.64	0.64
M4	0.83	0.91	0.12	0.99	0.97	0.89	0.44	1.00	0.96
	0.32	0.34		0.35	0.65	0.63		0.66	0.65

3.2.5.

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.

3.2.6.

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.

3.2.7.

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.

3.2.8.

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.

4.

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 ${\sigma }_1$ and ${\sigma }_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.