Open Access
Add to BibSonomy
Abstract
This work concerns a fluorescence optical projection tomography system for low scattering tissue, like lymph nodes, with angular-domain rejection of highly scattered photons. In this regime, filtered backprojection (FBP) image reconstruction has been shown to provide reasonable quality images, yet here a comparison of image quality between images obtained by FBP and iterative image reconstruction with a Monte Carlo generated system matrix, demonstrate measurable improvements with the iterative method. Through simulated and experimental phantoms, iterative algorithms consistently outperformed FBP in terms of contrast and spatial resolution. Moreover, when projection number was reduced, in order to reduce total imaging time, iterative reconstruction suppressed artifacts that hampered the performance of FBP reconstruction (structural similarity of the reconstructed images with “truth” was improved from 0.15 ± 1.2 × 10−3 to 0.66 ± 0.02); and although the system matrix was generated for homogenous optical properties, when heterogeneity (62.98 cm-1 variance in µs) was introduced to simulated phantoms, the results were still comparable (structural similarity homo: 0.67 ± 0.02 vs hetero: 0.66 ± 0.02).
© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. Introduction
Angular-domain fluorescence OPT is a method to reject scattered light in objects less than about 3 times the inverse of the reduced scattering coefficient of the tissue [1,2], and is a low cost solution to significantly improve spatial resolution of OPT to address the challenge of poor detection sensitivity in lymph node pathology [3,4]. While conventional OPT of lymph node sized tissues (∼5-10 mm diameter) requires optical clearing of the tissue to enable predominantly forward-travelling light propagation—and therefore simple back-projection reconstruction algorithms—the proposed method makes use of the low scattering properties of lymph nodes and angular domain imaging to reject highly scattered light. By reducing the detector angle of acceptance, only the preferentially straightest (ballistic and quasi-ballistic) traveling photons are collected. Although greater levels of scatter rejection can theoretically provide enhanced spatial resolution, there is a tradeoff with the number of photons detected. Through simulation [1,5] it was demonstrated that for the particular application of lymph node assessment (3-6 mm diameter samples with low scattering optical properties), strict angular restriction (NA = 0.005) and filtered back-projection (FBP) reconstruction were sufficient to detect and localize 200 µm “micrometastases” – the smallest clinically relevant [6]. Based on these simulations, a prototype system was developed following the same configuration. Experimentation with the system and fluorescent inclusions in porcine lymph nodes revealed similar promising findings [2,7].
While FBP can provide adequate results, it is of interest to explore other methods of reconstruction to determine the degree of improvement that more complex algorithms can offer. In addition, inherent to OPT, is a limited depth of focus. That is, depending on where the focal plane is positioned, different parts of each image contain focused and out-of-focus data; and although angular domain scatter rejection aids in reducing the numerical aperture—thereby increasing the depth of focus—it is not an idealized parallel beam system. Mesoscopic fluorescence tomography (MFT) [8] and fluorescence molecular tomography (FMT) [9] demonstrate how accurate modeling of photon propagation can be used to solve tomography as an inverse problem and provide high resolution images. Here, those techniques of model-based reconstruction employed in MFT and FMT were applied, and results were compared to standard FBP reconstruction used in conventional OPT, as well as the angular-domain system described above. Monte Carlo simulations were used to model the behavior of the imaging system and construct a system matrix. The aim in doing so was to be able to compare reconstruction performance with an iterative reconstruction algorithm to push the limits of spatial resolution, and to help in the design of the imaging system. While the methods used here are not novel, this work served as a first investigation as to whether imaging and reconstruction performance could be improved for the specific application of lymph node detection sensitivity.
2. Methods
2.1 Tomography as a linear problem
The inverse problem of tomographic image reconstruction can be modeled as a system of linear equations:
where g is an M-dimensional column vector representing the 2D measured image data; f is an N-dimensional column vector of the voxelized object where N has n × n elements; and H is an M × N system matrix that transforms the object data to image data. The measured data, g, is made up of m detectors for k different angles, such that M has length m × k. Elements of the system matrix, hij, each represent the contribution of the voxel j to detector element i. Thus, each row in H is the contribution of all voxels to a given detector element, and each column is the vectorized 2D image corresponding to a single voxel. The value or weight of each element is based on several aspects of the imaging system, which include source illumination geometry, detector response, aberrations from optical components, sample properties (geometry and optical scattering, attenuation, etc.) and specifications of the actual image acquisition [10].2.2 Forward model and Monte Carlo simulations
The detected fluorescence light field is given by Zhu et al. [9] as a pair of coupled equations:
Applying this theory, MC simulations were conducted using an open-source MC MATLAB program [14], which solves the radiative transfer equation based on the well-established mcxyz.c program developed by Jacques et al. [15]. Simulations were structured to match the configuration of the angular domain fluorescence OPT system (simplified schematic shown in Fig. 1(a) and more detailed description in [7]): the source was a 780 nm Gaussian LED emitter with 2.5 cm diameter spot size such that the entire sample was illuminated, and the detection sensitivity profiles of single detector elements within a sCMOS camera were modeled by a pencil beam “source,” (i.e., assuming angular restriction of all but the preferentially straightest photons). It was possible to model detectors in this way because of the isotropic emission of fluorescence, and actual detection of sufficient number of photons (especially with very small detector acceptance angle) would require long simulation time. While a point source illumination has the potential to improve resolution with more spatially nuanced system matrices, wide-field illumination was assumed for practical advantages (i.e., increased maximum output power while remaining below ANSI safety limits and reduced imaging time). As a preliminary test, 2.4 cm × 2.4 cm × 0.6 cm (L × W × H) volumes with lymph node matching optical properties throughout (µa = 0.3 cm-1, µs = 43 cm-1, g = 0.92, n = 1.4 [1]) were used [Fig. 2(a)]. The volume was divided into 500 × 500 × 125 voxels (L × W × H), 4× the size of the region of interest in the x and y dimensions to permit translation and rotation of the source-detector pairs during formation of the system matrix without loss of information. Source and single detector element simulations were run for 1 × 108 launched photon “packets” [∼90 min each on a laptop PC from early 2015 (MacBook Pro) with 64-bit macOS]. Although 3D photon propagation through the sample was modeled and volumetric reconstruction is possible, 2D analysis was carried out for simplicity since similar results would be expected [1].
Fig. 1. (a) Simplified angular domain optical projection tomography system. (b-d) Monte Carlo-generated normalized fluence rate maps plotted as a function of tissue thickness (z-axis, parallel to direction of illumination) and tissue radial distance (x-axis, perpendicular to the direction of illumination). (b) Source sensitivity profile. Representative (c) detector sensitivity profile and (d) source-detector pair sensitivity profile for a detector element located at x = 0 cm and y = 0 cm. Light is presented as traveling in the direction of the arrows – from left (source side, S) to right (detector side, D).
Fig. 2. Representative process of system matrix generation from Monte Carlo-produced source-detector sensitivity profiles shown for the (a) first and (b) last detector elements in a single array at one collection angle. Sensitivity profiles are discretized into N voxels of length n × n, and then vectorized, such that each row of the system matrix H is the contribution of all voxels to a single detector element m at collection angle k. The illustrated process is repeated for all source-detector pairs and collection angles. The red box indicates the row populated by the corresponding sensitivity profile.
Discretizing Eqs. (2) and (3) into N voxels replaces the integrals with summation, and the forward model can then be written in matrix form as shown above in Eq. (1), where H describes the probability for all source-detector pairs that at a specific position, a photon is emitted and detected. That is, vector f is then the unknown fluorophore distribution that we wish to solve for given the measurements g of detected fluorescence (see Section 2.4).
2.3 System matrix generation
The forward-adjoint MC method was used for system matrix generation [16]. The MC simulations described above were used to model the behavior and propagation of light detected by the angular domain OPT system, and more specifically, the probability of detecting fluorescence emission at each location of the object for each source-detector pair. For 2D reconstruction, a single slice (middle of the object) in the x-y plane, parallel to the optical axis and perpendicular to the axis of rotation was used. MC-generated representative source and detector sensitivity profiles are shown in Figs. 1(b) and 1(c), respectively, where light traveling through the system was simulated to propagate from left to right: from the source, through the medium, and then collection at the detector. The simulated detector profile was representative of a single detector element located at the center of the x-y plane, and the FOV was translated in one dimension to model an array of 125 detector elements. Source and detector sensitivity profiles were multiplied on an element-wise basis to construct individual source-detector pair probabilities of fluorescence emission detection [Fig. 1(d)] [16,17]. To construct the system matrix, H, individual source-detector sensitivity profiles were discretized into 125 × 125 voxels, then arranged into a vector, which then became one row of the matrix. Figure 2(a) illustrates the process for the first detector element (located at x = -0.3 cm, y = 0 cm) at collection angle 0°. For the next source-detector pair, the same source simulation was used, but the detector FOV was translated in the x-direction to model the sensitivity of the adjacent detector element. This was repeated for all k detector elements, continuing to populate H [Fig. 2(b)]. Then to model 72 different collection angles, the same discretization, vectorization, and ordering process was employed for each source-detector pair rotated over 360° in 5° intervals until the total system sensitivity matrix had H of size M × N where M consisted of k collection angles × m detector elements and N had n × n object voxels. Here, the size of H was 9,000 × 15,625 [(72 collection angles × 125 detector elements) × (125 x-axis voxels in the object × 125 y-axis voxels in the object)].
2.4 Image reconstruction
A non-regularized maximum-likelihood expectation-maximization (MLEM) algorithm was applied for image reconstruction. MLEM was chosen over other iterative reconstruction techniques because it is well-suited for emission tomography reconstruction where photon emission and detection follow Poisson noise distributions. Unlike MLEM, algebraic reconstruction techniques (ART) and their variations (MART: multiplicative ART, SIRT: simultaneous iterative reconstruction technique) do not consider the statistical nature of the data, and as such, are more susceptible to noise, and even more so as the number of iterations increase. While noise is also multiplicative for MLEM algorithms, the fact that each measurement (projection ray) in every iteration is updated in ART algorithms (as opposed to all pixels simultaneously in MLEM) causes noise amplification to be greater. In addition, MLEM is a well-established and proven technique for positron emission tomography and single photon emission computed tomography (similar linear inverse problems with Poisson noise), with simple implementation [10,18]. The iterative equation is given by [19]:
Typically, the initial estimate ${\hat{f}^0}$ is a uniform distribution; however, computation and convergence is inherently slow [9,20]. Here, the FBP reconstruction of the measured data was used as the initialization to investigate whether or not a usable solution could be achieved in fewer iterations. Since signal from fluorescence molecular imaging applications are usually localized to small regions, it was expected that such prior information would accelerate the algorithm. Reconstructions from the MLEM technique initialized with ones is referred to in this article as MLEMones, while results with the FBP reconstruction used as the first image estimate is designated MLEMFBP. FBP reconstruction was performed with the built-in iradon() function in MATLAB.
2.5 Test phantoms
2.5.1 Simulation phantom #1
To compare the performance of the two MLEM initialization techniques, a simple homogenous simulated phantom was created. A 0.28-mm-diameter inclusion of intensity 10 (arbitrary units) was simulated within a background of intensity 1, to represent embedded fluorescence (first image in Fig. 4). The system matrix (Section 2.3) was used as the forward model to simulate measured data; 1000 independent noise realizations were generated, simulating photon counting experiments, using Poisson statistical distributions.
2.5.2 Simulation phantom #2
To evaluate the ability of the algorithms to reconstruct objects of different size, a simulated phantom with four inclusions of increasing diameter (0.19, 0.38, 0.58, and 0.77 cm) were created [Fig. 3(c)]. Again, the objects were assigned a value of 10, while the background was set to 1. 1000 independent noise realizations were generated, simulating photon counting experiments, using Poisson statistical distributions. The same phantom geometry was used to assess reconstruction performance on heterogenous tissue. A new system matrix and simulated measurements were produced in the same way described above, except different lymph node tissue components were used to introduce heterogeneity. The associated geometry [Fig. 3(b)] and optical properties (Table 1) are found below. For simplicity, different tissue types were symmetric within the volume.
Fig. 3. (a) Homogenous and (b) heterogenous lymph node tissue volume geometry used in Monte Carlo simulations. Distinct tissue types are indicated by different colors with corresponding optical properties summarized in Table 1. (c) Variable-size inclusion phantom to test reconstruction algorithms.
2.5.3 Physical phantom #1
Solid lymph node-matching resin phantoms were used to test the iterative algorithm on experimental data, and to compare reconstruction quality to FBP alone. The phantom investigated was 6 mm in diameter and included a ∼1.15 mm inner diameter glass capillary tube filled with 10 µM of IRDye-800CW at the half the radius. The number of collected projections was also reduced from 72 to 36, 18, and 9, to investigate the robustness of each reconstruction technique against missing data. The same 72 projection data set was used, but it was resampled to match the number of desired views (e.g. every other projection was used for the 36-view data, every fourth projection for the 18-view data, etc.).
2.5.4 Tissue phantom #1
Freshly excised porcine lymph nodes (∼10 mm in diameter) were used to test the utility of the iterative algorithms on biological samples where optical properties are heterogenous and structure is more complex. To model metastases, the nodes were implanted as close to the center as possible, with ∼200 µm diameter breast cancer cell spheroids (MDA-MB-231), which were externally stained with 100 nM of fluorescent dye IRDye 800CW (LI-COR Biosciences).
2.6 Evaluation of image quality
Reconstruction performance was quantified with the following figures-of-merit: (1) region of interest (ROI) bias and variance, (2) mean-squared error (MSE) of the entire image, (3) contrast-to-noise ratio (CNR), (4) full width at half maximum (FWHM) of ROIs, and (5) structural similarity (SSIM). ROIs were defined as the “embedded fluorescence” regions in the phantoms. Ten noise realizations were performed for bias-variance analysis. Bias, Bi, and variance, si2, respectively, were calculated as [23]:
and where $I_i^{[r]}$, the relative error for noise realization r, is calculated as: and ${\bar{I}_i}$ is the mean relative error defined by: with Nr equal to the number of noise realizations. The means of Bi and ${s_i}^2$ were calculated for the ROI of each reconstruction iteration. For a given noise realization, r, the value of voxel i from the estimated image is denoted by $\hat{f}_i^{[r]}$, and ftrue,i is the true value from the corresponding voxel.The mean MSE between the estimated and whole true image was computed for each noise realization and reconstruction iteration as:
Owing to the discordance often found between values of MSE and actual perceived visual quality [24], SSIM was also determined for total images, and CNR and FWHM were measured for ROIs. For signals x and y, an SSIM index, $\textrm{SSIM}({\mathbf x},{\mathbf y}) = {[l({\mathbf x},{\mathbf y})]^\alpha } \cdot {[c({\mathbf x},{\mathbf y})]^\beta } \cdot {[s({\mathbf x},{\mathbf y})]^\gamma }$ given by Wang et al. [24], includes three components: luminance (l), contrast (c), and structure (s) that are defined respectively as:
CNR of the single inclusion phantom was calculated with the equation:
where Sobj and Sbkg, and $\sigma _{obj}^{}$ and $\sigma _{bkg}^{}$ are the mean signals and standard deviations of the object and background ROIs, respectively. Eight identically-sized regions located directly around the object ROI were used to compute a mean CNR value. Meanwhile, measures of FWHM were taken from the horizontal and vertical line profiles across the center of the inclusion.3. Results and discussion
3.1 Simulated phantom #1: Single inclusion reconstructions
The purpose of the single inclusion phantom studies was to compare the behavior of the two MLEM algorithms (first initialized with ones, and second with FBP images) as a function of iteration number. Results of both are shown in Fig. 4. From inspection, MLEMFBP appeared to converge toward a solution more rapidly than the uniformly initialized algorithm. The inclusion only became obvious in the 50th iteration image for MLEMones, whereas a discernable object was seen as early as the 5th iteration when initialized with FBP. These results were quantified through bias-variance analysis and MSE, summarized in Figs. 5(a) and 5(b).
Fig. 4. Comparison of iterative reconstructions of simulated phantom data with different initial image estimates for increasing number of iterations. The 0th iteration is the initialization image. MLEM: maximum-likelihood expectation-maximization, FBP: filtered backprojection.
Fig. 5. Quantitative results of MLEM reconstructions with increasing number of iterations for simulated data. MLEM: maximum-likelihood expectation-maximization, FBP: filtered backprojection, MSE: mean-squared error.
From observation of the plots, MLEMFBP had lower ROI bias than MLEMones for all iterations; and while both algorithms started with a relatively high bias (∼0.88), the FBP-initialized values dropped rather quickly in the first few iterations before following a slope similar to that of the ones-initialized [Fig. 5(a)]. This is consistent with what was seen in the image reconstructions, where rapid sharpening (decreased blur) occurred from the 5th to 20th iteration images using MLEMFBP, and after that, changes were rather constant (MSE percent difference between subsequent reconstructions from iterations 5 to 20: 0.18 ± 0.18% and iterations 21 to 100: 0.03 ± 5.4 × 10−3%). It should also be noted that similar levels of bias were achieved after approximately 30 iterations of MLEMFBP (0.81 ± 0.02) and 100 iterations of MLEMones (0.81 ± 0.01), again demonstrating the accelerated performance of the former approach. On the other hand, while variance increased for both algorithms with each iteration, FBP-initialized ROI variance grew more rapidly (∼3.5×), as evidenced by the steep curve in Fig. 5(b). This demonstrates a caveat of using the FBP reconstruction as the first guess – although it provides some desired prior information, FBP, especially with missing data (lower number of collected projections), is vulnerable to artifacts. Inherent to MLEM, high frequency components are amplified as iterations progress; thus, the noise artifacts from an initial FBP estimate are maintained or even increased. This bias-variance tradeoff is illustrated in Fig. 5(c) and based on that plot, MLEMones appeared to have better tradeoff compared to MLEMFBP. It is important to note here, however, that values of variance were about 3 orders of magnitude lower than the bias.
As a measure of the total image reconstruction quality, MSE was also computed. Figure 5(d) plots MSE versus iteration number, and the results demonstrate that FBP-initialized MLEM had lower MSE than the uniformly initialized algorithm after four iterations. The initial large error value of MLEMFBP can be attributed to artifacts carried over from the initial FBP estimate; however, this error was quickly diminished as iterations progressed. Interestingly, the MLEMones algorithm showed little change in MSE from the 1st to 100th iteration (variance of 1.98 × 10−5). This likely occurred because MSE is a measure of absolute error, and the true background value of the phantom, as well as the initial image estimate was ones.
Since resolution was also important in characterizing the algorithms, CNR and FWHM of the ROI were computed. Figure 6(a) shows greater CNR for all iterations (∼1.4×) when using FBP-initialized MLEM compared to ones-initialized MLEM. In addition, CNR remained rather constant for MLEMFBP as iteration number increased (variance of 2.8 × 103), suggesting that good contrast can be achieved with few iterations and maintained even with increasing variance. The plot of FWHM as a function of iteration number [Fig. 6(b)] illustrated the ability of MLEMFBP to provide more accurate spatial resolution than MLEMones after about 5 iterations. Similar to the MSE curve, the MLEMFBP algorithm corrected itself quickly after a large initial error and converged toward the expect FWHM value. MLEMones on the other hand, was slow to converge and was not able to reach the same level of approximation of the MLEMFBP algorithm within the first 100 iterations explored here – MLEMFBP achieved a minimum relative error to the expected FWHM of 0.57, while the same measure for MLEMones was two times greater at 1.14.
Fig. 6. (a) Contrast-to-noise ratio (CNR) and (b) full width at half maximum (FWHM) as a function of iteration number for different MLEM reconstructions of a simulated phantom. Shaded regions indicate standard deviation. MLEM: maximum-likelihood expectation-maximization, FBP: filtered backprojection.
3.2 Physical phantom #1: Reduced projection number reconstructions
Figure 7 displays the results of MLEM reconstructions of the lymph node-matching resin phantom embedded with fluorescence. The projection data was previously analyzed [2] using FBP, so the purpose here was to see if image quality could be improved with a more complex reconstruction algorithm. Results show that both methods of MLEM provided visually superior reconstructions compared to FBP alone (images in the 0th iteration column; rows 2 and 4). Image estimates generated from FBP exhibited ring and streak artifacts, as well as blurring around the inclusion edges – all of which became worse as the number of collected projections decreased. MLEM reconstructed estimates on the other hand, were able to suppress those image artifacts even as views were reduced. Between the two initialization approaches, MLEMFBP was able to achieve this more rapidly than MLEMones—about 5 iterations versus 50, respectively. Only at nine collected projections did the iterative algorithms fail to eliminate visible artifacts by the 100th iteration. It can be seen in both 100th iteration images that some parts of the ring and streak artifacts remained. In addition, the flaw of FBP-initialization was also revealed. Prominent streak artifacts near the inclusion were amplified, which caused the object to appear distorted and reduce shape fidelity. The ones-initialized estimate, however, did not have this problem, and the spherical inclusion shape was preserved. Nonetheless, overall the MLEM reconstructed images provided clearer image estimates compared to FBP reconstruction alone. Qualitatively, these differences were more subtle for the 72-projection data set; but for the severely reduced view of 9 total projections, the improvements were more noteworthy. These findings demonstrate the potential for acquisition of less projections, which can reduce overall tomographic imaging time—a favorable concept for clinical application.
Fig. 7. Ones- and FBP-initialized MLEM reconstructions of a fluorescence-embedded lymph node matching phantom for increasing number of iterations, and 72 and 9 projections. All images were auto-scaled independently for visualization. MLEM: maximum-likelihood expectation-maximization, FBP: filtered backprojection.
Because the ML criterion is slow to reach a solution, and variance is generally high at that point, an alternate premature stopping criterion was explored. The approach used by Ng et al. [25] was investigated where reconstructions were repeated until the maximal normalized pixel change in successive sinograms fell below a threshold. Based on preliminary simulations the percent change between consecutive reconstructions plateaued at approximately 1%, and so that threshold level was used for the following experiments.
3.3. Simulated phantom #2: homogenous and heterogenous tissue reconstruction
Results of the variable-sized four-inclusion phantom are shown in Fig. 8. The top panel [Fig. 8(a)], which presents the phantoms with homogenous optical properties, demonstrates the robustness of the MLEM algorithm against ring and streak artifacts that are present in the FBP-alone reconstructed images. The effects became more obvious as the number of projections was reduced; however, consistent with the results presented above, by nine projections, some of the artifacts from the FBP-initialization were compounded in progressive iterations. The shape of the largest inclusion was visibly distorted and it began to smear into the object next to it. While this was true, all four inclusions were still captured in the reconstructions. Conversely, the MLEMones image estimates provided rather poor delineation of the smallest inclusion even with all 72 projections, and by reduction to nine acquired views, the inclusion was completely blurred out. Therefore, it can be said that while image quality of FBP reconstruction alone suffers from artifacts specifically when data is missing, it is still able to reconstruct some small details.
Fig. 8. Filtered backprojection and iterative reconstructions of a simulated phantom with (a) homogenous and (b) heterogenous lymph node optical properties. Each case shows a plot of SSIM vs. iteration number for the different reconstruction methods and number of collected projections. MLEM: maximum-likelihood expectation-maximization, FBP: filtered backprojection, SSIM: structural similarity.
To quantify these results, mean SSIM was plotted against iteration number to provide an index of image quality more closely associated with visual perception [Fig. 8(a), right]. While metrics such as MSE (as used in section 3.1 above) can be used for a more direct quantitative comparison between methods, SSIM was deemed useful because the images will ultimately be evaluated by a human. For all iterations and number of projections, the MLEM algorithms had better SSIM index than FBP alone (an index of 1 indicates identical sets of data). Structural similarity was improved from 0.15 ± 1.4 × 10−3 with FBP alone to 0.67 ± 0.02 on average for MLEM initialized with both ones and FBP, over all acquired projection numbers. For higher number of projections (72 and 36), SSIM was greater for MLEMFBP compared to MLEMones after about 3-8 iterations. For less collected data (18 and 9 projections), FBP-initialized MLEM had lower SSIM than ones-initialized MLEM, and only appeared to meet it at around 80 iterations for the case of 18 projections, while never achieving the same performance for 9 projections . This was consistent with what was observed in the image reconstructions. With the exception of the 9-projection data set, FBP-initialized MLEM converged toward a solution with greater SSIM faster than ones-initialized MLEM, and was able to detect smaller inclusions with fewer projections, although noisier overall. In fact, visually the MLEMFBP estimates closely resembled the FBP-alone reconstructions, but without the image artifacts. Together these findings highlight the interplay that exists between noise suppression and image resolution, where the choice of reconstructive algorithm may depend on the clinical task. The superior performance of MLEMones at early iterations may be explained by the fact that the majority of the phantom is background space made up of intensity of one, and the initial estimate for the algorithm was unity.
Panel (b) of Fig. 8 displays the reconstructions of the phantom with heterogenous properties, and the results are very similar to that of the homogenous case. Again, the iterative algorithms were more robust against image artifacts and FBP initialization tended to provide visually superior contrast and resolution compared to ones-initializations (mean SSIM ± standard deviation over all projection numbers for FBP alone: 0.15 ± 6.45 × 10−4, MLEMones: 0.65 ± 0.03, MLEMFBP: 0.66 ± 0.02). These observations were supported by SSIM measures, shown in the plot on the right of Fig. 8(b). From inspection, the homogenous and heterogenous plots looked similar, and that was confirmed with comparable values of average SSIM for the MLEM algorithms (homogenous: 0.67 ± 0.02 vs. heterogenous: 0.66 ± 0.02). Visually, performance appeared better in the heterogenous reconstructions than the homogenous, which can be explained by the lower scattering properties of the heterogenous phantom. Nonetheless, the system matrix used in the iterative algorithm was not generated with various scattering properties, and the similar findings between the homogenous and heterogenous cases demonstrate the ability of the algorithm to overcome simplifications and inaccuracies in the system matrix.
3.3 Tissue phantom #1: porcine lymph node metastases model reconstruction
Figure 9 shows the reconstructions of a transverse slice through a lymph node implanted with fluorescent-stained cancer cell spheroids. As demonstrated here [Fig. 9(a)], and in our previous work [2], FBP was sufficient to reconstruct the fluorescent inclusions in ∼10 mm diameter lymph nodes. The MLEM reconstructions generated similar results where the ones-initialized method appeared slightly more blurry compared to the FBP-initialized approach. Similar to the simulated heterogeneous phantom above, the system matrix was not optimized for the lymph node sample because of its heterogeneity in optical properties, more complex structure and larger sample size (which may explain the loss in resolution); and yet, the MLEM algorithms were qualitatively comparable to that of the FBP. Standard lymph node pathology involves ∼2 mm-thick sectioning of the tissue such that <1% of entire node volumes are assessed and 30-60% of micrometastes go undetected. Comprehensive evaluation via tomography, and high resolution (200 µm sensitivity – the smallest clinically relevant) localization would help guide pathologists and eliminate blind gross-sectioning that contributes to high rates of false negatives. Thus, any improvements in resolution offered by MLEM reconstruction would be beneficial. While any significant differences in image performance among the methods tested is difficult to determine because of the nature of the metastases model and spheroid implantation (no true “gold standard” image, possibility of the spheroid breaking), the findings still support the potential for system-based reconstruction with the proposed angle-restricted imager to meet, if not improve the detection and localization capabilities of FBP alone. Therefore, selection of a reconstruction method with the goal of improving lymph node detection sensitivity requires evaluation with task-based metrics.
Fig. 9. Filtered backprojection and iterative slice reconstructions of fluorescence overlaid onto absorption of a porcine lymph node embedded with fluorescent-stained cancer cell spheroids. All images were scaled independently and thresholded for visualization, such that 90% of fluorescence signal above the background is shown. MLEM reconstructions were resampled to match that of the FBP reconstruction. Scale bars are 2 mm. MLEM: maximum-likelihood expectation-maximization, FBP: filtered backprojection
4. Conclusion
The work presented here demonstrated the feasibility of filtered backprojection (FBP) and iterative image reconstruction (with Monte Carlo generated system matrix) for the angular-domain fluorescence OPT system in the sub-diffuse photon propagation regime (> 2 scattering lengths, < 2 reduced scattering lengths [8]). Through simulated and experimental phantoms, iterative algorithms proved to consistently outperform FBP in terms of contrast and spatial resolution. Specifically, MLEM approaches where the FBP reconstruction was used as the initial image estimate instead of the uniform initialization, were able to achieve superior results in as few as five iterations. Moreover, when the number of collected projections was reduced, MLEM reconstruction was able to suppress ring and streak artifacts that hamper the performance of FBP reconstruction. These findings support the potential for reducing acquisition time utilizing limited angle tomography and therefore, more rapid imaging, which is favorable for clinical environments. In addition, although the system matrix was generated for very simple geometry and homogenous optical properties, when physiologically relevant heterogeneity was introduced to simulated phantoms, and porcine lymph node tissue samples were imaged, the results were still satisfactory. Thus, application of the proposed reconstruction methods to obtain images and improve image quality in low scattering biological tissue samples is promising.
Funding
National Institute of Biomedical Imaging and Bioengineering (EB02396); National Science Foundation (1653627); Illinois Institute of Technology (Nayar Prize); Pritzker Institute of Biomedical Science and Engineering.
Disclosures
The authors declare no conflicts of interest.
References
1. L. Sinha, F. Massanes, V. C. Torres, C. Li, K. M. Tichauer, and J. G. Brankov, “Comparison of time- and angular-domain scatter rejection in mesoscopic optical projection tomography: a simulation study,” Biomed. Opt. Express 10(2), 747–760 (2019). [CrossRef]
2. V. C. Torres, C. Li, Y. He, L. Sinha, G. Papavasiliou, H.A. Sattar, J. G. Brankov, and K. M. Tichauer, “Angular restriction fluorescence optical projection tomography to localize micrometastases in lymph nodes,” J. Biomed. Opt. 24(11), 1–4 (2019). [CrossRef]
3. S. M. Shiller, R. Weir, J. Pippen, M. Punar, and D. Savino, “The sensitivity and specificity of sentinel lymph node biopsy for breast cancer at Baylor University Medical Center at Dallas: a retrospective review of 488 cases,” Proc (Bayl Univ Med Cent) 24(2), 81–85 (2011). [CrossRef]
4. D. L. Weaver, “Pathology evaluation of sentinel lymph nodes in breast cancer: protocol recommendations and rationale,” Mod. Pathol. 23(S2), S26–S32 (2010). [CrossRef]
5. V. C. Torres, L. Sinha, C. Li, K. M. Tichauer, and J. G. Brankov, “Excised whole lymph node imaging for cancer staging with angular restriction dual fluorescent optical projection tomography,” in IEEE 16th ISBI (2019).
6. D. L. Weaver, U. P. Le, S. L. Dupuis, K. A. Weaver, S. P. Harlow, T. Ashikaga, and D. N. Krag, “Metastasis detection in sentinel lymph nodes: comparison of a limited widely spaced (NSABP protocol B-32) and a comprehensive narrowly spaced paraffin block sectioning strategy,” Am. J. Surg. Pathol. 33(11), 1583–1589 (2009). [CrossRef]
7. V. C. Torres, C. Li, W. Zhou, J. G. Brankov, and K. M. Tichauer, “Characterization of an angular domain fluorescence optical projection tomography system for mesoscopic lymph node imaging,” Appl. Opt. 60(1), 135–146 (2021). [CrossRef]
8. M. S. Ozturk, V. K. Lee, L. Zhao, G. Dai, and X. Intes, “Mesoscopic fluorescence molecular tomography of reporter genes in bioprinted thick tissue,” J. Biomed. Opt. 18(10), 100501 (2013). [CrossRef]
9. Y. Zhu, A. K. Jha, D. F. Wong, and A. Rahmim, “Image reconstruction in fluorescence molecular tomography with sparsity-initialized maximum-likelihood expectation maximization,” Biomed. Opt. Express 9(7), 3106–3121 (2018). [CrossRef]
10. D. S. Lalush and M. N. Wernick, “Iterative Image Reconstruction,” in Emission Tomography (Elsevier, 2004), ch. 21, sec. II. A., pp. 444–445.
11. F. Yang, R. Yao, M. Ozturk, D. Faulkner, Q. Qu, and X. Intes, “Improving mesoscopic fluorescence molecular tomography via preconditioning and regularization,” Biomed. Opt. Express 9(6), 2765–2778 (2018). [CrossRef]
12. F. Yang, D. Faulkner, R. Yao, M. S. Ozturk, Q. Qu, and X. Intes, “System configuration optimization for mesoscopic fluorescence molecular tomography,” Biomed. Opt. Express 10(11), 5660–5674 (2019). [CrossRef]
13. L. Wang and S. L. Jacques, “Hybrid model of Monte Carlo simulation and diffusion theory for light reflectance by turbid media,” J. Opt. Soc. Am. A 10(8), 1746–1752 (1993). [CrossRef]
14. D. Marti, R. N. Aasbjerg, P. E. Andersen, and A. K. Hansen, “MCmatlab: an open-source, user-friendly, MATLAB-integrated three-dimensional Monte Carlo light transport solver with heat diffusion and tissue damage,” J. Biomed. Opt. 23(12), 1–6 (2018). [CrossRef]
15. S. L. Jacques, T. Li, and S. Prahl, “mcxyz.c, a 3D Monte Carlo simulation of heterogeneous tissues,” omlc.org/software/mc/mcxyz.
16. J. Chen and X. Intes, “Comparison of Monte Carlo methods for fluorescence molecular tomography-computational efficiency,” Med. Phys. 38(10), 5788–5798 (2011). [CrossRef]
17. H. Dehghani, M. E. Eames, P. K. Yalavarthy, S. C. Davis, S. Srinivasan, C. M. Carpenter, B. W. Pogue, and K. D. Paulsen, “Near infrared optical tomography using NIRFAST: Algorithm for numerical model and image reconstruction,” Commun. Numer. Meth. Engng. 25(6), 711–732 (2009). [CrossRef]
18. N. Zeraatkar, A. Rahmim, S. Sarkar, and M. R. Ay, “Development and Evaluation of Image Reconstruction Algorithms for a Novel Desktop SPECT System,” Asia Ocean J. Nucl. Med. Biol. 5(2), 120–133 (2017). [CrossRef]
19. J. A. Fessler, “Statistical Image Reconstruction Methods,” in Handbook of Medical Imaging: Medical Image Processing and Analysis, vol. 2 (SPIE, 2000), ch. 1.9.1, sec. Emission Reconstruction, pp. 50–51.
20. D. S. Lalush and M. N. Wernick, “Iterative Image Reconstruction,” in Emission Tomography: Elsevier, 2004, ch. 21, pp. 443–472.
21. L. Scolaro, R. A. McLaughlin, B. R. Klyen, B. A. Wood, P. D. Robbins, C. M. Saunders, S. L. Jacques, and D. D. Sampson, “Parametric imaging of the local attenuation coefficient in human axillary lymph nodes assessed using optical coherence tomography,” Biomed. Opt. Express 3(2), 366–379 (2012). [CrossRef]
22. V. Torres, T. Wilson, C. Li, L. Sinha, J. Brankov, and K. Tichauer, Characterization of Lymph Node Optical Properties for Phantom Fabrication (SPIE, 2019).
23. T. Correia, M. Koch, A. Ale, V. Ntziachristos, and S. Arridge, “Patch-based anisotropic diffusion scheme for fluorescence diffuse optical tomography–part 2: image reconstruction,” Phys. Med. Biol. 61(4), 1452–1475 (2016). [CrossRef]
24. Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, “Image quality assessment: from error visibility to structural similarity,” IEEE Trans. on Image Process. 13(4), 600–612 (2004). [CrossRef]
25. E. Ng, F. Vasefi, B. Kaminska, G. H. Chapman, and J. J. Carson, “Contrast and resolution analysis of iterative angular domain optical projection tomography,” Opt. Express 18(19), 19444–19455 (2010). [CrossRef]
