1. Introduction

Accurate diagnosis of diseases, particularly cancers, requires more detailed, disease-specific and preferably non- or minimally-invasive imaging information. The morphological images, traditionally obtained from X-ray based imaging modalities, are becoming insufficient in many cases [1]. In the past several decades, increasing research efforts in combining multiple imaging modalities, particularly those providing complementary information, have formed the grounds and active frontiers of multi-modality imaging [2]. Dual- or triple-modality systems are becoming abundant in publications and gradually become part of the standard clinical practices. Among them, a successful example is the combined Positron Emission Tomography and Computed Tomography (PET/CT) [3]. In PET/CT, the low-resolution albeit disease-specific PET image is co-registered to a high-resolution 3D CT image providing the background tissue anatomy. By interpreting the functional PET and structural CT images simultaneously, one is able to show a diagnostic accuracy that is not achievable by each individual modality [4]. PET/MRI (magnetic resonance imaging) [5], SPECT(single photon emission computed tomography)/CT [6] and MicroCT-Bioluminescence [7,8] are among other examples of multi-modality imaging that are gaining popularity for human or animal disease model studies.

Diffuse optical tomography (DOT) is another emerging technique that can provide functional assessment about tissue status [9,10]. By measuring the scattered light at multiple wavelengths within the near-infrared frequency range, one can recover hemoglobin (oxygenated, HbO, or deoxygenated, HbR), water and lipid concentrations as well as the tissue scattering properties. These measured quantities were found to be clinically relevant and are particularly promising as biomarkers for breast cancer diagnosis [11–15]. However, similar to PET, standalone DOT images suffer from low resolution due to the diffusive nature of the photons inside the tissue and the subsequent ill-posedness in solving the inverse model. Combining DOT with a high-resolution structural-oriented imaging modality can be highly beneficial [2]. Various dual-modality breast imaging systems involving DOT have been developed in the past decade and yielded different levels of success. Zhu et al. [16] built a combined DOT/Ultrasound hand-held probe for breast tumor diagnosis. Ntziachristos et al. [17], Brooksby et al. [18] and Carpenter et al. [19] have reported efforts combining DOT with MRI by either simultaneous measurement or post-registration. At Massachusetts General Hospital, we have developed a combined DOT and digital breast tomosynthesis (DBT) system for breast cancer screening [20,21]. In the past 4 years, we have been focusing on the clinical evaluation of the system by measuring healthy or tumor-bearing breasts from nearly 200 volunteers; we have shown promising diagnostic efficiency for malignant tumors using the reconstructed optical parameters from this pilot clinical study [22].

In most multi-modality DOT studies as mentioned above, the structural imaging modalities were only used to provide the exterior boundary of the breast, or overlaid on the reconstructed DOT images and used for image interpretation. This appears to be an under-use of the available information as the internal tissue structures are not fully exploited. In fact, the under-used internal structure is the key ingredient that is missing to produce DOT images of higher spatial resolution. To incorporate this ingredient into the optical image construction, novel algorithms were proposed and tested by combining a spatial prior to “inform” the optical image recovery. The simplest approach to bind a spatial structure to the DOT reconstruction is to segment the structural image and represent nodes within the same tissue segment by a single set of optical parameters, i.e. the so-called “hard-prior” approach [18,23]. Several limitations were found for this approach. The implied piecewise-constant optical properties within each tissue type is likely inaccurate in realistic biological tissue. The highly over-determined reconstruction equation effectively cancels the information from most optical measurements, which are usually precious to get. The resulting images are also intolerable to segmentation error, as it is essentially identical to the segmentation itself. To balance the contributions from the structural prior and optical measurements, “softer” constraints, in the form of regularizations, have been investigated. Li et al. has applied different regularization coefficients for different tissue types based on a segmented structural image [24]. Brooksby et al. [23] and Yalavarthy et al. [25] used a regularization matrix, defined by a Laplacian operator for the nodes sharing the same tissue type. Although these algorithms allow certain variations from the structural prior, the underlying piecewise-constant assumption remains to be arguable in realistic tissues. Only recently, researchers begin to exploit more general prior-guided DOT reconstruction by treating a target domain as statistical mixtures of different tissue types [26,27].

Here, we present a structural-prior guided DOT reconstruction algorithm using a compositional segmentation. Different from the binary approach, in a compositional segmentation we assume each pixel/node in the structural image is a combination of two or more tissue components. Using a weighted-graph representation in the compositional space, we are able to construct a regularization matrix, being the modified Laplacian operator of the weighted-graph, and use it in the DOT image reconstruction. We investigate the control of the “softness” of the structural constraints by changing the weight of the regularization. To test this algorithm, we use a clinical data set collected in the past few years with both DOT and DBT measurements. With these clinical data, we are able to identify the potential improvement in image resolution and quality using the proposed algorithm. We also calculate the estimates for the optical properties for each component of the breast tissue, i.e. adipose and fibroglandular tissue, using the compositional priors and compare them to our previous results using prior-free reconstruction and expert segmentation.

In the Methods section of the paper, we first report the theory for compositional segmentation of a structural image and then the formulations of the regularization matrix. In the Results section, we use a subset of our healthy breasts measurements (N = 31) to establish an empirical correlation between the X-ray image intensities for adipose and fibroglandular tissue. Using this relationship, we apply the compositional segmentation to another 58 healthy breasts from bilateral measurements of 29 subjects and study the improvement in image resolution and accuracy. Our findings are summarized in the Discussion section.

2. Methods

The proposed algorithm includes two essential steps. In the first step, we decompose a structural image into a compositional map using a given set of components and their pre-defined characteristics. After the decomposition, each pixel or node in the structural image is represented by a vector recording the concentration of each component. In the second step, we build a regularization matrix based on the previously generated compositional maps and run a DOT image reconstruction to produce a solution informed by the structural constraints. These steps are explained in detail in the following sub-sections.

2.1 Compositional segmentation and applications to DBT images

We consider that the imaging target is composed of a finite number of components, where each component is independent of the others and has differentiable and additive contrast when measured in the structural imaging modality. We assume the image intensity of the target in the structural image is a linear combination of the separate contrast from each of its components, which is proportional to its concentration. We define a compositional vector at a given location r in the target as

where C_i(r) denotes the concentration (0≤C_i(r)≤1) of the i-th component at location r, Ω denotes the target domain and N_c is the total number of components. By definition, Eq. (1) also implies ${\displaystyle {\sum }_i{C}_i(r)}=1$ for $\forall r\in \Omega$. If we assume the imaging intensity produced by the components can be represented by a transfer function f(·) governed by the biophysics of the tissue, then we have

where I_s(r) is the image intensity in the measured structural image and f(P_i) is the intensity anticipated if the target at r is occupied by the pure i-th component (P _i). The image segmentation problem is now identified as the task to uncover C = {C_i(r)} from I_s(r) for $r\in \Omega$. If there are only two components, i.e. N_c = 2, from Eq. (2) and ${\displaystyle {\sum }_i{C}_i(r)}=1$, we are able to uniquely solve for C_i(r) if all f(P_i) are known. When there are more than 2 components, additional knowledge needs to be acquired in order to uniquely solve for C_i(r).

For a healthy breast, two major tissue types are concerned: fibroglandular (denoted with subscript “f”) and adipose (denoted with subscript “a”) tissue. From a 3D DBT image, it is possible to identify a small region where considered as pure adipose (or fibroglandular) tissue, from which, we can estimate f(P _a) (or f(P _f)), and the compositional vector C(r) = {C_f,(r), C_a(r)} can be calculated as

The above approach can be extended to other types of structural images, for example, anatomical MRI. One may use the Gaussian-Mixture-Model (GMM) [28] or other data-clustering algorithms to identify the compositional vectors inside the target domain [27]. For algorithm and simulation studies, one can also derive such a compositional map from a published segmented breast data archive [29].

When such a compositional map is available for a given target, for any reconstructed functional image within the same domain, we can compute the functional values for each “pure” component by a least-square estimation as

where matrix H is a horizontal concatenation of the vertical compositional vectors for all parameter nodes and H ⁺ = (H ^T H)⁻¹ is the Moore-Penrose pseudo-inversion of H; μ represents a vertical vector recording the reconstructed functional values at all parameter nodes.

2.2 Regularized image reconstruction

In this subsection, we illustrate the compositional-map-guided regularization using our dual-modality data processing pipeline. This approach is general and can be easily extended to other pipelines that share a similar reconstruction scheme.

In a prior-guided DOT image reconstruction, we first discretize the target domain by generating a pair of meshes, a fine one for the forward solution (referred to as the forward mesh), and a coarse one for the inversion (referred to as the parameter mesh) [30]. For a co-registered dual-modality experiment, we register the meshes to the structural image space. Subsequently, we map the nodes on the forward mesh to the structural image volume, and assign a compositional vector, {C_f(r),C_a(r)} for each forward node. Using the mapping between the forward and parameter meshes, we calculate the compositional vectors at each node on the parameter mesh by weighted averages.

In our particular pipeline, we compute the forward solutions using an FE solver and perform an iterative Gauss-Newton reconstruction to recover the 3D profiles of the optical properties on the parameter mesh. For each iteration, the general form of the update equation looks like [31]

where $J_k\in C^{M\times N}$ is the complex Jacobian matrix at the k-th iteration, M is the number of measurements and N is the number of parameters; λ is the regularization parameter and L is the regularization matrix; $\Delta {\mu }_k$is the update to the parameters; A is the forward matrix and $A{\mu }_{k-1}$ is the forward solution using the estimated parameters from the k-1 iteration. Note that matrix $(J_k^T{J}_k+\lambda L^{T}L)$ has dimension of NxN. In practice, N is typically substantially bigger than the number of measurements M. In this situation, an under-determined form [32] of the update equation is more computationally efficient:

where $\Theta =L^{T}L\in R^{NxN}$ and can be pre-computed for each parameter mesh, and I is an identify matrix. Very often, when there are multiple (S>1) optical properties associated with each node, for example, HbO, HbR, scattering amplitude and scattering power, the calculations in Eq. (6) in the full-matrix form can be slow as it involves the multiplications of large matrices. This can be accelerated by the block-diagonal nature of the system. We simply expand Eq. (6) by

where the sub-matrix $J_{i,k}$ represents the Jacobian of the i-th optical property at iteration k. Using the above block-diagonal formula leads to a significant saving in memory and computational time.

2.3. Compositional-map-guided regularization

A non-trivial regularization matrix L ($L\ne I$ ) encodes the prior knowledge about the solution. When L is defined as a spatial Laplacian or Helmholtz operator [23] based on the mesh connectivity, it essentially penalizes large gradients and ensures smoothness of the solution. In the binary “soft-prior” approach [25], it penalizes the differences between the nodes falling into the same piecewise-constant segment. Here, we further extend this notion by taking advantage of the refined compositional segmentation scheme. To illustrate the method, we first construct a weighted graph in the composition space. In Fig. 1 , we plot all nodes in the parameter mesh using its compositional vector. In fact, it is a hyper-plane in the $R^{N_c-1}$ space defined by ${\displaystyle {\sum }_i{C}_i(r)}=1$. For each node, we connect it to the nodes sharing similar compositions within a predefined limit, i.e $\left|\right|C_i-C_j\left|\right|<\alpha N_c$ where $\left|\right|\cdot \left|\right|$ is the L₂ norm. Scalar $\alpha \in (0,1)$ is a regularization parameter and was heuristically set to 0.2 in our simulations. Smaller α values appear to produce sharper images Then we add an edge between the nodes satisfying the above condition and assign a weight to this edge proportional as $u_{i,j}=\alpha -\left|\right|C_i-C_j\left|\right|/N_c$ where N_c is the number of total components. To this end, we have constructed a weighted graph in the compositional space.

 

Fig. 1 Weighted graph representation of the target domain in the compositional space. Each solid dot is a node in the parameter mesh and u_i,j is the weight of the edge between nodes i and j.

Download Full Size | PDF

Based on [33], the (normalized) Laplacian operator, L = [l_i,j], for a weighted graph is defined by

where d_i is the “weighted” degree of the i-th node, defined by $d_i={\displaystyle {\sum }_{n,\left|\right|C_i-C_n\left|\right|<\alpha N_c}{u}_{i,n}}$. To avoid degeneracy and ensure the existence of ${\left(L^{T}L\right)}^{-1}$, we slightly modified L by

where β>1 is a constant. We choose β = 1.2 although other values may also be used. Using the L defined in Eq. (9) effectively penalizes the differences between nodes sharing similar compositions. One can combine it with the spatial Laplacian to impose additional smoothness constraints.

2.4. Clinical measurements for algorithm validation

We test the proposed regularization algorithm with a set of healthy subjects measured with our combined DOT/DBT imaging system. The details of the experiments and protocols are described in [21] and [22]. Both radio-frequency modulated (RF) and continuous-wave (CW) lasers at 685 nm and 830nm were used for the data acquisition. From these data, we recover the 3D profiles of total hemoglobin concentration (HbT = HbO + HbR), oxygen saturation (SO₂ = HbO/HbT), scattering amplitude and scattering power of the target breast. Water and lipids volume fractions are assumed to be constants at 23% and 58%, respectively, based on our previous studies [21].

3. Results

In this section, we show the recovered images for healthy breasts using the above described algorithm and compare the output with the results from unconstrained image reconstructions. We also estimated the optical properties for adipose and fibroglandular tissue and compare the findings with the previous results computed based on expert segmentation.

First, using a set of 31 DBT images of healthy breast as a learning set, we study the X-ray contrasts of different tissue types and seek possibilities for an automatic compositional segmentation approach. The intensity values in all DBT images were first rescaled to the range of 0-255. For each DBT image, we manually selected two regions-of-interest (ROI) corresponding to “pure” adipose and fibroglandular tissue. Each ROI is a 2D rectangular region about 1cm x 1cm in size from a selected image slice. From these ROIs, we estimated f(P_f) and f(P _a) values by averaging the DBT image intensities inside the ROIs. We plot the estimated f(P_f) and f(P _a) pairs for all selected breasts in Fig. 2 .

 

Fig. 2 X-ray intensity values (0-255) for “pure” adipose, f(P _a), and fibroglandular tissue, f(P_f), extracted from 31 DBT healthy breast images. The red-line denotes the linear regression line with the equation labeled in the figure.

Download Full Size | PDF

From Fig. 2, we assume a linear relationship between f(P_f) and f(P _a), i.e.

and the coefficients a = 1.226 and b = 13.76 are determined by linear regression from the plot in Fig. 2. The R ² of the linear fitting is 0.92. Thus, to compute the compositional map from Eq. (3), we only need to select a single ROI for adipose (or fibroglandular) tissue and compute the other from Eq. (10).

Using the above correlation, we apply the compositional segmentation for another 29 pairs of bilateral DOT/DBT measurements. In this case, we only manually select a single ROI corresponding to the “pure” adipose tissue, estimate f(P _a) from the ROI and calculate the value for f(P_f) from Eq. (10). The compositional maps of the 58 breasts are then obtained by Eq. (3).

To use the compositional maps in our FE-based image reconstructions, we first generated forward/parameter meshes using the “iso2mesh” mesh generator [34]. In all reconstructions, we first estimated a set of bulk optical properties which served as the homogeneous initial guess for the full image reconstruction. We used both RF and CW measurements at 685 nm and 830 nm in all the reconstructions. A multi-spectral algorithm [35,21] was used to estimate HbO, HbR and scattering properties from multi-wavelength measurements simultaneously. All image reconstructions ran for 5 Gauss-Newton iterations. All calculations were performed on a single core of an Intel Xeon E5530 (2.4 GHz) CPU.

In Fig. 3 , we show sample reconstructed HbT images (in μM) from the right breast of a 45-year-old healthy volunteer. The node numbers for forward and parameter meshes are 14824 and 1934, respectively. We solved for the forward solutions at 13/6 RF source/detector and 26/19 CW source/detector locations at two wavelengths and built the Jacobian matrix using the adjoint method [36]. The DBT image and the compositional map slice for fibroglandular tissue are shown in Figs. 3 (a) and (b). From (g) to (j), we show the image slices reconstructed with the compositional-prior-guided algorithm with different λ values, (λ = 2.5 to 0.039 with a multiplicative step factor of 0.25). In (c), we show the reconstruction output from the “binary soft-prior” algorithm (the binary segmentation of fibroglandular tissue was achieved by thresholding at C_f >0.5) and (d) the result without structural prior (i.e. L = I). All images were extracted at a horizontal plane z = 2.6 cm from the bottom of the compressed breast. The prior-guided image reconstructions took about 102 seconds per iteration, while that for the prior-free reconstruction is around 73 seconds. The relative residuals (normalized by the residual from a common initial guess) for the prior-guided reconstruction range from 0.27 (λ = 0.039) to 0.34 (λ = 2.5); that for the prior-free reconstruction is 0.32.

 

Fig. 3 Sample reconstructed images for a healthy breast: (a) the DBT image, (b) the fibroglandular composition map (C_f), and the reconstructed HbT image slices using (c) binary “soft-prior”, (d) prior-free and (g)-(j) compositional-prior-guided recon. with different regularizations. We also show reconstructed SO₂ images (e) without and (k) with the priors, similarly for μ’_s at 830nm (f) without and (l) with priors. The thin lines are segmentations using ITK-SNAP [37].

Download Full Size | PDF

Notice that the compositional map in Fig. 3(b) contains an edge artifact due to the application of an edge-enhancing algorithm in the DBT image processing. We deliberately choose this case as we want to see how the algorithm performs when the compositional map contains errors.

Based on Fig. 3, we chose λ around 1 as the default regularization parameter. We show SO₂ and μ’_s (at 830nm) for the same breast using compositional priors (λ = 1); in comparison, those without priors are shown in Figs. 3(e) and 3(f), respectively. Then we ran reconstructions for 58 healthy breasts. Similar to the previous cases, we ran 5 Gauss-Newton iterations to recover 3D volumetric images for HbO, HbR, scattering amplitude and power, respectively. Images from 4 additional healthy subjects are shown in Fig. 4 .

 

Fig. 4 Sample reconstructed images from 4 additional healthy breasts: the DBT images, the reconstructed HbT images using compositional-priors and prior-free reconstructions are shown in the 1st to the 3rd columns, respectively. The corresponding SO₂ and μ’_s (830nm) images reconstructed using the priors are shown as the 4th and 5th columns for each case. The dashed vertical lines mark the extent of the optical source/detector coverage.

Download Full Size | PDF

From the recovered images, we estimate the optical properties by Eq. (4) for adipose and fibroglandular tissue. In Figs. 5 (a)-(c) , we plot the estimated adipose HbT, SO₂ and reduced scattering coefficient (μ _s’) at 830nm for left vs. right breasts (circles). We also overlapped the results estimated by the prior-free algorithm using manually created ROIs [22] (dots). We also used these same manually generated ROIs to extract the optical properties for adipose and fibroglandular tissue from the prior-guided images. These were plotted as “pluses” in Fig. 5. In Table 1 , we summarize the correlation coefficients between the optical properties of the adipose tissue between left and right breasts, as well as the parameter values recovered from the 3 approaches. Those for the fibroglandular tissue are similar (not shown).

 

Fig. 5 Bilateral correlations between the adipose tissue for left and right breasts from 29 healthy subjects: (a) HbT, (b) SO₂ and (c) scattering coefficients μ _s’.

Download Full Size | PDF

Table 1. Correlation coefficients and the corresponding P-values between the bilateral adipose tissue properties extracted from 29 pairs of bilateral measurements are reported in columns entitled “R (P)”. We also summarize the means and standard deviations of the adipose physiological parameters estimated from all 58 breasts.

View Table

4. Discussion

The reconstructed images using the compositional priors in Figs. 3 (g)-(j) are very encouraging and exhibit several merits over the conventional image reconstructions. First of all, for most of the images, we can clearly identify the presence of the compositional map patterns in the recovered optical images, indicating that it has indeed informed the image reconstruction. Another finding from Figs. 3 is the robustness of the algorithm with respect to segmentation error: the edge artifacts in Fig. 3(b) are effectively removed in (g)-(j) with the incorporation of the optical data. In comparison, the result from the “binary soft-prior” approach seems more sensitive to segmentation error. Comparing the images from (g)-(j) and (d), we find that by adjusting the regularization parameter λ we are able to tune the strength of the priors: with an increasing λ, the reconstructed images blend toward the compositional map as we are giving more weight to the priors; with a small λ, on the other hand, the result approaches that of the prior-free reconstruction. By comparing (g)-(j) to (c), we find that the compositional prior approach gives more spatial details than the “binary” soft-prior approach, which is essentially a smoothed version of the binary segmentation. The images for SO₂ and μ _s’ show similar improvement in resolution. In (k), the SO₂ for the adipose tissue is generally greater than that of the fibroglandular tissue. This is consistent with our previous findings [21]. Additional images from 4 other subjects show similar improvement as we observed in Fig. 3, except for a few localized SO₂ contrasts in (i) and (s), which may be a result of coupling differences from the nearby optodes.

As we discussed in [22], the bilateral correlations serve as a metric to test the robustness of the system. From Fig. 5 and Table 1, we found that the optical properties estimated by the compositional-map-guided reconstructions generally match those estimated from prior-free reconstruction with expert segmentations. Particularly, the least-square estimates using the compositional map (circles) show a great similarity to those using the same image but with ROIs (pluses). These results render the prior-guided approach a very attractive choice as it is almost entirely automated, while the conventional method requires a manual segmentation to create the adipose and fibroglandular ROIs. Moreover, the compositional-prior-guided approach offers more objectivity than expert segmentation as the latter approach depends on the experiences of the operator. We also noticed that the standard deviations for HbT and SO2 values became larger when using priors in the reconstructions. This may be potentially a result of different regularization mechanisms in these methods and deserves further investigation. The optical properties estimated using priors are generally consistent with literature findings [11,18,21–23].

In the next stage of the study, we will further investigate the use of this approach for accurate quantification of breast cancers. To add cancer as the third component in the segmentation, we can simply segment the tumor and superimpose it on the adipose/fibroglandular tissue compositional maps. However, in cases where the size/location of the tumor is unknown, the application of this scheme can be more challenging. We will consider a two-step reconstruction where in the first step, we treat the breast as a healthy breast and reconstruct optical images as described here; in the second step, we use the previous result as the initial guess and reduce the strengths of the prior in order to recover the contrast from the tumor. Alternatively, we can use the HbT contrast from a prior-free reconstruction to estimate the compositional map for tumor tissue as we know tumors usually have high HbT values [11,12,22]. Combining the optically derived tumor tissue composition with the DBT-derived adipose/fibroglandular compositions, we can run the prior-guided reconstruction. In either case, we will use over 50 available tumor breast measurements in our current data set to assess the improvement using the statistical test similar to our previous work [22]. In addition, we will also explore strategies for optimal selection of the regularization parameters, λ, α and β, which are currently determined in a heuristic fashion. A deterministic selection scheme is expected to further improve the robustness of the algorithm.