Novel Regularized Image Reconstruction Algorithms for Sparse Photoacoustic Tomography
Among all tissue imaging modalities, photo-acoustic tomography (PAT), has been getting increasing attention in the recent past due to the fact that it has high contrast, high penetrability, and has capability of retrieving high resolution. By using the combination of optical absorption and acoustic wave propagation, PAT has been able to image tissues at a relatively large depths with high resolution compared to purely optical modalities. Upon shining with a laser pulse, the substance under investigation absorbs optical energy and undergoes thermoelastic expansion; as a result, the spatial distribution of the concentration of the substance gets translated into the distribution of pressure-rise. This initial pressure rise travels outwards as ultrasound waves which are collected by ultrasound transducers placed at the boundary. From the ultrasound signal measured by the transducers as a function of time, a PAT reconstruction method recovers an estimate of the initial pressure-rise by solving the associated inverse problem. The inverse problem is however challenging. It is challenging because the image has to be recovered for the entire cross-sectional plane, whereas the samples of the acoustic pressure are available only from the points lying in the periphery of the imaging specimen where the transducers are located. In this thesis, we make contributions in two widely used types of reconstructions methods known as the time-reversal method, and the model based method. We summarize our contributions in three parts in the following. In the first part, we develop an improved model based method. Model-based reconstruction methods in PAT express the measured pressure samples as a linear transformation on the initial pressure-rise and perform a regularized reconstruction. Model based methods yield superior image quality even in the situation where measured data size is small. We propose a model-based image reconstruction method for PAT involving a novel form of regularization and demonstrate its ability to recover good quality images from datasets of significantly reduced size. The regularization is constructed to suit the physical structure of typical PAT images. We construct it by combining second-order derivatives and intensity into a non-convex form to exploit a structural property of PAT images that we observe: in PAT images, high intensities and high second-order derivatives are jointly sparse. This regularization is combined with a data fidelity cost, and the required image is obtained as the minimizer of this cost. As this regularization is non-convex, the efficiency of the minimization method is crucial in obtaining artefact-free reconstructions. We develop a custom minimization method for efficiently handling this non-convex minimization problem. Further, as non-convex minimization requires a large number of iterations and the PAT forward model in the data-fidelity term has to be applied in the iterations, we propose a computational structure for efficient implementation of the forward model with reduced memory requirements. We evaluate the proposed method on both simulated and real measured data sets and compare them with a recent reconstruction method that is based on well-known total variation regularization. Appropriate tuning of the regularization weight, lambda , plays a crucial role in determining the quality of reconstructed images in PAT. To make any regularization method practicable, we need to have a way to determine the lambda from the measured data. Unfortunately, an appropriately tuned value of the regularization weight varies significantly with the variation in the noise level, as well as, with the variation in the high resolution contents of the image, in a way that has not been well understood. In the part of the work described above, we did not address this problem as the focus has been to demonstrate the suitability of the intensity augmented regularization for PAT image recovery; in the experimental demonstration, we determined the required regularization weight by using the models that generated data. In the second part of the thesis, we develop a semi-automatic method for determining the regularization weight from measured data. As a first step, we introduce a relative smoothness constraint with a parameter; this parameter computationally maps into the actual regularization parameter, but, its tuning does not vary significantly with variation in the noise level, as well as with the variation in the high resolution contents of the image. Next, we construct an algorithm that integrates the task of determining this mapping along with obtaining the reconstruction. Finally we demonstrate experimentally that we can run this algorithm with a nominal value of the relative smoothness parameter—a value independent of the noise level and the structure of the underlying image—to obtain good quality reconstructions. We compare the structural similarity (SSIM) scores of reconstruction obtained this way to that of reconstructions in which the regularization weight was determined using the models themselves; we show that the SSIM scores are comparable. This means that, in a practical point of view, our work solves the problem of determining the required regularization weight from measured images. In the first two parts, we assumed that the forward model that measures the signal from target object to be ideal. In particular, we assumed that the excitation pulse and transducers impulse response are Dirac deltas. We focused only on the non-ideality of the transducer configuration, i.e., we handled the case where the transducer locations do not densely sample the detection surface as required by well-known back-projection method to work. In reality, both excitation pulse and transducer impulse response have finite width and this leads to some distortions in the reconstructed image. In the last part of the thesis, we propose a pre-processing method for correcting the distortions in the context of using time-reversal methods which are similar to back-projection method. To this end, we formulate the broadening of the PA signals as a convolution between the impulse response of the system and the input excitation pulse. A deconvolution method using Tikhonov regularization is proposed to correct the PA signals before applying the time-reversal method. This resulted in improved resolution in the reconstructed images. A two level deconvolution with Tikhonov regularization method is also proposed to remove the blurring caused by the finite bandwidth of transducers and by the broad excitation pulses. We evaluate the usefulness of our method using numerical simulations and demonstrate that the reconstructed images from the deconvolved PA signals remain unaffected by change in pulse widths or pulse shapes, as well as by the limited bandwidth of the ultrasound detectors.