1. Introduction

Non-invasive medical imaging is needed for early detection of breast cancer. Presently the “gold standard” for screening is x-ray mammography. Other imaging modalities, like optical mammography, are investigated for their diagnostic ability, to improve acceptance by not using ionizing radiation, and to reduce costs. When using continuous light the fundamental problem of optical mammography is the inability to separate scattering and absorption features resulting in blurred images [1]. In the past decade frequency domain and time-resolved techniques have been developed, to overcome this limitation [2–16]. Optical mammography has the potential to derive relevant physiological information, such as hemoglobin concentration, blood oxygenation, water and fat content, from the spatial distributions of the absorption coefficients measured at various selected near-infrared wavelengths.

In tomographic arrangements of source and detector positions the optical measurements can be analyzed for an approximate reconstruction of the said spatial distributions [17–19]. Alternatively, prior knowledge from other imaging modalities like e.g. MRI can be used for that purpose [20, 21]. In contrast, scanning a compressed breast in parallel-plane geometry yields projection images that can be compared for instance with x-ray mammograms without the need of time-consuming reconstruction algorithms [9,13]. Present models used to analyze such mammograms consider the lesion of interest as being imbedded in an otherwise optically homogeneous tissue slab [12,22–24]. In the present paper we employ perturbational calculations to derive the difference between the optical properties of the lesion from the surrounding tissue slab (Fig. 1).

 

Fig 1. Geometry for perturbation analysis. In our calculations a sphere is assumed although arbitrary shapes can be considered as well. The point model accounts only for the influence of a perturbation in the center of the sphere with radius 10 mm. The full calculation accounts for the whole volume V of the lesion (black area).

Download Full Size | PDF

Perturbation theory was reported by several authors as a valuable tool to derive optical properties of a given tumor [25–28]. Early work covered forward calculations that were compared with experimental data on phantoms and Monte-Carlo calculations, considering an infinite slab of tissue. Though the slab geometry is a reasonable approximation for a slightly compressed breast, specific corrections and assumptions have to be made to apply perturbation theory to derive maps of optical properties from in-vivo measurements. For instance, the finite size of the breast leads to lateral photon losses close to the edge of the breast, besides higher transmittance because of its reduced thickness. Such edge effects influence the measurement in a substantial part of the optical mammogram. The edge corrections introduced in the following will not only allow to generate maps of optical properties as a complete mammogram, but will also provide guidelines to modify the data analysis for calculating optical properties of lesions in the peripheral part of the optical mammogram. Our correction method is similar to that already tested successfully on the homogeneous slab model [9]. In addition, a procedure to choose a reference area in the mammogram for deriving effective background optical properties is described. Time-domain optical mammograms analyzed in the following were recorded by a laser-pulse scanning mammograph built at the PTB Berlin. The instrument was operated at two wavelengths, i.e. 785 nm and 670 nm. Details of this apparatus were published elsewhere [9].

Recently, a modified perturbation approach, based on the application of “Padé approximants”, was reported to interpret in-vivo optical data [22,23]. It was the aim of the authors of those papers to empirically extend the linear perturbation method to large variations in the optical properties, particularly for lesions with strongly increased absorption. In the present article we will not use Padé approximants although the edge corrections procedures outlined in this paper can be applied to this method as well.

2.Linear Perturbation Analysis of Optical Mammograms and Tumor Optical Properties

2.1 Linear Perturbation Theory

The linear perturbation model of Hebden and Arridge [26, 29] was used as forward model to analyze measured point spread functions, i.e. distributions of times of flight (DTOF) of photons through the slightly compressed breast, measured at a large number of scan positions.

The situation assumed in the perturbation model is sketched in Fig. 1. The tissue is considered as an infinitely extended partially homogeneous slab of thickness d(x,y), background absorption coefficient μ _a0, and photon diffusion coefficient D ₀ = 1/3μ′_s0 (μ′_s0 is the reduced scattering coefficient of the background). The inclusion of a homogeneous perturbation of volume V modifies local optical properties as μ_a = μ _a0 + δμ_a and D = D ₀ + δD. In linear perturbation theory the total transmitted time-resolved photon flux T_tot (t) at the site of the inhomogeneity is calculated as sum of the time-resolved transmitted flux of the bulk and of contributions from the inhomogeneity according to

where the so called shape functions are calculated from integrals over the volume of the perturbation

Explicit formulas for the integral kernels K^a,D are given in Ref. [28]. We emphasize that the integral kernels depend solely on μ _a0 and μ′_s0 of the homogeneous bulk, on the position $\overrightarrow{r}$ (relative to the source position $\overrightarrow{r}$ _s) and size (shape) of the impurity, as well as on the detector-source separation $\overrightarrow{r}$ _d-$\overrightarrow{r}$ _s. The in-vivo measurements analyzed here, were taken without lateral offset between source and detector, i.e. the dependence on the detector-source separation reduces to the dependence on breast thickness d(x, y).

2.2 Point Model for Generation of Absorption and Scattering Mammograms

Calculation of the shape functions (2) is time consuming. For a quick initial data analysis to generate optical absorption and scattering mammograms a simplified approximation of the exact integrals suffices, which will be called “point model” in the following. We note in passing that for a quantitative analysis of tumor optical properties the point model was not applied but the integrals in Eq. (2) were evaluated exactly.

In the point model the shape functions are calculated as the product of the integral kernel taken at the middle between source and detector and the volume of a test sphere with radius r¯ = 15 mm:

The exact numerical value of r¯ is not relevant, as the point model is used merely for the initial qualitative data analysis step as detailed below. It should be emphasized, however, that the value δμ_a scales inversely with the volume of the test sphere in the point model. Thus more accurate tumor optical properties can be obtained by the point model using the actual tumor radius rather than r¯.

2.3 Reference Area and Reference Point Spread Function for Background Optical Properties

For deriving the background optical properties μ _a0 and μ′_s0 a reference area Ω _ref which was assumed to be representative of the homogenous background was selected in each projection optical mammogram. To this end distributions of times of flight N(x, y, t) recorded at scan positions (x,y) were analyzed to find their first moment 〈t(x,y)〉. The median value 〈t〉_median of the first moments 〈t(x, y)〉 is used to define the limits of the reference area Ω_ref ,comprising all scan positions at which 〈t(x, y)〉 is larger than 〈t〉_median. Since the first moment 〈t(x, y)〉 is known to scale with the true thickness d(x, y) of the breast tissue at the scan position (x,y) [9], the reference area corresponds to that part of the optical mammogram where the breast filled most of the gap d ₀ between both compression (glass) plates of the optical mammograph. Background optical properties μ _a0 and μ′_s0 were derived from a reference point spread function N_ave (t), obtained by averaging all distributions of times of flight N(x,y,t) within the reference area Ω_ref. Subsequently, N_ave (t) was normalized to unit amplitude deriving the normalization factor $N_{\mathit{\text{exp}}}^{\mathit{\text{norm}}}$ = max[N_ave (t]

2.4 Background Optical Properties

For deriving the background absorption coefficient μ _a0 and reduced scattering coefficient μ′_s0 we calculated the time-resolved transmitted flux T ₀(d ₀,t) of a homogeneous infinite tissue slab of thickness d ₀ with extrapolated boundary conditions. Subsequently T ₀(d ₀,t) was convolved with the instrumental response function R(t) and normalized to unit amplitude deriving the normalization factor $N_{\mathit{\text{theo}}}^{\mathit{\text{norm}}}$ (d ₀) = max [T ₀(d ₀,t)*R(t)], (where * indicates convolution). The resulting curve was chosen as theoretical reference of normalized photon flux through a homogeneous tissue slab of thickness d ₀. By fitting this reference curve to the normalized reference point spread function N_ave (t) the desired optical coefficients μ _a0 and μ′_s0 were deduced.

2.5 Absorption and Scattering Mammograms from Time-Domain Data

In order to generate absorption and scattering mammogams, measured DTOFs N(x, y, t) were analyzed for absorption coefficients μ_a (x,y) = μ _a0 +δμ_a (x,y) and reduced scattering coefficients μ′_s (x, y) = μ′_s0 + δμ′_s (x, y) using the point model of linear perturbation theory. To this end all experimental curves N(x, y, t) were scaled by 1/$N_{\text{exp}}^{\mathit{\text{norm}}}$, independent of the real thickness d(x, y) of the breast. Whereas measured distributions of times of flight N(x, y, t) are normalized to approximately unit amplitude for scan positions within Ω_ref, distributions of times of flight of photons at locations (x,y) close to the edges of the breast and normalized in this way will exhibit amplitudes significantly larger than 1, because of the reduced thickness of the breast and hence increased transmittance.

In order to calculate the normalized, convolved total transmitted photon flux within the reference area Ω_ref the shape functions ${\stackrel{\mathit{~}}{f}}_{\mathit{\text{shape}}}^{a\mathit{,}D}$(d ₀, t) are evaluated for a breast thickness d ₀, subsequently convolved with the instrumental response R(t) and normalized by $N_{\mathit{\text{theo}}}^{\mathit{\text{norm}}}$ (d ₀). According to Eq.(1) the results weighted by δμ_a (x,y) and δD(x,y) are then added to the normalized, convolved homogeneous flux. For each scan position (x,y) ∉ Ω_ref we analyzed the first moments 〈t(x,y)〉 to estimate the thickness of breast designated by d^FM (x,y). The area of the mammogram outside the reference area Ω_ref where the breast is no longer in contact with both glass plates is then divided into a certain number of zones having approximately equal thickness $d_n^{\mathit{\text{FM}}}$ where n is the index of the zone. This reduces the continuous set d^FM (x,y) to discrete thickness values $d_n^{\mathit{\text{FM}}}$ . Subsequently the homogeneous photon flux T ₀($d_n^{\mathit{\text{FM}}}$ ,t) was calculated using the background optical properties (μ _a0 and μ′_s0) but taking the reduced breast thickness $d_n^{\mathit{\text{FM}}}$ < d ₀ into account. The resulting curve was convolved with the instrumental response and normalized obtaining the new normalization factor $N_{\mathit{\text{theo}}}^{\mathit{\text{norm}}}$ ($d_n^{\mathit{\text{FM}}}$ ) = max⌊T ₀($d_n^{\mathit{\text{FM}}}$ , t)* R(t)⌋.

For each zone shape functions ${\stackrel{\mathit{~}}{f}}_{\mathit{\text{shape}}}^{a\mathit{,}D}$ ($d_n^{\mathit{\text{FM}}}$ , t) were calculated for the appropriate thickness, subsequently convolved with the instrumental response and normalized by $N_{\mathit{\text{theo}}}^{\mathit{\text{norm}}}$ ($d_n^{\mathit{\text{FM}}}$ ), in analogy to the procedure described above for scan positions located within the reference area. The weighted results were added to the homogeneous contribution T ₀($d_n^{\mathit{\text{FM}}}$ , t)*R(t)/ $N_{\mathit{\text{theo}}}^{\mathit{\text{norm}}}$ ($d_n^{\mathit{\text{FM}}}$ ) to yield the total transmitted flux, normalized to the amplitude of the homogeneous reference and convolved with the instrumental response function. More generally

where d = d ₀ for (x,y) ∊ Ω_ref and d = $d_n^{\mathit{\text{FM}}}$ (x,y) falling into the n^th zone, and the star indicates the convolution.

Generally, scanning optical mammograms exhibit “edge effects” to various degrees. Such edge effects are artifacts caused by the reduced thickness of the breast close to the edge and by lateral photon losses. Actually, the most important correction for “edge effects” is to use the thickness d^FM (x,y) or $d_n^{\mathit{\text{FM}}}$ rather than constant thickness d ₀. However, if mammograms were to be extended further to the edge two additional corrections were applied, which were also used previously with the homogeneous model [9].

The first correction amounts to an amplitude scaling of the normalized convolved total transmitted photon flux $N_{\mathit{\text{fit}}}^{\mathit{\text{theo}}}$ (x,y,t;d) by the exponential factor exp(W η(d ₀-d)), where W = k ₁ μ _a0, and k ₁ amounts to about 15, typically but can be as large as k ₁ = 25. It follows from the discussion given above that measured distributions of times of flight N(x, y, t) normalized by $N_{\mathit{\text{theo}}}^{\mathit{\text{norm}}}$ (d ₀) exhibit amplitudes larger than 1 at locations (x,y) close to the edge of the breast whereas $N_{\mathit{\text{fit}}}^{\mathit{\text{theo}}}$ (x,y,t;d) (cf. Eq. (4)) is normalized essentially to unit amplitude even at those locations. The amplitude scaling factor exp(W η (d ₀ - d)) corrects for the different amplitude normalization. Furthermore, lateral photon losses have to be considered at the edge of the breast, causing the average path of transmitted photons to be longer than expected from the geometrical thickness d^FM (x, y) between source and detector. This effect is taken into account by correcting d^FM (x,y) according to the following empirical relation (“shape correction”)

and by using d_n rather than $d_n^{\mathit{\text{FM}}}$ in Eq.(4). The parameter S varies between 0 ≤ S ≤1, typically. In this way one obtains for the corrected, convolved, normalized total photon flux

where d = d ₀ for (x,y)∊ Ω_ref and d = d_n for (x,y) falling into the n^th zone. Finally, ${\mathrm{N}}_{\text{fit}}^{\mathit{\text{theo}}\mathit{,}\mathit{\text{corr}}}$(x,y,t;d) was fitted to the normalized experimental data N(x,y,t)/$N_{\text{exp}}^{\mathit{\text{norm}}}$to deduce μ_a (x, y) and μ′_s (x, y), ie. absorption and scattering mammograms.

Throughout this paper a standard simplex method was utilized for the curve fitting to minimize the residuals weighted to account for photon counting statistics [31].

Presently, an analytical model to quantify the effect of lateral phonon losses is not available for the geometry of a compressed breast. The empirical correction rules reported here were derived from measurements on more than 100 patients. We optimized the correction parameters W and S for each measured mammogram, using the same parameters to display absorption and scattering mammograms. The parameter S was adjusted to obtain a visually uniform μ′_s distribution pattern in the mammograms, free of steps and kinks caused by the separation into thickness zones. Fig. 2 shows a histogram of the optimal parameter for 107 mammograms. For those mammograms that needed a shape correction, a value of about S = 0.65 was found to be typical. This result agrees well with the range 0.63–0.7 reported in Ref [9] for the model of the homogeneous slab. This indicates that correction (5) is not caused by the limitations of the linear perturbation or point model but results from edge effects influencing the light propagation in breast tissue. The values S =0.65 and k ₁ = 15 can be used in general for the mammograms recorded within the collection of patients in our studies. Further improvements can be achieved in some cases using freely adjustable parameters S and W.

 

Fig. 2. Number of mammograms corrected by a shape correction factor S according to Eq. (5). Most mammograms do not need this correction (S = 1), since the scan did not extend far enough to the edge of the breast.

Download Full Size | PDF

We note in passing that we tried also to choose the homogeneous reference for every zone separately, by summing and averaging DTOF curves over all points in the zone. This choice turned out to be advantageous for some selected patients, but in general the reconstruction of optical mammograms was poor due to strong intensity variations in the narrow zones of constant thickness.

2.6 Tumor Optical Properties

In the absorption mammograms obtained by the point model of linear perturbation theory the location of the lesions was determined by comparison with other imaging modalities (mostly x-ray mammography). This step is necessary since the absorption of tumor tissue dominates the absorption mammograms only in a limited number of cases [12], if no contrast agents are applied. At the selected tumor position the experimental DTOF of photons is fitted with the true perturbation model (1) and (2). Tumor volume was inferred from pathology. The depth position of the lesion was not known in general. For data analysis the lesion was assumed in the middle between source and detector. In order to check the validity of this assumption we repeated calculations by moving the perturbation closely to the surface. To our estimate the error produced by assuming a wrong depth position is less 20 %. Data analysis by linear perturbation theory is thus robust and does not depend heavily on the knowledge of the exact depth location of the tumor. To our experience the choice of reference area described above allows to deduce quite reasonable optical properties of tumors in most cases.

3. Analysis of in-vivo Measurements and Discussion.

3.1 Optical Mammograms from Linear Perturbation Method

The effect of edge correction on absorption mammograms is demonstrated in Fig. 3 and in Fig. 4. Taking into account the variable breast thickness obtained from the first moment calculation as discussed above the tumor lying in the upper part can easily be located in Fig. 3(a,c). In contrast, assuming a constant breast thickness d ₀, this feature is considerably masked by edge effects (Fig. 3(b,d)).

 

Fig. 3. Effect of thickness correction on absorption mammograms derived by perturbation theory (patient #75, invasive ductal carcinoma, λ = 785 nm); (a, b) craniocaudal projection, tumor at x = -0.25 cm and y = 7.25 cm, μ _a0 = 0.032 cm^-1; (a) with thickness correction, edge correction factor W = 0.49 cm^-1, S = 0.47, tumor absorption μ_a = 0.088 cm^-1 from exact calculation; (b) without edge thickness correction, (c, d) mediolateral projection, same tumor at x = -2.0 cm, y = 8.0 cm; μ _a0 = 0.032 cm^-1; (c) with thickness correction, edge correction factor W= 0.8 cm^-1, S = 1.07 , tumor absorption μ_a = 0.072 cm^-1 from exact calculation; (d) without edge correction. Absorption coefficients exceeding the gray scale are marked black.

Download Full Size | PDF

Figure 4 shows another example where amplitude scaling and shape correction (cf. Eqs. (5, 6)) is needed. Fig. 4(a) displays the mammogram obtained using the thickness correction only (cf. Eq.(4)). After applying amplitude scaling as described above the range of applicability for perturbation analysis is extended as shown in Fig. 4(b). The nominal position of the tumor inferred from radiological information is at the scan position x = 3 cm, y = 5.25 cm. Fig. 4(b) clearly indicates that additional corrections (cf. Eq.(5,6)) permit to detect tumors located close to the edge of breasts.

 

Fig. 4. Comparison of absorption mammograms (a, b) generated by perturbation calculations (point model) and imaging photon counts in the late-time window for absorption contrast (patient #38, λ = 785nm, craniocaudal projection, μ _a0 = 0.032 cm^-1). (a) map of the absorption coefficient without amplitude and shape corrections (cf.Eq.(4)); (b) same mammogram but with amplitude scaling and shape correction (W = 0.848 cm^-1, S = 0.825) applied, tumor at x = 3 cm, y = 5.25 cm, tumor absorption μ_a = 0.10 cm^-1 from exact calculation; (c) reciprocal number ($N_8^{-1}$) of photons in the late-time window in relative units.

Download Full Size | PDF

As was shown previously [9,12,13,24], a convenient and efficient way to achieve online imaging in scanning optical mammography is to display (inverse) photon counts in selected time windows. Fig. 4(c) illustrates the same mammogram based on inverse photon counts 1/N ₈ accumulated in the eighth out of ten time windows, a model free approach known to be very robust against edge effects. Such mammograms display essentially absorption contrast. Imaging based on absorption coefficients derived by the point model of linear perturbation theory need to be corrected for edge effects, but provide comparable information as can be seen from Fig. (4a,c). A systematic analysis of lesion contrast in both kinds of mammograms confirms this finding [32].

Scattering features can be mapped qualitatively displaying inverse photon counts of the early-time window 1/N ₁. Fig. 5 compares mammograms with a cyst displaying reduced scattering coefficients derived by the point model of linear perturbation theory and based on (inverse) photon counts in an early time window. The lower scattering in the cyst leads to the bright area in both images. However, since mammograms based on inverse photon counts in early time windows provide absorption contrast to some extend besides scattering contrast, there are additional features clearly visible in Fig. 5(b) which are absent in Fig. 5(a). Fig. 5(c) demonstrates that these features coincide with the absorption of blood vessels, that are not supposed to display substantially increased scattering.

We note that for a PC (with AMD Athlon XP Processor, 1.5 GHz ) the calculation of a complete mammogram, applying the point model of linear perturbation theory takes about 1 min. to 2 min. This time is acceptable for online data analysis, allowing this method to be integrated into optical mammographs to enhance separation of absorption and scattering in raw data imaging.

 

Fig. 5. Comparison of optical mammogram generated by perturbation calculations (point model) and time-window imaging for scattering contrast (patient #98, cyst at x = 2.5 cm and y = 2 cm, λ = 785 nm, craniocaudal projection, μ _a0 = 0.059 cm^-1, μ′_s0 = 11.7 cm^-1). (a) Map of the reduced scattering coefficient employing amplitude and shape corrections (W = 0.88 cm^-1, S = 0.63), reduced scattering coefficient μ_s ′ = 7.2 cm^-1 of cyst from exact perturbation calculation; (b) reciprocal number (arbitrary units) of photons in the early time window $N_1^{-1}$ showing artifacts from absorption features; (c) map of absorption coefficient calculated with point model as in (a).

Download Full Size | PDF

3.2 Optical Properties of Tumors

A comparative study of optical properties of breast like phantoms, derived by linear perturbation theory using Eq.(1,2) reveals systematic deviations compared with the results obtained from nonlinear models such as diffraction of photon density waves by spherical inhomogeneities, or finite element calculations [34]. Despite such deficiencies we have analyzed tumor optical properties by linear perturbation theory. Fig. 6 displays examples of the change in the absorption and reduced scattering coefficients of tumors compared to the corresponding bulk values. Tumor optical properties were derived using Eq. (1,2) and normalized to the bulk values. Cases shown were selected to give comparable tumor optical properties as obtained from the model of diffraction of diffuse photon density waves [32,33]. The absorption coefficients of the tumors investigated are always enhanced with respect to bulk values (Fig. 6(a)), a result well known from previous investigations [33, 35–37]. Error bars for absorption coefficients are not given, because linear perturbation analysis systematically underestimates an increased absorption for tumors [23]. To overcome such limitations higher order perturbation analysis has to be used. The reduced scattering coefficient may be enhanced or reduced relative to that of the bulk, a result also reported previously [33]. Several tumors were located close to the edges of breasts. In those cases amplitude and shape corrections using the procedures outlined above were applied for data analysis.

The error bars in Fig. 6(b) were estimated from a linear perturbation analysis of simulated data generated by finite element calculations for breast-like phantoms. The error bars do not include contributions originating from uncertainties on the actual size, shape and depth of the tumors.

 

Fig. 6. Distribution of tumor optical properties (linear perturbation theory, cf. Eq.(1,2), λ = 785 nm). Relative change (a) of absorption coefficient and (b) of reduced scattering coefficient with respect to background optical properties for a number of patients (n=18). No clear trend is observed for the scattering whereas absorption is always increased.

Download Full Size | PDF

4. Conclusion

An efficient method to generate optical mammograms in the frame of a linear perturbation approach of diffusion theory is presented. The procedure to correct for edge effects suggested extends the applicability of perturbation analysis to those parts of the optical mammogram, where the breast is no longer in contact with both compression plates. Linear perturbation theory is shown here to be a powerful method for imaging absorption and scattering features, and is robust against missing information on geometrical parameters of lesions. The derived maps of the reduced scattering coefficient are showing well resolved cyst patterns and no obvious cross talk with absorption features. Quantitative data analysis of more than 18 patients has been presented as well. Further work is needed to evaluate the limits of linear perturbation theory. It is expected that methods for edge correction outlined here for linear perturbation analysis of compressed breast optical mammography can be used also to extend the applicability of genuine nonlinear perturbation analysis and as well as the method of Padé approximants.