CC BY
27
Инженерия в медико-биологической, образовательной и спортивной практике
SPACE-VARYING RESTORATION OF DIFFUSE OPTICAL TOMOGRAMS RECONSTRUCTED BY THE FILTERED BACKPROJECTION ALGORITHM
A.B. Konovalov *, D.V. Mogilenskikh *, V.V. Lyubimov **
* Russian Federal Nuclear Center- Institute of Technical Physics, Snezhinsk, Chelyabinsk Region, 456770 Russia,
** Research Institute for Laser Physics, St-Petersburg, 199034 Russia
Possibility is investigated to enhance spatial resolution of diffuse optical tomograms reconstructed by the photon average trajectories (PAT) method. The PAT method is based on a concept of average statistical trajectory of light energy transfer from point source to point detector. The inverse problem of diffuse optical tomography reduces to solution of an integral equation with integration by conventional PAT. In the result for reconstruction of diffuse optical images the conventional algorithms of projection tomography can be applied, including filtered backprojection algorithm. The shortcoming of the PAT method is that it reconstructs I mages blurred in the result of averaging by photons spatial distribution contributing into the signal measured by a detector. To enhance resolution we apply a spatially variant blur model based on interpolation of spatially invariant point spread functions simulated for different image regions. To restore tomograms two iterative algorithms for solution of system of linear algebraic equations are used: conjugate gradient algorithm for least squares problems and modified residual norm steepest descent algorithm. It is shown that one can achieve 30% enhancement of spatial resolution.
1. Introduction
To reconstruct the diffuse optical tomograms with best spatial resolution «multi-step» techniques such as Newton-like and gradient-like ones [1-3], which use exact forward models, are generally applied. These techniques require a few minutes for reconstruction and are inapplicable for real-time medical exploration. Recently, new method [4-6] based on a probabilistic interpretation of the photon energy transport from a source to a detector has been proposed by us for reconstruction of optical macroinhomogeneities. By this method, the inverse problem of diffuse optical tomography (DOT) is reduced to a solution of an integral equation with integration along an average statistical trajectory for transfer of light energy, photon average trajectory (PAT), which is curvilinear in the common case. As a result the PAT method may be implemented as a «one-step» technique with the use of the algebraic reconstruction algorithms [5]. Moreover within the internal zone of the object, well away from the boundaries, PATs tend to a strait line and the linear integral transforms such as the Radon inversion can be successfully used for real-time reconstruction of the images [6].
The main problem is that the PAT method recon-
structs the images blurred due to averaging over spatial distributions of photons, which form the signal measured by the detector. Two ways are available to improve the spatial resolution of the PAT tomograms. The first approach is to apply the backprojection algorithms with special filtration of shadows (Vaingerg or hybrid Vainberg-Butterworth filtration [6, 7]). This approach allows a 20%-gain in spatial resolution to be obtained but does not restore the inhomogeneity profile. In the second way, the blurred images are restored with the use of an available blur model. As the resolution is spatially variant due to multiple light scattering, the blur model can not be described by a convolution and the linear inverse filtering as well as nonlinear deconvolution algorithms can not be used for deblurring the PAT tomograms in a strict sense. In [8] we have assumed the convolution blur model to restore not a full image, but its local internal regions containing the optical inhomogeneities, wherein the point spread function (PSF) is weakly space-variant. This approach is not handy for restoring an image with a few inhomogeneities, as each local region containing the inhomogeneity uses its own spatially invariant PSF, and the restoration results must then sew together to obtain the full restored image. Moreover we
have observed the considerable profile distortions for the case of big size macroinhomogeneity. In present paper we use the spatially variant blur model developed by J. Nagy with his colleagues [9]. According to it the blurred image is partitioned into subregions. However, rather than deblurring the individual subregions locally and then sewing the individual results together, this method interpolates the individual PSFs, and restores the image globally. To study the efficiency of such approach, the numerical experiment is conducted, wherein cross-sections of cylindrical strongly scattering objects with absorbing inhomogeneities are reconstructed using the PAT method with the standard backprojecton algorithm, the spatially invariant PSFs for all individual subregions are simulated, and two well-known algorithms for solving a system of linear algebraic equations are applied to calculate the original image. Those are the conjugate gradient algorithm for least squares problems (CGLS) [10] and modified residual norm steepest descent algorithm (MRNSD) [11,12].
2. Theoretical backgrounds
Let us the photons migrate in strongly scattering media from a source space-time point (0,0) to a receiver space-time point (r,t) . A relative contribution of photons located at an intermediate space-time point (rj. r) to the value of photon density at (r,t) can be characterized by a conditional probability density:
P(r,t)
(1)
where P(r,t) is a probability density of the photon migration from (0,0) to (r,t). If the photon density <p(r,t) satisfy the time-dependent diffusion equation for a volume V with a limited piecewise-closed smooth surface for an instantaneous point source and the Robin boundary condition, the probability density P(r\,T\r,t) are expressed as
P(rbr;r,0 =
<p(ri,T)G(r-rut-T) jV(ri,т)G(r -r\,t-T)d3i\
(2)
where G(r,t) is the Green function. Let us define a relative shadow g as a logarithm of the relation between the value of the signal intensity I caused by presence of the inhomogeneities and the value of unperturbed signal intensity Iq , measured at the object
surface at the time moment t. The fundamental equation of the PAT method in the case of the time-domain measurement technique are written as
g(A0= \S{rx,r)P(rvr,r,t)d3rx
L nV\4 Vk
dh (3)
where c is a light velocity in vacuum, n is a refraction index of a media, I is a PAT from a source to a receiver, v(l) is a velocity of the mass centre of photon distribution, which moves along the PAT, S(r,t) is a inhomogeneity distribution function. In the case of an absorbing inhomogeneity S(r,t) = r), where
8fxa{ r) is a local disturbance of the absorption coefficient fiair) • Using the backprojection algorithm, equation (3) may be directly inverted in relation to the function
(S(r,t)) = |S’(r1,T)/>(r1jr;r,0A , (4)
v
viz. the function blurred due to averaging over the spatial distribution of photons, which form the signal measured by the receiver at the moment t. Let f be a vector representing the original image of an absorbing inhomogeneity Sjua(r) and ^f) be a vector representing the image blurred due to averaging (4). The spatially variant blur model are written in a matrix form
as
f =A-f
У
A = ^rD/A; ,
/=1
(5)
(6)
where A is a large ill-conditioned matrix that models the blurring operator. If the image is partitioned into p subregions, the matrix A has the following structure where A, are the block Toeplitz matrices with Toeplitz blocks [13] and D, are diagonal matrices satisfying ^ D, = I,
Piecewise constant interpolation implemented implies that the k th diagonal entry of D, is one if the k th point is in subregion i, and zero otherwise.
3. Results
To demonstrate the effect of improving the spatial resolution of diffuse optical images reconstructed by the PAT method, a numerical experiment was conducted, wherein cross-sections of cylindrical strongly scattering objects were reconstructed. The diameter of objects was equal to 6.8 cm. The refraction index, the reduced scattering and absorption coefficients of the objects were equal to 1.4, 5.0 and 0.05 cm'1, correspondingly. We considered two sets of phantoms. Each phantom of the first set contained a cylindrical absorbing inhomogeneity with the diameter equal to 1.4 cm (absorption coefficient was equal to 0.075 cm-1). In one of the cases, the inhomogeneity was located in the center of the object, in the two others it was displaced from the center by 1.4 and 2.5 cm, correspondingly. The second set of phantoms was designated to define the modulation transfer function (MTF) that is commonly used for spatial resolution estimation. We used
Konovalov A.В., Mogilenskikh D.V., Lyubimov V.V.
Space-varying restoration of diffuse optical tomograms reconstructed by the filtered back projection algorithm
four cylindrical strongly scattering objects, each containing two cylindrical absorbing inhomogeneities equal in diameter. The diameter and optical pa-
x 1(T3
gain in spatial resolution with improvement of image profiles can be obtained due to the space-varying restoration.
x 10~2 x 1СГ2
(c)
The example of numerical experiment results: (a) the blur tomogram, (b) the restored tomogram, (c) the image profiles:
1 - ideal, 2 - blurred, 3 - restored. The points in the tomograms represent the positions of the sources on the object boundary. The coordinate axes are graduated in centimeters and the intensity scale in reverse centimeters
rameters of these objects, as well as the absorption coefficient of inhomogeneities, were identical to those for the phantoms of the first set. The distance between inhomogeneities was equal to their diameter. Diameters of inhomogeneities of different objects were equal to 1.4, 1.2, 1.0, and 0.8 cm. Sources (32) and receivers (32) were installed along the perimeter of the objects at equal step angles (11.25°). In so doing, the angular distance between the nearest-neighbor source and receiver constituted 5.625°. The relative shadows caused by the absorbing inhomogeneity were simulated via numerical solution of time-dependent diffusion equation for the instantaneous point source (the case of time-domain measurement technique) with the use of the FEM method. Reconstruction of cross-sections for each phantom with the use of the backprojection algorithm was realized onto rectangular grid 63*63.
To create the matrix A the PSFs for all subregions were simulated by the FEM method (shadows) and the backprojection algorithm (reconstructions). The CGLS and MRNSD algorithms were applied for image restorations. The restoration result obtained by the MRNSD algorithm for the case of one inhomogeneity displaced from the center by 1.4 cm is presented in Figure as an example. The image at that was partitioned into 16 (4x4) subregions.
The results of the MTF estimations show that a 30%-gain in spatial resolution can be obtained due to the space-varying restoration.
4. Conclusion
The paper investigates possibility to apply iterative algorithms of image restoration to improve spatial resolution of diffuse optical tomograms reconstructed by the photon average trajectories method. To describe a blurring operator we used a spatially variant blur model. The obtained results testify the undoubted efficiency for proposed approach. A 30%-
5. Acknowledgments
The authors would like to thank professor J. G. Nagy and his colleagues at Emory University for their code package «Restore Tools» provided for calculations. This work was supported in part by the International Science and Technology Center under Grants 2151 and 2184.
References
1. Arridge S.R., «Optical tomography in medical imaging», Inverse Probl. 15, R41-R93 (1999).
2. Yodh A. and Chance B., «Spectroscopy and imaging with diffusing light», Phys. Today 48, 34-40
(1995).
3. A.H. Hielscher, A.D. Klose and KM. Hanson, Gradient-based iterative image reconstruction scheme for timeresolved optical tomography, IEEE Trans. Med. Imaging 18, 262-271 (1999).
4. V. V. Lyubimov, «Optical tomography of highly scattering media using first transmitted photons of ultrashort pulses», Opt. Spectrosc. 80, 616-619
(1996).
5. V.V. Lyubimov, A.G. Kalintsev, A.B. Konovalov, O.V. Lyamtsev, O.V. Kravtsenyuk, A.G. Murzin, O. V. Golubkina, G.B. Mordvinov, L.N. Soms and L.M. Yavorskaya, «Application of photon average trajectories method to real-time reconstruction of tissue inhomogeneities in diffuse optical tomography of strongly scattering media», Phys. Med. Biol. 47, 2109-2128 (2002).
6. A.B. Konovalov, V. V. Lyubimov, 1.1. Kutuzov, O.V. Kravtsenyuk, A.G. Murzin, G.B. Mordvinov, L. N. Soms and L.M. Yavorskaya, «Application of the transform algorithms to high-resolution image reconstruction in optical diffusion tomography of strongly scattering media», J. Electron. Imaging 12, 602-612 (2003).
7. V.V Lyubimov, O. V. Kravtsenyuk, A. G. Kalintsev, A.G. Murzin, L.N. Soms, A.B. Konovalov, 1.1. Ku-
tuzov, O.V. Golubkina and L.M. Yavorskaya, «The possibility of increasing the spatial resolution in diffusion optical tomography», J. Opt. Technol. 70, 715— 720 (2003).
8. A.B. Konovalov and V. V. Lyubimov, «High-resolution restoration of diffuse optical images reconstructed by the photon average trajectories method» in Optical Technologies in Biophysics and Medicine V, V.V. Tuchin, ed., Proc. SPIE 5474, 66-79 (2004).
9. K.P. Lee, J. G. Nagy and L. Perrone, «Iterative methods for image restoration: a Matlab object oriented approach», see http:/ /-www.mathcs.emory.edu/ ~nagy/RestoreTools/.
10. A. Bjdrck, Numerical Methods for Least Squares Problems (SIAM, Philadelphia, PA, 1996).
11. L. Kaufman, «Maximum likelihood, least squares, and penalized least squares for PET», IEEE Trans. Med. Imag. 12, 200-214 (1993).
12. J. G. Nagy and Z. Strakos, «Enforcing nonnegativity in image reconstruction algorithms», in Mathematical Modeling, Estimating, and Imaging, D. C. Wilson et. al., eds., Proc. SPIE 4121, 182-190 (2000).
13. M. Hanke and J.G. Nagy, «Restoration of atmospherically blurred images by symmetric indefinite conjugate gradient techniques» Inverse Probl. 12, 157-173 (1996).
