## 1.

## Introduction

The investigation of the human body by diffusely propagated light is a highly appealing quest due to the intrinsic noninvasiveness of light at low power and the wealth of information carried by absorption (related to chemical composition) and scattering (linked to microstructure). Thus, a range of different applications have been devised and initially tested in clinics, such as optical mammography for the detection of breast lesions,^{1}2.^{–}^{3} monitoring of neoadjuvant chemotherapy,^{4} assessment of breast density as a cancer risk factor,^{5} functional imaging of the brain,^{6}7.^{–}^{8} monitoring of blood perfusion and oxygenation in the injured brain,^{9} muscle oximetry,^{10} studies of epilepsy,^{11}^{,}^{12} and investigation of bone and joint pathologies,^{13} just to cite same major examples.

Recently, the compelling need for standardization and quality assessment of diffuse optics instruments has been pointed out as a key requirement for the translation of the optical tools to effective clinical use.^{14} Indeed, the large variety of techniques and instruments, the diverse nature of clinical problems as well as the need of advanced theoretical modeling to extract clinically relevant parameters from the measurements demand some clearly identified and shared procedures to test and validate the optical systems. Common standardization tools help in comparing results obtained in different clinical studies, permit the grading of different instruments or subsequent upgrades of the same system, address the research—particularly at the industrial level—toward those parameters that best match the clinical needs. The work toward standardization requires, on one hand, the definition of shared protocols, identifying the key physical quantities to be measured and the procedures for their assessment^{14}15.16.17.^{–}^{18} and, on the other hand, the description and production of phantoms mimicking some paradigmatic *in vivo* problems. For those problems that concern the protocols, several successful multilaboratory attempts have been initiated at the European Union (EU) level—fostered by EU-funded projects—such as the MEDPHOT protocol for the assessment of diffuse optics instruments on homogeneous media^{15} or, more recently, for characterization of time-domain brain imagers, the Basic Instrumental Performance protocol^{16} as well as the nEUROPt protocol,^{17} the latter being related to localized absorption changes.

Several studies have been focused on the design of tissue phantoms over the last decades, with the contribution^{19}20.21.22.23.24.25.^{–}^{26} of many groups (see Ref. 27 for a recent review). For homogeneous liquid phantoms, a high accuracy in characterization of Intralipid and India ink has been achieved.^{28}^{,}^{29} A recent multicenter study of nine different laboratories^{30} resulted in common mean values for the intrinsic scattering of Intralipid and the intrinsic absorption of ink with an associated uncertainty of 1% and 2%, respectively. Solid homogenenous phantoms can also be accurately characterized and are commercially available with a stated uncertainty of about 11% for the absorption coefficient, ${\mu}_{a}$, and about 7% for the reduced scattering coefficient, ${\mu}_{s}^{\prime}$ (Ref. 31). Conversely, the production of inhomogeneous phantoms—mimicking, for instance, a focal activation in the brain cortex or a tumor mass within the breast—is more challenging, due to the added complexity of fabricating the structure by solid–solid, solid–liquid, or liquid–liquid phases. The possible artifacts due to the solid–liquid interface, or to the material layer separating two liquid phases, raise concern. Also, the large number of possible combinations, involving not only the host medium optical properties but also the location of source and detector, size, position, and properties of the local inhomogeneities, together with the source–detector arrangement, add to the overall complexity of the problem. Finally, the construction of dynamic phantoms, with the chance to easily and quickly vary phantom parameters, is a further need.

Just to quote a few examples to understand the complexity of inhomogeneous phantoms, the use of agar or agarose was proposed to create solid inclusions starting from Intralipid and ink dilutions,^{21}^{,}^{32} yet with some criticalities, like scattering reduction produced by agar requiring proper compensation, diffusion of chromophores into and from the inhomogeneity, complex fabrication procedure, and fragility in handling and conservation. The use of transparent tubes, even with thin walls, can produce artifacts because of light-guiding effects. Solid inhomogeneities made, for example, of resin are definitely practical for use in liquid environments, yet refractive-index mismatch and potentially not-perfect matching of scattering properties could be an issue.^{33}

In this paper, we will show that the effect of a localized variation of the absorption coefficient within a finite volume can be adequately reproduced using a small totally absorbing object with a proper volume (equivalence relation). This concept can be useful in two ways. First, it can be exploited to design realistic inhomogeneous phantoms using a set of black objects with different volumes. Further, it permits one to grade physiological or pathological changes on a reproducible scale of black-equivalent volumes, thus facilitating the performance assessment of clinical instruments.

Although totally absorbing objects (e.g., black PVC cylinders) have already been used in diffuse spectroscopy experiments,^{34}35.36.37.^{–}^{38} their far-ranging equivalence to realistic absorption inhomogeneities has never been thoroughly investigated. Also, the effect of a totally absorbing volume has never been quantified. The present study encompasses both continuous wave (CW) and time-resolved measurements, with an obvious extension to the frequency-domain case.

In the following, we will first present the equivalence relation from basic theoretical inferences (Sec. 2) and the methods used (Sec. 3). Then we will investigate its applicability with Monte Carlo (MC) simulations spanning a wide range of possible scenarios (Sec. 4). We will provide a set of look-up plots and interpolating functions to derive the equivalent black volume (EBV) yielding the same effect as a realistic absorption object (Sec. 5). Finally, an example of the use of EBV for grading different optical perturbations will be shown (Sec. 6).

This work will be followed by a companion paper where the EBV approach is actually implemented for the construction of simple and highly versatile inhomogeneous phantoms. All specifications on the phantom construction will be provided, together with an extensive numerical and experimental validation. It is worth mentioning that the EBV concept and the related phantom forms an essential basis of the nEUROPt protocol^{17} for the Performance Assessment of Time-Domain Optical Brain Imagers. This protocol has been agreed to by a consortium of 17 partners all over Europe and has been applied to the clinical brain imagers and additional laboratory instruments developed by four groups of this consortium.

## 2.

## Equivalence Relation Between Black and Realistic Inclusions: Formulation within the Diffusion Theory

Biological tissue is a heterogeneous medium and several imaging applications to tissue, such as the case of a focal activation in the brain cortex or the case of a tumor mass within the breast, are related to the detection of an absorption variation confined within a certain volume in such a medium. In this section, we refer to a diffusive medium in which an absorption variation takes place inside a localized volume of the medium (see Fig. 1 for a schematic).

When a small single absorbing object is placed inside a diffusive medium (homogeneous or heterogeneous), it is possible to formulate an equivalence relation so that a small and completely black object can give the same result for the perturbation as that for a larger object with only moderate absorption change. By small object we mean an absorbing inclusion of volume ${V}_{i}$ placed in a position ${\mathbf{r}}_{\mathbf{0}}$ inside the diffusive medium so that the unperturbed Green function is nearly constant over ${V}_{i}$. The time-dependent fluence rate can be expressed as^{39}

## (1)

$$\mathrm{\Phi}(\mathbf{r},t)={\mathrm{\Phi}}^{0}(\mathbf{r},t)-{\int}_{{V}_{i}}\mathrm{d}{\mathbf{r}}^{\prime}\mathrm{\Delta}{\mu}_{a}({\mathbf{r}}^{\prime})\phantom{\rule{0ex}{0ex}}{\int}_{-\infty}^{+\infty}\mathrm{d}{t}^{\prime}\mathrm{\Phi}({\mathbf{r}}^{\prime},{t}^{\prime})G(\mathbf{r},{\mathbf{r}}^{\prime},t,{t}^{\prime}),$$## (2)

$$\mathrm{\Phi}(\mathbf{r},t)={\mathrm{\Phi}}^{0}(\mathbf{r},t)-{\int}_{{V}_{i}}\mathrm{d}{\mathbf{r}}^{\prime}\mathrm{\Delta}{\mu}_{a}({\mathbf{r}}^{\prime})\phantom{\rule{0ex}{0ex}}{\int}_{-\infty}^{+\infty}\mathrm{d}{t}^{\prime}{\mathrm{\Phi}}^{0}({\mathbf{r}}^{\prime},{t}^{\prime})G(\mathbf{r},{\mathbf{r}}^{\prime},t,{t}^{\prime}).$$## (3)

$$\mathrm{\Phi}(\mathbf{r},t)={\mathrm{\Phi}}^{0}(\mathbf{r},t)[1-P(\mathbf{r},{\mathbf{r}}_{\mathbf{0}},{\mu}_{s}^{\prime},t)Q(\mathrm{\Delta}{\mu}_{a},{V}_{i})],$$## (4)

$$P(\mathbf{r},{\mathbf{r}}_{\mathbf{0}},{\mu}_{s}^{\prime},t)={\mathrm{\Phi}}^{0}{(\mathbf{r},t)}^{-1}{\int}_{-\infty}^{+\infty}\mathrm{d}{t}^{\prime}{\mathrm{\Phi}}^{0}({\mathbf{r}}_{\mathbf{0}},{t}^{\prime})G(\mathbf{r},{\mathbf{r}}_{\mathbf{0}},t,{t}^{\prime})$$## (5)

$$Q(\mathrm{\Delta}{\mu}_{a},{V}_{i})={\int}_{{V}_{i}}\mathrm{d}{\mathbf{r}}^{\prime}\mathrm{\Delta}{\mu}_{a}({\mathbf{r}}^{\prime}).$$## (6)

$$C(t)=[\mathrm{\Phi}(t)-{\mathrm{\Phi}}^{0}(t)]/{\mathrm{\Phi}}^{0}(t)=P(\mathbf{r},{\mathbf{r}}_{\mathbf{0}},{\mu}_{s}^{\prime},t)Q(\mathrm{\Delta}{\mu}_{a},{V}_{i}),$$## (7)

$${C}_{\mathrm{CW}}=({\mathrm{\Phi}}_{\mathrm{CW}}-{\mathrm{\Phi}}_{\mathrm{CW}}^{0})/{\mathrm{\Phi}}_{\mathrm{CW}}^{0}={P}_{\mathrm{CW}}(\mathbf{r},{\mathbf{r}}_{\mathbf{0}},{\mu}_{s}^{\prime},{\mu}_{a})Q(\mathrm{\Delta}{\mu}_{a},{V}_{i}),$$Paasschens and Hooft,^{35} making use of an iterative procedure and by expressing $\mathrm{\Phi}$ for the frequency domain and for an infinite medium as a Neumann series, extended the validity of Eq. (3) beyond the Born approximation. In their proof, Eq. (3) holds exactly within the diffusion approximation and for small absorbing objects. Yet, it can accommodate even large values for the absorption of the inhomogeneities.

Within the range of validity of this approximation, we find that the relative perturbation depends only on the position of the object and on its strength, while the shape of the object and the precise distribution of the absorption inside the object are not important. Accordingly, a small and completely black object can give the same relative perturbation as a larger object, with only moderate absorption. The equivalence relation can also be summarized noting that, according to Eqs. (6) and (7), the shapes of the time-dependent relative perturbation, $C(t)$, and of the CW relative perturbation, ${C}_{\mathrm{CW}}$, are independent of the intensity of the perturbation $Q$.

When the diffusion equation cannot be applied because the inclusion is too close to the boundary, to the source, or to the detector, the validity of the equivalence relation needs to be checked. With this purpose, by using numerical solutions of the radiative transfer equation reconstructed with an MC code, in Sec. 4, we have investigated the validity of the equivalence relation beyond the diffusive regime. It is also worth stressing that a validation of this relation is also necessary because in those regions, such as black inclusions, where absorption is much higher than scattering, the diffusion approximation is not expected to hold. Although Paasschens and Hooft have derived an expression for the strength of the perturbation due to a small black sphere, their theory shows discrepancies compared with the results of MC simulations. For this reason, a validation of this relation is necessary in order to extend its use to general conditions of photon migration: Sec. 4 will be devoted to this purpose for the time domain and for the CW domain. The results presented in this section cover a wide range of possible geometrical scenarios.

## 3.

## Methods

## 3.1.

### Monte Carlo Code

A perturbation MC code was used to simulate photon migration through absorbing objects placed inside turbid media,^{39}40.^{–}^{41} in order to have reference data for our investigations. The MC method provides a physical simulation of photon migration inside the medium. Thus solutions of the radiative transfer equation can be obtained without any need for simplifying assumptions. Therefore the solutions obtained can be used as a standard reference for validation of approximate theories. In particular, in Sec. 4, MC data have been used to investigate the validity of the equivalence relation beyond the diffusive regime of light propagation and in Sec. 5 MC data have been used for measuring the EBV of an absorbing object. The perturbation MC used to generate the data shown in Secs. 4 and 5, strictly speaking, is not an elementary MC code since it makes use of scaling relationships so that the temporal response, when scattering or absorbing inhomogeneities are introduced, can be evaluated in a short time from the results of only one simulation carried out for the homogeneous medium.^{39}^{,}^{40} We stress that an elementary MC is not suitable for calculating, with enough good statistics, the effects of absorbing inhomogeneities like those considered in our investigation.

## 3.2.

### Determination of the Equivalent Black Volume

The procedure we adopted to derive the volume of the totally absorbing sphere that is equivalent to a given realistic inclusion is based on the comparison of the temporal contrast produced by any target object (realistic or totally absorbing) to the one yielded by an arbitrarily chosen reference object (e.g., a totally absorbing sphere of $1\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}$ volume). The procedure is explained in Fig. 2. First, the temporal contrast $C(t)$ is calculated, both for the target and the reference objects, for a large set of geometries (locations of the inclusion), and for a given value of the background reduced scattering coefficient [$C(t)$ does not depend on the background absorption]. The figure refers to a semi-infinite medium with ${\mu}_{s}^{\prime}=1\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$; the source is at (0,0,0) and the receiver at (30,0,0) mm. The target object is a volume of $1000\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}$ with $\mathrm{\Delta}{\mu}_{a}=0.02\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$. The reference object is a totally absorbing sphere of $1\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}$. The two plots in Fig. 2 display the contrast function $C(t)$ for the target (upper graph) and the reference (lower graph) object in a given geometry. It can be observed that the shape of the $C(t)$ functions is much the same in the two cases. Then, the area $A$ under the $C(t)$ distribution, $A={\int}_{0}^{+\infty}C(t)\mathrm{d}t$, is calculated for both objects and for all the source–detector positions with respect to the inclusion, obtaining a set of $A$ and ${A}_{0}$ values for the target and reference objects, respectively. Figure 2 shows an example of the set of $A$ (top table) and ${A}_{0}$ (middle table) values calculated for different choices of the lateral position ${x}_{i}$ (heading row of the tables) and of the depth ${z}_{i}$ of the object (left column of the tables). It is worth pointing out that the value of $A$ is not the contrast produced by the perturbation in a CW case. Although there is an integral over the photon traveling time $t$, this is applied after dividing $\mathrm{\Delta}\varphi (t)$ for $\varphi (t)$. Finally, the measure $m$ of the target object is derived as the average of the ratios $A/{A}_{0}$ (bottom table) of the two datasets. Thus, a realistic inclusion and a totally absorbing inclusion are equivalent if they yield the same measure $m$ with respect to the reference totally absorbing sphere. Other procedures to find the set of equivalent objects can be designed, yet all results presented in the following sections were obtained using the method described above.

## 4.

## Validation of the Equivalence Relation

## 4.1.

### Time-Resolved Case

The validity of the equivalence relation was tested on MC simulations by comparing the effect on the photon time-distribution produced by a realistic optical inhomogeneity or by its equivalent totally absorbing object determined with the method described in Sec. 3.2. Figure 3 (reflectance) and Fig. 4 (transmittance) show the plots of the relative temporal contrast $|C(t)|$ for a wide range of possible scenarios. In particular, four different realistic inhomogeneities with a volume of $V=1000\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}$, and four different values of $\mathrm{\Delta}{\mu}_{a}$ (rows), are set at different depths (columns) (${z}_{i}/\mathrm{mm}=10,15,20$ for reflectance and ${z}_{i}/\mathrm{mm}=30,25,20$ for transmittance) and compared to the corresponding totally black inclusions, yielding a similar perturbation. The refractive index of the medium and of the external was 1.33. The plots display the contrast $C(t)$ for six different arrangements of source, inclusion, and detector. In Fig. 3, we have reflectance measurements at $\rho =5\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$ (R5) and $\rho =30\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$ (R30 and R30S): the inclusion is set midway for R5 and R30 (${x}_{i}/\mathrm{mm}=2.5$, 15; ${y}_{i}=0$) between source [placed at (0,0,0)] and detector [placed at (5,0,0) mm for R5 and at (30,0,0) mm for R30] and with an $x$-offset with respect to the midpoint of 15 mm for R30S (that has ${x}_{i}=30\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$, ${y}_{i}=0$). In Fig. 4, we have transmittance measurements through a 4-cm thick slab (for three configurations of the inclusion), with the detector at $x=y=0$ and $z=40\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$: the inclusion is set on axis for (${x}_{i}={y}_{i}=0$) and off-axis with an $x$-offset of 10 and 20 mm, respectively, i.e., ${x}_{i}=10\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$, ${y}_{i}=0$ and ${x}_{i}=20\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$, ${y}_{i}=0$. Generally, there is a good correspondence between the shapes of the temporal profiles for the entangled couples of objects, which change consistently upon varying the measurement geometry of the depth of the inclusion, matching, in most cases, also the absolute value of the contrast. The main discrepancies arise for low absorbing perturbations (first row) or a shallower depth (first column). The discrepancy is easily explained for geometrical reasons. When the absorption contrast is low, the totally absorbing inclusion is small (e.g., $5\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}$); thus fewer photon trajectories are intercepted by the black object as compared to the much larger ($1000\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}$) realistic perturbation. Moreover, this difference is more marked for shallower objects. Transmittance measurements with the inhomogeneity near the axis are more affected by this discrepancy since, for objects nearer to the source or to the detector, the strong reduction in size of the perturbation has a marked effect on the interception of photon trajectories.

## 4.2.

### Continuous Wave Case

In Sec. 2, we pointed out that $C(t)$ is rigorously independent of the background absorption. Conversely, in the CW domain, both the ${P}_{\mathrm{CW}}$ function and the CW contrast depend on ${\mu}_{a}$, and the equivalence relation must also be validated against different choices of ${\mu}_{a}$. Figures 5 and 6 show the plots of the relative contrast ${C}_{\mathrm{CW}}$ as a function of the absorption coefficient of the background ${\mu}_{a}$ for the same four couples of realistic/black inhomogeneities (rows) set at several depths $z$ (columns) and the six geometries considered in Figs. 3 and 4. Overall, there is a generally good agreement between the black and the realistic perturbations, with largest discrepancies for the less-absorbing and shallower objects. Discrepancies are larger for larger values of ${\mu}_{a}$.

All simulations shown above were obtained for a background scattering coefficient ${\mu}_{s}^{\prime}=1\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$. For a different choice of ${\mu}_{s}^{\prime}=1\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$, a similar pattern is obtained, yet with different values of the totally absorbing inclusion.

Both for the time-resolved and the CW case, the discrepancy between the couples of objects is reduced for a lower volume (e.g., $500\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}$) of the realistic inhomogeneity (data not shown). Generally speaking, for a given ${\mu}_{s}^{\prime}$, the equivalence relation identifies—within the set of all possible (confined) absorption inhomogeneities—subclasses of equivalent objects producing the same effect. Each subclass contains only one totally absorbing spherical inclusion but an infinite variety of absorption inhomogeneities with different combinations of $\mathrm{\Delta}{\mu}_{a}$, volume, and shape. Thus, the equivalence relation can also be used to find matches between absorption inhomogeneities of different size.

## 5.

## Equivalent Black Volume

Since any realistic absorption perturbation—as far as the equivalence relation holds true—can be described by its equivalent totally absorbing inclusion, it is possible to use the volume of this inclusion as a measure of the perturbation. It is important to note that, within the same assumption, the equivalent volume is approximately the same, independent of the source–detector location, of the position of the perturbation, of ${\mu}_{a}$, while it changes for a different ${\mu}_{s}^{\prime}$. We first define a unitary EBV as

## (8)

$$\mathrm{\Delta}{\mu}_{a}({V}_{\text{black}})=\alpha {({V}_{\text{black}})}^{\beta}\phantom{\rule[-0.0ex]{14.8em}{0.0ex}}({V}_{\text{black}})<30\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}\phantom{\rule{0ex}{0ex}}\mathrm{\Delta}{\mu}_{a}({V}_{\text{black}})=\sum _{n=0}^{6}{C}_{n}{({V}_{\text{black}})}^{n}\text{\hspace{0.17em}}30\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}<({V}_{\text{black}})<600\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}.\phantom{\rule{0ex}{0ex}}$$## Table 1

Coefficients of the interpolating function of Eq. (8) retrieved using MC data for μs′=0.5, 1, and 2 mm−1. The correlation coefficient R is also shown.

Coefficients | μs′=0.5 mm−1 | μs′=1 mm−1 | μs′=2 mm−1 | R |
---|---|---|---|---|

α | 1.07E-3 | 9.34E-4 | 7.03E-4 | 0.999 |

β | 6.59E-1 | 6.05E-1 | 5.56E-1 | 0.999 |

C0 | 4.15E-3 | 2.91E-3 | 1.93E-3 | 1 |

C1 | 2.14E-4 | 1.58E-4 | 9.58E-5 | 1 |

C2 | −5.92E-7 | −4.01E-7 | −1.80E-7 | 1 |

C3 | 2.66E-9 | 1.85E-9 | 7.62E-10 | 1 |

C4 | −5.28E-12 | −3.36E-12 | −9.28E-13 | 1 |

C5 | 5.67E-15 | 3.26E-15 | 4.58E-16 | 1 |

C6 | −1.74E-18 | −6.14E-19 | 4.20E-19 | 1 |

The concept of EBV makes possible the construction of realistic inhomogeneous phantoms for imaging purposes using a set of black objects with different volumes, By changing the volume, we change the intensity of the perturbation simulated in the optical phantom. Making use of this approach, the calibration of the optical properties of the inclusion is not necessary and the characterization of the absorption properties is based only on the equivalence relation. We stress that with this concept, the kind of perturbation that can be simulated is restricted to the condition $\mathrm{\Delta}{\mu}_{a}>0$. When $\mathrm{\Delta}{\mu}_{a}<0$, a totally absorbing object cannot be used to reproduce the same perturbation of a realistic inclusion. However, for small values of $\mathrm{\Delta}{\mu}_{a}$, the case $\mathrm{\Delta}{\mu}_{a}<0$ is complementary to the case $\mathrm{\Delta}{\mu}_{a}>0$. Within the Born approximation, it is trivial to show that for an absorbing inclusion of fixed volume $C(\mathrm{\Delta}{\mu}_{a},t)=-C(-\mathrm{\Delta}{\mu}_{a},t)$. For large values of $\mathrm{\Delta}{\mu}_{a}$, this property does not exactly hold true and some corrections need to be introduced. Finally, we also point out that within the range of validity of the equivalence relation, the results presented in this section are applicable both to the time domain and to the CW domain.

## 6.

## Examples of EBV Assessment

To better envision the potential use of the EBV measure, we report in Table 2 some examples of quantification of different perturbations for the particular applications of functional imaging of the brain and of optical mammography. Most of the results used referred to time-domain experiments. In particular, the table compares the measure of different optical perturbations, namely (a) a totally absorbing sphere of $1\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}$ volume, (b) an increase in absorption of $0.01\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}$ over a volume of $1000\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}$, (c) an increase of oxy-(${\mathrm{HbO}}_{2}$) or deoxy-(Hb) hemoglobin of 100 nmol assessed at 800 nm, (d) the effect of a typical motor cortex activation assessed *in vivo*,^{42}^{,}^{43} (f) the effect of a breast malignant lesion assessed *in vivo*.^{4}^{,}^{5} All these perturbation effects can be compared in a quantitative way, since all of them can be graded using the EBV. In general, it is not necessary to specify the measurement geometry or the background optical properties, since the equivalence holds true within any combination of these parameters, yet within the limitations discussed above. However, in order to calculate the EBV corresponding to our target perturbation, the information on geometry and optical properties is required. This information is necessary for implementing the forward calculations used to generate the perturbations of the total absorbing objects to be matched to the measured perturbation. Therefore, for cases b, d, and f, the realistic volume of $1000\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}$ was assumed for calculating the EBV according to the previous section and a ${\mu}_{s}^{\prime}=1\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$ was assumed for determining the EBV corresponding to cases b and c. Finally, for the last two cases, the calculation was carried out assuming a 15-mm depth (case d) and the lesion at midway in a thickness of 40 mm (case f). Truly, the assessment of the EBV for the *in vivo* data requires a rough estimate of both the inclusion depth and the background reduced scattering coefficient. In particular, the depth dependence is quite steep, leading to a substantial dependence of the estimated EBV on the assumed depth. Yet, once the EBV measure is derived, it can be compared to other *in vivo* results independent of the measurement system or geometry.

## Table 2

Examples of absorption perturbations expressed as equivalent black volume.

Perturbation | EBV (mm3) | Wavelength (nm) |
---|---|---|

mm3 totally absorbing sphere | 1 | — |

Δμa=0.01 mm−1 | 50 | — |

nmol Hb | 120 | 800 |

Motor task | ≈10 | 830 |

Malignant breast lesion | ≈100 | 635 |

The evaluation of the EBV units of a certain perturbation requires a forward solver for identifying the EBV. For producing the results presented in the previous sections, we have largely used an MC forward solver. The main limitation of the MC is the computation time since the evaluation of a certain perturbation should be preferably done in real time. For this reason, as an alternative to the MC, we have also developed an analytical forward solver based on a high-order perturbation theory^{44}45.46.^{–}^{47} and on the equivalence relation. It is beyond the scope of this work, which is dedicated to the basic concepts, to provide details on this solver. What we want to stress is that the evaluation of the EBV units using the above-mentioned solver can be done in real time with any standard PC and with negligible differences compared to the results obtained with an MC solver.

The *in vivo* result presented above is just an example taken from the dataset collected in our laboratory and has no claim to represent an average or typical response. Rather, it is a realistic example presented to show the applicability of the approach and a rough idea on the quantification of an *in vivo* exercise. It would be worth applying this method in a systematic way to a collection of *in vivo* data to provide a quantitative description of *in vivo* activation to be used for tailoring phantom construction, or even to compare the results of different studies. These studies could also be greatly beneficial for the definition of standards, with the construction of an easily reproducible scale of “universal absorption inhomogeneities.” Similarly, the same exercise could be engaged for the quantification of breast lesions in optical mammography. Clearly, all this goes beyond the scope of the present paper.

## 7.

## Conclusions

In summary, we demonstrated that the equivalence relation, formulated in Sec. 2 within the diffusion theory, between a realistic and a totally absorbing inclusion has a rather broad range of validity (see Sec. 4). Its only limitation arises whenever there is a large discrepancy in size (as compared to the distance from the source–detector pair) between the realistic and the black inclusions. Practically this happens when the $\mathrm{\Delta}{\mu}_{a}$ (between realistic inclusion and background) and the EBV are small and the inclusion is located quite shallowly and close either to the source or the detector. From Figs. 3, 4, 5, and 6 we notice that the discrepancies are larger in the transmittance when the inclusion is in front of the beam and in shallow positions. Except for these cases, discrepancies remain of the order of few percents. The complete overview provided with the presented figures can be used to select the optimal range where the equivalence relation holds with negligible error, or alternatively to account for the error due to its use. The range of absorption values where the equivalence relation shows larger approximations is also typical of some applications like breast tumors and brain functional imaging. Therefore, the use of the proposed approach with imaging instrumentation has to be carefully evaluated when the inclusion is in a shallow position. Finally, the equivalence relation can be used indifferently in the time domain and in the CW domain and, therefore, also in the frequency domain.

Based on the equivalence relation, we propose the construction of well-reproducible inhomogeneous phantoms where the defects are represented by totally absorbing objects, like black PVC cylinders, that produce the same effect of a realistic absorption inhomogeneity. Therefore, inhomogeneous phantoms for imaging purposes can be designed using a set of black objects with different volumes: the intensity of the perturbation simulated in the optical phantom can be changed by varying the volume of the black object used. Making use of this approach, the calibration of the optical properties of the defect is not necessary and the characterization of the absorption properties of the inclusion is based only on the equivalence relation. This approach is addressed to absorbing objects embedded in diffusive media and can be applied to design either liquid or solid phantoms.

Based on the equivalence relation, we have also introduced the EBV as a measure of a realistic perturbation defined as the volume of the equivalent totally absorbing inclusion that gives the same perturbation. In this way, physiological or pathological changes can be graded on a reproducible scale, thus facilitating the performance assessment of clinical instruments. Some interpolating functions derived from MC data have been provided in order to facilitate this task. Therefore, by using the functions of Eq. (8), it is possible to calculate the absorption variation $\mathrm{\Delta}{\mu}_{a}$ of a realistic inclusion of volume $1000\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}$ that gives the same perturbation of a black volume. Finally, as an example of the measure in EBV units, we have expressed some perturbations for the particular applications of functional imaging of the brain and of optical mammography in this proposed units. These two realistic examples of perturbations (Table 2) can be represented in the range of about 10 to 100 EBV units and thus can be potentially reproduced with optical phantoms based on totally absorbing objects.

## Acknowledgments

The research leading to these results has received funding from the European Community’s Seventh Framework Programme [FP7/2007-2013] under grant agreement No. FP7-HEALTH-F5-2008-201076. The research leading to these results has also received funding from LASERLAB-EUROPE (grant agreement No. 284464, European Community’s Seventh Framework Programme). This research is also supported by NIH Grant: R01-CA154774 and NSF Award: IIS 1065154.

## References

A. Torricelliet al., “Use of a nonlinear perturbation approach for *in vivo* breast lesion characterization by multiwavelength time-resolved optical mammography,” Opt. Express 11(8), 853–867 (2003).OPEXFF1094-4087http://dx.doi.org/10.1364/OE.11.000853Google Scholar

P. Taroniet al., “Seven-wavelength time-resolved optical mammography extending beyond 1000 nm for breast collagen quantification,” Opt. Express 17(18), 15932–15946 (2009).OPEXFF1094-4087http://dx.doi.org/10.1364/OE.17.015932Google Scholar

D. R. Leffet al., “Diffuse optical imaging of the healthy and diseased breast: a systematic review,” Breast Cancer Res. Treat. 108(1), 9–22 (2008).BCTRD60167-6806http://dx.doi.org/10.1007/s10549-007-9582-zGoogle Scholar

R. Choeet al., “Diffuse optical tomography of breast cancer during neoadjuvant chemotherapy: a case study with comparison to MRI,” Breast Cancer Res. Treat. 32(4), 1128–1139 (2005).BCTRD60167-6806http://dx.doi.org/10.1118/1.1869612Google Scholar

P. Taroniet al., “Effects of tissue heterogeneity on the optical estimate of breast density,” Biomed. Opt. Express 3(10), 2411–2418 (2012).BOEICL2156-7085http://dx.doi.org/10.1364/BOE.3.002411Google Scholar

E. Molteniet al., “Load-dependent brain activation assessed by time-domain functional near-infrared spectroscopy during a working memory task with graded levels of difficulty,” J. Biomed. Opt. 17(5), 056005 (2012).JBOPFO1083-3668http://dx.doi.org/10.1117/1.JBO.17.5.056005Google Scholar

M. Buttiet al., “Effect of prolonged stimulation on cerebral hemodynamic: a time-resolved fNIRS study,” Med. Phys. 36(9), 4103–4114 (2009).MPHYA60094-2405http://dx.doi.org/10.1118/1.3190557Google Scholar

T. Durduranet al., “Diffuse optics for tissue monitoring and tomography,” Rep. Prog. Phys. 73(7), 076701 (2010).RPPHAG0034-4885http://dx.doi.org/10.1088/0034-4885/73/7/076701Google Scholar

M. N. Kimet al., “Noninvasive measurement of cerebral blood flow and blood oxygenation using near-infrared and diffuse correlation spectroscopies in critically brain-injured adults,” Neurocrit. Care 12(2), 173–180 (2010).NCEACB1541-6933http://dx.doi.org/10.1007/s12028-009-9305-xGoogle Scholar

R. A. De Blasiet al., “Cerebral and muscle oxygen saturation measurement by frequency-domain near-infra-red spectrometer,” Med. Biol. Eng. Comput. 33(2), 228–230 (1995).MBECDY0140-0118http://dx.doi.org/10.1007/BF02523048Google Scholar

B. J. SteinhoffG. HerrendorfC. Kurth, “Ictal near infrared spectroscopy in temporal lobe epilepsy: a pilot study,” Eur. J. Epilepsy 5(2), 97–101 (1996)1059-1311http://dx.doi.org/10.1016/S1059-1311(96)80101-4.Google Scholar

L. Spinelliet al., “Multimodality fNIRS-EEG, fMRI-EEG and TMS clinical study on cortical response during motor task in adult volunteers and epileptic patients with movement disorders,” in Biomedical Optics, OSA Technical Digest, p. BM4A.1, Optical Society of America, Washington, DC (2012).Google Scholar

A. H. Hielscheret al., “Frequency-domain optical tomographic imaging of arthritic finger joints,” IEEE Trans. Instrum. Meas. 30(10), 1725–1736 (2011).IEIMAO0018-9456http://dx.doi.org/10.1109/TMI.2011.2135374Google Scholar

J. HwangJ. C. Ramella-RomanR. Nordstrom, “Introduction: feature issue on phantoms for the performance evaluation and validation of optical medical imaging devices,” Biomed. Opt. Express 3(6), 1399–1403 (2012).BOEICL2156-7085http://dx.doi.org/10.1364/BOE.3.001399Google Scholar

A. Pifferiet al., “Performance assessment of photon migration instruments: the MEDPHOT protocol,” Appl. Opt. 44(11), 210414 (2005).APOPAI0003-6935http://dx.doi.org/10.1364/AO.44.002104Google Scholar

A. Wabnitzet al., “Assessment of basic instrumental performance of time-domain optical brain imagers,” Proc. SPIE 7896, 789602 (2011).PSISDG0277-786Xhttp://dx.doi.org/10.1117/12.874654Google Scholar

A. Wabnitzet al., “Performance assessment of time-domain optical brain imagers: the nEUROPt protocol,” in Biomedical Optics, OSA Technical Digest, p. BSu2A.4, Optical Society of America, Washington, DC (2012).Google Scholar

C. H. Schmitzet al., “Instrumentation and calibration protocol for imaging dynamic features in dense-scattering media by optical tomography,” Appl. Opt. 39(4), 6466–6486 (2000).APOPAI0003-6935http://dx.doi.org/10.1364/AO.39.006466Google Scholar

M. FirbankD. T. Delpy, “A design for a stable and reproducible phantom for use in near infra-red imaging and spectroscopy,” Phys. Med. Biol. 38, 847–853 (1993).PHMBA70031-9155http://dx.doi.org/10.1088/0031-9155/38/6/015Google Scholar

U. Sukowskiet al., “Preparation of solid phantoms with defined scattering and absorption properties for optical tomography,” Phys. Med. Biol. 41(9), 1823–1844 (1996).PHMBA70031-9155http://dx.doi.org/10.1088/0031-9155/41/9/017Google Scholar

R. Cubedduet al., “A solid tissue phantom for photon migration studies,” Phys. Med. Biol. 42(10), 1971–1979 (1997).PHMBA70031-9155http://dx.doi.org/10.1088/0031-9155/42/10/011Google Scholar

M. Lualdiet al., “A phantom with tissue-like optical properties in the visible and near infrared for use in photomedicine,” Lasers Surg. Med. 28(3), 237–243 (2001).LSMEDI0196-8092http://dx.doi.org/10.1002/(ISSN)1096-9101Google Scholar

J. C. Hebdenet al., “A soft deformable tissue-equivalent phantom for diffuse optical tomography,” Phys. Med. Biol. 51(21), 5581–5590 (2006).PHMBA70031-9155http://dx.doi.org/10.1088/0031-9155/51/21/013Google Scholar

P. Di NinniF. MartelliG. Zaccanti, “Intralipid: towards a diffusive reference standard for optical tissue phantoms,” Phys. Med. Biol. 56(2), N21–N28 (2011).PHMBA70031-9155http://dx.doi.org/10.1088/0031-9155/56/2/N01Google Scholar

A. E. Cerussiet al., “Tissue phantoms in multicenter clinical trials for diffuse optical technologies,” Biomed. Opt. Express 3(5), 966–971 (2012).BOEICL2156-7085http://dx.doi.org/10.1364/BOE.3.000966Google Scholar

G. Lamoucheet al., “Review of tissue simulating phantoms with controllable optical, mechanical and structural properties for use in optical coherence tomography,” Biomed. Opt. Express 3(6), 1381–1398 (2012).BOEICL2156-7085http://dx.doi.org/10.1364/BOE.3.001381Google Scholar

B. W. PogueM. S. Patterson, “Review of tissue simulating phantoms for optical spectroscopy, imaging and dosimetry,” J. Biomed. Opt. 11(4), 041102 (2006).JBOPFO1083-3668http://dx.doi.org/10.1117/1.2335429Google Scholar

F. MartelliG. Zaccanti, “Calibration of scattering and absorption properties of a liquid diffusive medium at NIR wavelengths. CW method,” Opt. Express 15(2), 486–500 (2007).OPEXFF1094-4087http://dx.doi.org/10.1364/OE.15.000486Google Scholar

L. Spinelliet al., “Calibration of scattering and absorption properties of a liquid diffusive medium at NIR wavelengths. Time-resolved method,” Opt. Express 15(11), 6589–6604 (2007).OPEXFF1094-4087http://dx.doi.org/10.1364/OE.15.006589Google Scholar

L. Spinelliet al., “Inter-laboratory comparison of optical properties performed on Intralipid and India Ink,” in Biomedical Optics, OSA Technical Digest, p. BW1A.6, Optical Society of America, Washington, DC (2012).Google Scholar

J. P. Bouchardet al., “Reference optical phantoms for diffuse optical spectroscopy. Part 1–Error analysis of a time resolved transmittance characterization method,” Opt. Express 18(11), 11495–11507 (2010).OPEXFF1094-4087http://dx.doi.org/10.1364/OE.18.011495Google Scholar

L. Spinelliet al., “Experimental test of a perturbation model for time-resolved imaging in diffusive media,” Appl. Opt. 42(16), 3145–3153 (2003).APOPAI0003-6935http://dx.doi.org/10.1364/AO.42.003145Google Scholar

A. Puszkaet al., “Time-domain reflectance diffuse optical tomography with Mellin-Laplace transform for experimental detection and depth localization of a single absorbing inclusion,” Biomed. Opt. Express 4(4), 569–583 (2013).BOEICL2156-7085http://dx.doi.org/10.1364/BOE.4.000569Google Scholar

D. A. Boaset al., “Scattering of diffuse photon density waves by spherical inhomogeneities within turbid media: analytic solution and application,” Proc. Natl. Acad. Sci. 91(11), 4887–4891 (1994).PMASAX0096-9206http://dx.doi.org/10.1073/pnas.91.11.4887Google Scholar

J. C. J. PaasschensG. W. Hooft, “Influence of boundaries on the imaging of objects in turbid media,” J. Opt. Soc. Am. A 15(7), 1797–1812 (1998).JOAOD60740-3232http://dx.doi.org/10.1364/JOSAA.15.001797Google Scholar

A. Takatsukiet al., “Absorber’s effect projected directly above improves spatial resolution in near infrared backscattered imaging,” Jpn. J. Physiol. 54(1), 79–86 (2004).JJPHAM0021-521Xhttp://dx.doi.org/10.2170/jjphysiol.54.79Google Scholar

A. Pifferiet al., “Time-resolved diffuse reflectance using small source-detector separation and fast single-photon gating,” Phys. Rev. Lett. 100(3), 138101 (2008).PRLTAO0031-9007http://dx.doi.org/10.1103/PhysRevLett.100.138101Google Scholar

Q. Zhaoet al., “Functional tomography using a time-gated ICCD camera,” Biomed. Opt. Express 2(3), 705–716 (2011).BOEICL2156-7085http://dx.doi.org/10.1364/BOE.2.000705Google Scholar

F. Martelliet al., Light Propagation Through Biological Tissue and Other Diffusive Media, SPIE Press, Bellingham, Washington (2010).Google Scholar

A. Sassaroliet al., “Monte Carlo procedure for investigating light propagation and imaging of highly scattering media,” Appl. Opt. 37(31), 7392–7400 (1998).APOPAI0003-6935http://dx.doi.org/10.1364/AO.37.007392Google Scholar

C. K. Hayakawaet al., “Perturbation Monte Carlo methods to solve inverse photon migration problems in heterogeneous tissues,” Opt. Lett. 26(17), 1335–1337 (2001).OPLEDP0146-9592http://dx.doi.org/10.1364/OL.26.001335Google Scholar

D. Continiet al., “Multi-channel time-resolved system for functional near infrared spectroscopy,” Opt. Express 14, 5418–5432 (2006).OPEXFF1094-4087http://dx.doi.org/10.1364/OE.14.005418Google Scholar

R. Reet al., “A compact time-resolved system for near infrared spectroscopy based on wavelength space multiplexing,” Rev. Sci. Instrum. 81(11), 113101 (2010).RSINAK0034-6748http://dx.doi.org/10.1063/1.3495957Google Scholar

A. SassaroliF. MartelliS. Fantini, “Perturbation theory for the diffusion equation by use of the moments of the generalized temporal point-spread function. I. Theory,” J. Opt. Soc. Am. A 23(9), 2105–2118 (2006).JOAOD60740-3232http://dx.doi.org/10.1364/JOSAA.23.002105Google Scholar

A. SassaroliF. MartelliS. Fantini, “Perturbation theory for the diffusion equation by use of the moments of the generalized temporal point-spread function. II. Continuous-wave results,” J. Opt. Soc. Am. A 23(9), 2119–2131 (2006).JOAOD60740-3232http://dx.doi.org/10.1364/JOSAA.23.002119Google Scholar

A. SassaroliF. MartelliS. Fantini, “Higher-order perturbation theory for the diffusion equation in heterogeneous media: application to layered and slab geometries,” Appl. Opt. 48(10), D62–D73 (2009).APOPAI0003-6935http://dx.doi.org/10.1364/AO.48.000D62Google Scholar

A. SassaroliF. MartelliS. Fantini, “Perturbation theory for the diffusion equation by use of the moments of the generalized temporal point-spread function. III. Frequency-domain and time-domain results,” J. Opt. Soc. Am. A 27(7), 1723–1742 (2010).JOAOD60740-3232http://dx.doi.org/10.1364/JOSAA.27.001723Google Scholar