## 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\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}1000\phantom{\rule{0.3em}{0ex}}\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}\phantom{\rule{0.3em}{0ex}}\mathrm{mm})$
, but low intrinsic contrast
$(\sim {10}^{-1})$
, NIR imaging possesses exceptionally high intrinsic contrast
$(\sim {10}^{1\u20132})$
, but exhibits inferior spatial resolutions
$(\sim {10}^{1}\phantom{\rule{0.3em}{0ex}}\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^{7} 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^{8} applied the Markov-chain technique to enhance optical image resolution. Jiang and Paulsen^{9} 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 Vasus^{25} 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^{26} 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 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,

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

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^{37}as follows:

## Eq. 3

$${H}_{\mathrm{hp}}\left(0\right)={A}_{1}-{A}_{2}\phantom{\rule{1em}{0ex}}\text{and}\phantom{\rule{0.3em}{0ex}}{\left({H}_{\mathrm{hp}}\right)}_{\mathrm{max}}\u2a7d{A}_{1}.$$## 2.1.1.

#### LPF-combined HP filter $({\sigma}_{1}<{\sigma}_{2})$

This filter is usually determined with a smaller ${\sigma}_{1}(\u2aa11)$ , 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}})\mathrm{exp}(-\frac{{r}^{2}}{2{a}^{2}}),$$## 2.1.3.

#### Negative-oriented Laplacian HP filter

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

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

### 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 algorithm^{38, 39, 40} is given as

## Eq. 7

$${\widehat{f}}_{n}={\widehat{f}}_{n-1}\phantom{\rule{0.2em}{0ex}}\mathrm{exp}[(\frac{g}{{\widehat{f}}_{n-1}\otimes h}-1)\ast h]\equiv {\widehat{f}}_{n-1}C,\phantom{\rule{1em}{0ex}}n=1,2,\dots ,N,$$## Eq. 9

$${\widehat{f}}_{n}\sim {\widehat{f}}_{n-1}[1+(\frac{g}{{\widehat{f}}_{n-1}\otimes h}-1)\ast h]={\widehat{f}}_{n-1}+{\widehat{f}}_{n-1}[(\frac{g}{{\widehat{f}}_{n-1}\otimes h}-1)\ast h]\equiv {\widehat{f}}_{n-1}+\Delta {\widehat{f}}_{n-1},$$## Eq. 10

$$\Delta {\widehat{f}}_{n-1}={\widehat{f}}_{n-1}[(\frac{g}{{\widehat{f}}_{n-1}\otimes h}-1)\ast h]={\widehat{f}}_{n-1}[\left(\frac{g-{\widehat{f}}_{n-1}\otimes h}{{\widehat{f}}_{n-1}\otimes h}\right)\ast h],$$## Eq. 11

$$\Delta {\widehat{f}}_{n-1}\sim \frac{{\widehat{f}}_{n-1}}{{\widehat{f}}_{n-1}\otimes h}(\Delta {\widehat{f}}_{n-1}\ast h).$$## Eq. 12

$$\left(\Delta {\widehat{f}}_{n-1}\right)\left(\Delta {\widehat{f}}_{n-1}\right)=\frac{{\widehat{f}}_{n-1}}{{\widehat{f}}_{n-1}\otimes h}(\Delta {\widehat{f}}_{n-1}\ast h)\Delta {\widehat{f}}_{n-1},$$## Eq. 13

$${\widehat{f}}_{n-1}\otimes h={\widehat{f}}_{n-1}\left[\frac{(\Delta {\widehat{f}}_{n-1}\ast h)\Delta {\widehat{f}}_{n-1}}{\left(\Delta {\widehat{f}}_{n-1}\right)\left(\Delta {\widehat{f}}_{n-1}\right)}\right].$$## Eq. 14

$${\widehat{f}}_{n-1}\otimes {h}_{\mathrm{hp}}={\widehat{f}}_{n-1}\otimes {h}_{\mathrm{lp}1}-{\widehat{f}}_{n-1}\otimes {h}_{\mathrm{lp}2}={\widehat{f}}_{n-1}\left[\frac{(\Delta {\widehat{f}}_{n-1}\otimes {h}_{\mathrm{hp}})\Delta {\widehat{f}}_{n-1}}{\left(\Delta {\widehat{f}}_{n-1}\right)\left(\Delta {\widehat{f}}_{n-1}\right)}\right]=\frac{(\Delta {\widehat{f}}_{n-1}\otimes {h}_{\mathrm{hp}})\Delta {\widehat{f}}_{n-1}}{\left(\Delta {\widehat{f}}_{n-1}\right)\left(\Delta {\widehat{f}}_{n-1}\right)\u2215{\widehat{f}}_{n-1}}.$$## Eq. 16

$$\Delta {\widehat{f}}_{n-1}=\frac{(\Delta {\widehat{f}}_{n-1}\otimes {h}_{\mathrm{hp}})\Delta {\widehat{f}}_{n-1}}{\left(\Delta {\widehat{f}}_{n-1}\right)\left(\Delta {\widehat{f}}_{n-1}\right)\u2215{\widehat{f}}_{n-1}}.$$## Eq. 17

$$\Delta {\widehat{f}}_{n-1}=\frac{\u27e8\Delta {\widehat{f}}_{n-1}\mid {h}_{\mathrm{hp}}\mid \Delta {\widehat{f}}_{n-1}\u27e9}{\u27e8\Delta {\widehat{f}}_{n-1}\mid \Delta {\widehat{f}}_{n-1}\u27e9\u2215{\widehat{f}}_{n-1}}\sim \frac{\u27e8\Delta {\widehat{f}}_{n-1}\mid {h}_{\mathrm{hp}}\mid \Delta {\widehat{f}}_{n-1}\u27e9}{w\Vert \Delta {\widehat{f}}_{n-1}\Vert},$$**⋅**∣ and ∣

**⋅**⟩ stand for the state and ⟨ $z\mid x$ and $x\mid y$ ⟩ 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.

## 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 ),$$## 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\},$$## 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\},$$## Eq. 22

$$({\mathsf{J}}^{T}\mathsf{J}+\lambda I)\Delta \chi ={\mathsf{J}}^{T}({\Phi}^{o}-{\Phi}^{c})={\mathsf{J}}^{T}\Delta \Phi ,$$## 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\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$ , respectively, while the maximum absorption and reduced scattering for the inclusion are 0.025 and $2.5\phantom{\rule{0.3em}{0ex}}{\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.

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\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$ , ${\mu}_{s}^{\prime}=0.25\phantom{\rule{0.3em}{0ex}}{\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\phantom{\rule{0.3em}{0ex}}\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})\u2215{\Phi}_{n}\mid}^{2}$ , reaches less than ${10}^{-2}$ , where each iteration takes about $2\phantom{\rule{0.3em}{0ex}}\mathrm{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.

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

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.

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

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
^{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.

## 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{\mathrm{max}}}^{\Delta \mathrm{incl.}}\u2215{\overline{\mathrm{min}}}^{\Delta \overline{\mathrm{incl.}}})}_{\text{reconstruction}}}{{({\overline{\mathrm{max}}}^{\Delta \mathrm{incl.}}\u2215{\overline{\mathrm{min}}}^{\Delta \overline{\mathrm{incl.}}})}_{\text{original}}},$$## Eq. 24

$${R}_{\mathrm{cont.}}^{1\mathrm{D},2\mathrm{D}}=2-{R}_{\mathrm{cont.}}^{1\mathrm{D},2\mathrm{D}},\phantom{\rule{1em}{0ex}}\text{if}\phantom{\rule{0.3em}{0ex}}1<{R}_{\mathrm{cont.}}^{1\mathrm{D},2\mathrm{D}}<2,$$## 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\u22152}\equiv {({\mathit{Ro}}_{\mathrm{sep}}^{1\mathrm{D},2\mathrm{D}}\times {R}_{\mathrm{cont.}}^{1\mathrm{D},2\mathrm{D}})}^{1\u22152},$$## 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\u22152}\equiv {({\mathit{Ro}}_{\text{size}}^{1\mathrm{D},2\mathrm{D}}\times {R}_{\mathrm{cont.}}^{1\mathrm{D},2\mathrm{D}})}^{1\u22152},$$## 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\u22152}\equiv {({\mathit{Ro}}_{\mathrm{locat}}^{1\mathrm{D},2\mathrm{D}}\times {R}_{\mathrm{cont.}}^{1\mathrm{D},2\mathrm{D}})}^{1\u22152},$$## Eq. 28

$${\mathit{Ro}}_{\mathrm{locat}}^{1\mathrm{D},2\mathrm{D}}=2-{\mathit{Ro}}_{\mathrm{locat}}^{1\mathrm{D},2\mathrm{D}},\phantom{\rule{1em}{0ex}}\text{if}\phantom{\rule{0.3em}{0ex}}1<{\mathit{Ro}}_{\mathrm{locat}}^{1\mathrm{D},2\mathrm{D}}<2,$$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.

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

## 3.3.

### 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

$${R}_{\mathrm{cont.}}^{1\mathrm{D},2\mathrm{D}}=\frac{{({\text{mean}}^{\mathrm{incl.}}\u2215{\text{mean}}^{\overline{\mathrm{incl.}}})}_{\text{Reconstruction}}}{{({\text{mean}}^{\mathrm{incl.}}\u2215{\text{mean}}^{\overline{\mathrm{incl.}}})}_{\text{Original}}},$$## Eq. 30

$${R}_{\mathrm{cont.}}^{1\mathrm{D},2\mathrm{D}}\{\begin{array}{ll}>c& \text{normal}\phantom{\rule{0.3em}{0ex}}\text{situation}\\ =c& \text{no}\phantom{\rule{0.3em}{0ex}}\text{contrast}\\ <c& \text{abnormal}\phantom{\rule{0.3em}{0ex}}\text{situation}\end{array}\phantom{\}},$$## Eq. 31

$${\mathit{Ro}}_{\mathrm{sep};\text{size}}^{1\mathrm{D},2\mathrm{D}}=\frac{{\left({\mathrm{MSE}}^{\overline{\mathrm{incl.}};\mathrm{incl.}}\right)}_{\mathrm{Recon.}2\mathrm{B}\left(\mathrm{b}\right)\mathrm{aseline}}}{{\left({\mathrm{MSE}}^{\overline{\mathrm{incl.}};\mathrm{incl.}}\right)}_{\mathrm{Ori.}2\mathrm{B}\left(\mathrm{b}\right)\mathrm{aseline}}}.$$## 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.

## Acknowledgments

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.

## References

*in vivo*breast tissue: implementation of a Laplacian-type regularization to incorporate magnetic resonance structure,” J. Biomed. Opt., 10 051504 (2005). https://doi.org/10.1117/1.2098627 1083-3668 Google Scholar

*a-priori*knowledge,” 74 –77 (2002). Google Scholar

*in vivo*breast tissue: implementation of a Laplacian-type regularization to incorporate magnetic resonance structure,” J. Biomed. Opt., 10 051504 (2005). https://doi.org/10.1117/1.2098627 1083-3668 Google Scholar