New Approach for Limited-Angle Problems in Electron Microscope Based on Compressed Sensing

New advances within the recently rediscovered field of Compressed Sensing (CS) have opened for a great variety of new possibilities in the field of image reconstruction and more specifically in medical image reconstruction. In this work, a new approach using a CS-based algorithm is proposed and used in order to solve limited-angle problems (LAPs), like the ones that typically occur in computed tomography or electron microscope. This approach is based on a variant of the Robbins-Monro stochastic approximation procedure, developed by Egaziarian, using regularization by a spatially adaptive filter. This proposal consists on filling the gaps of missing or unobserved data with random noise and enabling a spatially adaptive denoising filter to regularize the data and reveal the underlying topology. This method was tested on different 3D transmission electron microscope datasets that presented different missing data artifacts (e.g, wedge or cone shape). The test results show a great potential for solving LAPs using the proposed technique.


Introduction
A common problem in Transmission Electron Microscope (TEM) is the limited-angle view, where large gaps of missing data are present in the acquired volume.An LAP occurs when having missing or corrupted data from a certain angle, which makes the acquisition unreliable and therefore you cannot trust the results.This group of LAPs is due to various technical and fundamental limitations on the minimum and maximum attainable tilt angles of the instrumentation that are used to acquire the data.Therefore, the data is confined to a limited-angularrange acquisition which prevents its appropriate visualization and rendering.
The goal of our work is to estimate the missing or corrupted data from TEM datasets by applying a new strategy instead of the classically used techniques, such as POCS (Projection Onto Convex Sets), [1][2][3], which is widely used due to its easy implementation but offers poor results, or the more sophisticated ones like PICCS (Prior Image Constrained Compressed Sensing) or TV (Total Variation) [4], which are more complicated solutions.
In this paper, we study a new variant of the already mentioned method of Egiazarian, [5].We present and explain how the technique has been adapted and applied, and present the obtained results and a brief discussion along with the conclusions extracted from the results.

Method
In [6][7][8], it is shown that under CS assumptions, stable reconstruction of partly unknown data is possible and that in some cases the reconstruction can be exact.These techniques typically rely on convex optimization with a penalty expressed by the l 0 or l 1 norm, which is exploited to enable the assumed sparsity [9].Therefore, it results in parametric modelling of the solution and in problems that are commonly solved by mathematical algorithms.In [5], it is proposed to replace the traditional parametric modelling used in CS by a non-parametric one.This non-parametric modelling is implemented by the use of a spatially adaptive filter.The regularization imposed by the l 0 or l 1 norm is essentially only a tool to design some non-linear filtering.This implicit regularization is replaced by explicit filtering, exploiting a spatially adaptive filter, which is sensitive to the features and details of the image.If this filter is properly designed, there is a chance to achieve better results than those obtained by the common method based on formulation of imaging.In imaging, the regularization with global sparsity penalties (e.g, l p norms in some domain) often results in an inefficient filtering.It is known that a higher quality result can be achieved when the regularization criteria are local and adaptive.
The applied method to reconstruct the "dead zones" (e.g. with corrupted or missing data) is carried out by a recursive algorithm based on spatially adaptive denoising filtering [5].Each iteration consists of providing data to a block-matching and a spatially adaptive 3D filtering algorithm, called BM3D, by the injection of random noise in the unobserved portion of the data in frequency domain.To carry out block-matching and block-filtering, the filter implemented in [10] was used with some modifications such as the removal of the Wiener filter in order to speed up the algorithm and using Haar wavelet for the third-dimension transform.In addition, some parameters of the block hard-thresholding (HT) were optimized to speed up the computations, such as N 1 = 4, N 2 = 12, Ns = 31 and tau_match = 1800.This peculiar filter works in the image domain.It attenuates the noise and reveals new details and features of the corrupted image at each iteration.The process ends when a specific number of iterations are performed.
The algorithm is ruled by the following recursive system: where: Φ ≡ filtering block Image data (denoted as y) is divided into a known portion (denoted as y 1 ) and an unknown portion (denoted as y 2 ) as follows: ( ) . ? y S y S y y y The procedure of [5] was applied to reconstruct 2D images (e.g. the data in Figure 1(b)) and achieve a perfect reconstruction, while we apply the method to 3D TEM volume images (e.g. the data in Figure 2) by splitting the 3D problem into a number of different 2D problems and solve them separately using a different configuration for each case due to the variability of the dead zone from slice to slice.For example, for the case presented in Figure 2, when taking horisontal slices in the xy-plane in Fourier domain.

Datasets
Three different datasets were used in this work.

Hansandrey Crystallography Dataset
It consists of a 3D artificial crystallography of size 100 × 100 × 100 voxels in.mrc format, [11].It was created by a simulator-software to evaluate the performance of TEM reconstruction algorithms [12].In this case, we simulated a missing wedge and a missing cone in Fourier domain (as shown Figure 3) to apply the reconstruction procedure and evaluate the quality of the result.

Viral DNA Gatekeeper Dataset
It consists of a 3D model of a viral DNA gatekeeper, of size 100 × 100 × 100 voxels in.mrc format, obtained by cryoelectron microscopy technique.For this case, we only simulated a missing cone in the data to be able to apply the reconstruction method and evaluate its performance of filling up the missing cone.

Philip's Crystallography Dataset
It is a 3D model of a protein molecule acquired by a TEM but suffering from an LAP since it was not possible to tilt the specimen more than 45 degrees.Therefore, the 3D model has (in Fourier domain) a missing/corrupted cone with a vertex angle of 90 degrees.The size of this 3D model is 80 × 80 × 80 voxels in.mrc format.

Experimental Results
The peak signal-to-noise ratio (PSNR) is used as an evaluation parameter to compare the test results.The PSNR is defined as follows: ( ) Copyright © 2013 SciRes.ENG where max I is the maximum voxel value of the 3D model, I is the original 3D model, and Î is the reconstructed 3D model.
PSNR values of 47.6 dB, 36.5 dB and 46.2 dB were obtained for the 3D reconstruction of the Hansandrey dataset in the cases of missing wedge, missing cone when considering xy-plane slices (i.e.horizontal slices in Figure 3) and missing cone when considering xz-plane slices (i.e.radial slices in Figure 3), respectively.We can visually evaluate the quality of the best reconstructed result in Figure 4 (an open-source 3D visualization software called Chimera was used).Regarding the result of the case of missing cone when using xy-plane slices, we attribute this poor PSNR value to the large "dead zones" (i.e.empty discs) present in the slices at the edges of the 3D model.The algorithm has a limitation and fails in filling these large gaps of missing data.For the viral DNA gatekeeper dataset, the reconstruction ran as expected and a PSNR value of 45.9 dB was obtained.Visual comparison between the original 3D model and the reconstructed one, shown in Figure 5, assures that this high PSNR value corresponds to good visual similarity.The differences between the original 3D model and the reconstruction are almost negligible.
In the case of Philip's crystallography dataset, the reconstruction result and the original 3D acquisition are presented in Figure 6.In this case, there is no reference 3D model that can be used to verify whether the reconstruction result is satisfactory or not.
Visually speaking, we cannot detect large differences between the original model and the reconstructed one, but the difference in terms of a calculated PSNR value of 21 dB is quite considerable.But this doesn't mean that the method doesn't work properly.It reflects that new details appear in the reconstructed volume which could be closer to the real 3D model without missing data.

Discussion and Conclusions
The purpose of our work was to apply a new method to solve LPAs (e.g. with a missing cone or a missing wedge of the acquired image-data volume) in TEM data.The method was applied to different datasets and the resulting reconstructed image volumes were evaluated.Good 3D reconstruction results were obtained.The PSNR values were calculated for all resulting reconstructed images.These PSNR values were better than those obtained using existing commonly used techniques, such as POCS.
The approach proposed in our work can be considered as a great breakthrough, because for data acquisitions limited to [45˚, −45˚], POCS results in an error-rate around 40%, while our approach achieves an error-rate lower than 1% for the Hansandrey and the viral DNA gatekeeper cases when 50% of the acquired data is missing.However, different acquisition technique and procedures will produce data with different sparsity characteristics in frequency domain, which in its turn will affect the performance of our method.For example, if a large portion of the high frequency zones of the acquired data is missing or corrupted, then it gets much more difficult to reconstruct the missing part of the 3D model because the algorithm doesn't have enough prior information.
Therefore, we have to take under consideration that a test measure is needed to determine which datasets, with the presence of missing data, are valid or not to apply the proposed method and get good 3D reconstruction results.If such a test is performed, it will be possible to know if the obtained 3D reconstruction result is supposed to be similar to the real original model (i.e.without missing data) or not.Then it will be possible to know if the 3D reconstruction of Philip's crystallography dataset (presented in Figure 6) is correct or not.
One of the most exciting research projects that could emerge from our work is the possibility to develop a new optimized acquisition technique or procedure for TEM.In addition, achieving a considerable reduction of radiation dose applied to the specimen.Another possibility is to adapt the proposed method and apply it to other kinds of modalities like Computed Tomography (CT), Magnetic Resonance Imaging (MRI), astronomy, geophysical exploration or other type of electron microscope techniques.Since a complete reconstruction of the Hansandrey model took 32 hours in a common laptop (4-core 2.0 GHz and 4 Gb RAM), it would be necessary to speed up the algorithm by implementing it using GPU techniques (e.g.CUDA, OpenCL).

2 ŷ
≡ estimation of unknown data in Fourier domain y 1 ≡ known data in Fourier domain γ ≡ speed step of the algorithm (1 − S) ≡ mask to select the region of the unknown data ϒ ≡ Fast Fourier Transform 1 − ϒ ≡ Inverse Fast Fourier Transform k η ≡ Gaussian noise

Figure 2 .
Figure 2. 3D limited-angle example in the image domain.