FOURIER DOMAIN ITERATIVE APPROACH TO OPTICAL SECTIONING OF 3 D TRANSLUCENT OBJECTS FOR OPHTHALMOLOGY PURPOSES

In this paper we deal with the problem of optical sectioning. This is a post processing step while investigating of 3D translucent medical objects based on rapid refocusing of the imaging system by the adaptive optics technique. Each image, captured in focal plane, can be represented as the sum of in-focus true section and out-of-focus images of the neighboring sections of the depth that are undesirable in the subsequent reconstruction of 3D object. The problem of optical sectioning under consideration is to elaborate a robust approach capable of obtaining a stack of cross section images purified from such distortions. For a typical sectioning statement arising in ophthalmology we propose a local iterative method in Fourier spectral plane. Compared to the non-local constant parameter selection for the whole spectral domain, the method demonstrates both improved sectioning results and a good level of scalability when implemented on multi-core CPUs.


INTRODUCTION
Recently, a promising tool for investigation of 3D translucent medical objects was suggested which is based on the rapid refocusing of the incoherent imaging system by adaptive optics technique (see (Larichev et al., 2002)).By properly varying of the focal length it provides a stack of images each corresponding to different depths of cross section.In principle, having a sufficient amount of sections, one can try to reconstruct the original 3D object by known methods.However, a special step called sectioning is necessary for image processing of the stack obtained.Actually, images, captured in each focal plane, represent both infocus image of true section and a sum of out-of-focus images of the neighboring sections of the depth.Various aberrations of an optical system, fixation fluctuations, distortions of light-sensitive sensors, etc. should be taken into account as well.Thus, the problem is to elaborate a robust approach capable of obtaining a stack of cross section images that are purified from such distortions.
In this paper, we elaborate a 3D model of the incoherent imaging system which takes into account defocus between sections depending on their depth and high order phase aberrations.The problem consists in solving a 3D convolution equation.Direct deconvolution is an ill-posed problem.Analysis of contemporary 3D deconvolution methods (Wu et al., 2008) shows that iterative methods, compared with direct ones, are more accurate, but usually are accompanied by an increase in the complexity of computing and additional requirements for tuning of parameters.
The approach we propose utilizes the iterative algorithm in transversal Fourier domain, allowing parallelization of computations by modern multi core CPU/GPU software and hardware, in such a way that an increase in the complexity is not a matter of principle.In (Razgulin et al., 2015) the particular case of the approach was presented with a constant selection of parameters for all spectral components.It was found that for different values of the parameters the quality of recovery for high and low Fourier harmonics is different and can not be equally good for both cases.In this paper we utilize the effective local Fourier domain approach, based on the localization of the regularization parameter selection method.The results of computer simulations demonstrate the efficiency of the approach.The performance obtained by parallelization allows sectioning with transversal resolution up to 4K and acceptable time which is promising in noninvasive restoration of the 3D structure of human eye fundus in vivo for a consequent diagnosis of diseases.

3D MODEL OF IMAGING SYSTEM
We consider the following principal scheme of imaging system for optical sectioning (see Figure 1).A thin translucent 3D object The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLII-2/W4, 2017 2nd International ISPRS Workshop on PSBB, 15-17 May 2017, Moscow, Russia at a distance d1 behind the lens.Here, for the sake of simplicity, we suppose the parameters d0, d1 to be fixed and the focal length of the lens f (z) being a function of the variable z: By altering z from 0 to w one can obtain, at the image plane, a stack of translucent object cross-section images.Each crosssection image represents the true cross-section superposed with blurred neighboring sections.Further we discuss this fact in detail.
To this end, we consider infinitely thin translucent object illuminated by incoherent uniform plane wave and characterized by amplitude transmittance (optical density) OT (x , y , z ).If we place it at the plane z , 0 ≤ z ≤ w, and focus imaging system at the plane z with focal length f (z) according to (1), then, at the image plane, we observe the following distribution where is the image of the object in geometrical optics approximation, M (z) = d1 • (d0 + w − z) −1 -magnification.Point spread function h(ξ, η, z, z ) takes into account incoherent illumination and defocus caused by the misalignment of focal and object plane locations in the case z = z under the Fresnel scalar diffraction theory approach (Easton Jr., 2010, Gaskill, 1978): Strength of defocus is conducted by the function When w d0, d1 there is a spatially invariant approximation (4) In the case of thick N -layered translucent object composed of ∆z = w/N width layers placed at z n = n∆z, n = 0, 1, . . ., N with plane optical density OT (x , y , z n ) the amplitude transmittance of n-th layer is OT (x , y , z n )∆z .In the image plane, as mentioned in (Gu, 2000, p. 64), the total impact of all layers may be considered under linear superposition principle, which is valid for translucent thick objects where secondary diffraction is neglected.Thus the resulted distribution appears as a sum of corresponding images of each layer: I(x, y, z, zn)∆z.
When N → +∞ we obtain integral form Note that, under approximation (4), this formula can be easily represented in the form of 3D convolution integral, but this is not essential for us.From the relation ( 5) one can easily see that, due to the integration along z axis, the image of z-section of 3D object contains true cross-section and defocused images of neighboring sections of the depth.These distortions are undesirable in the problems of 3D object reconstruction from captured crosssections arisen in ophthalmology applications, bio-microscopy, etc. (Castleman, 2007, Wu et al., 2008).The problem is to elaborate robust approach capable to "purify" stack of captured crosssection images from such distortions.This procedure is further called "optical sectioning".

DISCRETE APPROXIMATION IN FOURIER DOMAIN
Equation ( 5) is a 3D integral Fredholm equation of the 1st kind.The problem of its solving is an ill-posed unstable problem.In addition, from the optical point of view, this is also related to the fact that corresponding optical transfer function (OTF), a Fourier transform of PSF, is close to zero in some regions of the Fourier frequency domain.A review of the advantages and disadvantages of the available direct and iterative methods for solving similar inverse problems of three-dimensional biomicroscopy, which differ from our statement in the form of a PSF function, is given, for example, in (Castleman, 2007, ch. 22.2) and (Wu et al., 2008, ch. 14).An analysis of these methods shows that the iteration methods are more accurate in comparison with the direct ones, but are accompanied by an increase in the computational complexity.However, with the proper choice of the iterative algorithm, which allows the parallelization of computations with modern software and hardware, such an increase in computing resources is not essential.
Typically, the sections of a three-dimensional image of an ophthalmic object are available for analysis and processing only in a finite set of cross-section planes z ∈ {z1, z2, . . ., zN }, equally located at a distance ∆z, where the number N is in orders of magnitude smaller than the number of points in each horizontal section.As we saw above, the same situation appears for a thick N -layered translucent object composed of ∆z-width layers.Thus it is reasonable to turn from equation ( 5) to its discrete analogue in depth variable z: where * is a sign of 2D convolution, Im(x, y) = I(x, y, zm), The Following to (Razgulin et al., 2015), we propose the iteration approach allowing effective parallelization.To start with, we use compactness of support property both for image Im(x, y) and PSF function h(x , y , zm, zn), which allows effective sampling of these functions and transforming continuous convolution integral to discrete convolution.Using discrete Fourier transform (DFT) and discrete convolution theorem we write equations ( 6) in Fourier spectral domain (u, v) using notations Im, Om, hmn for discrete Fourier coefficients of corresponding functions Im, Om, hmn: At each point (u, v) equations ( 7 Direct methods are unapplicable for solving (8) since H is a degenerate or ill-conditioned matrix in most interesting cases.That is why we use iterated Tikhonov regularization method (Bakushinsky and Goncharsky, 1989) where k = 0, 1, . . ., K, O 0 = (0, 0, . . ., 0).As it is mentioned in (Razgulin et al., 2015), in case of distortion present in observed data I, the parameter µ > 0 and the number of iterations K should be chosen according to the trade-off between a rapid draft sectioning with the smallest number of iterations for a rough visualization of layers at relatively large values µ and thorough restoration with fine identifying the texture due to additional iterations applied for small µ.

Phantom three-dimensional object
We have considered a phantom translucent three-dimensional object consisting of several vessels located on different layers.The following system parameters specific to the living three-dimensional structure of the eye were used: the depth of the object along z coordinate was 800 microns, the number of layers N = 20.Also a small defocus was introduced to each layer as a kind of stationary aberration.The specific form of three-dimensional object and corresponding true cross-sections are shown in Figure 2.
Since Poisson noise is typical for light-sensitive sensors used in imaging system at the image plane, sectioning was performed for the case of 2% Poisson noise was added to captured images, which corresponds to 100, 000 photons per unit of pixel intensity.
A typical view of noisy images for layers no. 3, 10, 19 is shown in Figure 3.

Restoration results for a constant µ
We utilize vector frequency criterion (VFC) to quantify the quality of image restoration in each cross-section focusing plane.VFC is defined by the function where Õ -recovered image, shifted to nonnegative area, Ooriginal object, D(X) = X/D(X), Ē(X) = X/E(X), D(•)dispersion, E(•) -average, F -Fourier transform, m(r) -number of points (u, v), satisfying the inequality r 2 ≤ u 2 + v 2 ≤ (r + 1) 2 .VFC allows one to control the restoration features in the spectral domain.
Figure 4 presents specific restoration results and graphs of the corresponding values of VFC for the section with the number of 3 for different values of the parameter µ for a fixed number of iterations K = 40 of the method (9).
As demonstrated in Figure 4, the value µ = 0.1(a) is too large for noisy images due to the high-frequency noise added to the restored image.Decrease of the parameter down to the value 0.01 improves the quality of restoration insignificantly.The VFC graph still shows degradation of the restored image in the highfrequency area.Further reduction of the parameter µ to the value 0.001 makes it possible to significantly reduce the noise level, but leads to the sectioning of the layers being insufficient.This can also be seen on the VFC graph, where a dip in the low frequency area is noticeable, although the high-frequency graph demonstrates a sufficient quality of restoration.

Restoration results for the parameter µ depending on the harmonic number
A significant improvement in the quality of the restored image can be achieved through choosing the parameter µ depending on the number of the harmonic to be restored.Since in the case of observable images without noise the sufficiently large values of the parameter µ are optimal, an analysis of the restored images was carried out for cases in which the low-frequency restoration was fulfilled with a larger value of the parameter µ, but with the same number of iterations.
Figure 5(a) shows the results of the reconstruction for the case of K = 40 iterations, when µ was equal to 0.1 for harmonics and equal to 0.001 for the remaining harmonics.Despite good sectioning, the restored image is noised by the low-frequency noise, which is confirmed by a jump in the corresponding frequency region of the VFC graph.By reducing the calculation area for low frequencies to (u, v) : √ u 2 + v 2 ≤ R/6, a significant improvement in the quality of the reconstructed image can be achieved.As we see from Figure 5(b), the reconstructed image is well sectioned from adjacent layers, and the VFC graph is close to unity throughout the region, including the high-frequency region.The above results show that, by fine-tuning of the parameter µ depending on the choice of the frequency domain, it is possible to achieve a significant improvement in the quality of the reconstruction of noisy images.If, in addition, a-priori information about the aberrations of spectral components of captured images is available then suboptimal filtration method (Sizikov, 1999) is applicable.

SOFTWARE IMPLEMENTATION
The developed iterative method of three-dimensional sectioning was implemented in C++ using the FFTW DFT library, the offi-cial release of BLAS (Basic Linear Algebra Subprograms) from Netlib for Windows within the LAPACK library.
OpenMP standard in C++ was used to parallel the implementation of the iterative method for several points on selected number of processor cores.Maximum size of the RAM used by each parallel process of the iteration method for some point (u, v) was 128 × (2N 2 + 3N ) bytes, which allowed to increase the maximum possible size of the processed data.
To measure the performance of the implemented algorithm, calculations for tasks with dimensions of 256 × 256, 512 × 512, 1024×1024 and 2048×2048 for 20 recoverable layers were performed on a personal computer with an Intel Core i7-4790K processor having 4 processor cores supporting Intel Hyper-Threading Technology, 16 GB of RAM under control of 64-bit operating system Windows 8.1.
With typical image sizes, as can be seen in Figure 7, the running time of the program changes little with the number of iterations and is acceptable for clinical practice when restoring medical image stacks.Note that one of the important features of the implemented algorithm is its ability of effective parallelizing using the multi-core architecture of modern CPUs.Furthermore, it allows to transfer a part of computations to the GPU, that, taking into account the specifics of modern multi-core GPUs, will potentially lead to a significant reduction of runtime for both DFT and SLAE.

CONCLUSION
A local iterative method in Fourier spectral plane is considered for a typical statement of three-dimensional sectioning problem arising in ophthalmology.Compared to non-local parameter selection, the method demonstrates both improved sectioning results and a good level of scalability when implemented on multicore CPUs.Additional quality improvements of sectioning can be achieved while utilizing grid warping techniques that provides both denoising and sharpening of images (Krylov et al., 2016).

Figure 2 .
Figure 2. The image of the 3D object under consideration and of true layers no. 3, 10.

Figure 6 Figure 6 .
Figure 6 compares the image restoration quality for cases the parameter µ (a), (b) constant to all harmonics numbers and case (c) described above on Figure 5(b).The recovery result for µ, depending on the harmonic number, is characterized by VFC graph closed to unity, in contrast to the VFC graphs for constant µ, which have jumps and irregular changes in the mid and high frequencies.Moreover, for the case under consideration, it was possible to achieve a good sectioning of the renewable layers from adjacent layers in combination with a sufficiently low level of input noise.

Figure 8 .
Figure 8. Performance acceleration graph for different image sizes depending on the number of threads.