Linear Reconstruction Methods for Large Thick Aperture Imaging

Large thick aperture imaging method is proposed to measure the radiation intensity distribution of radiation source whose size is several centimetres. The new method contains two steps which are coded imaging and image reconstruction. Geant4 is used to simulate the coded imaging process. Wiener filter, Lucy-Richardson, Regularized filter, blind deconvolution and CT reconstruction method are realized using Matlab software. Wiener filter and Regularized filter produce artifacts. CT reconstruction method can’t reveal the right size of 60Co-γ source because of the overlapping of knife edge images. Both L-R and blind deconvolution are effective. If the aperture is moved 0 mm, 2 mm, 4 mm, 6 mm along x-axis, the reconstruction result using L-R method has a larger correlation coefficient with source image. However, LR method gets bad result if the aperture rotates more than 1.4°. Under this condition, blind deconvolution is the better choice.


Introduction
Radiography technology of high energy γ ray or neutron can reveal the shape and size information of the radiation source. Since 1969, Los Alamos starts developing pin hole imaging [1]. The typical size of pinhole is 0.5 mm [2]. Image acquired through small size pinhole has high spatial resolution. But the radiation collection efficiency is low that makes it difficult to get an image of low intensity radiation source. Compared with pinhole, penumbral aperture has higher radiation collection efficiency. Penumbral imaging [3][4][5] has been successfully used in micrometre size ICF target plasmas diagnosis. However, the penumbral aperture size must be larger than the radiation source. For centimetre size [6] radiation source diagnosis, penumbral imaging faces some problems such as large area image detector design and radiation shield.
Large thick aperture imaging [7] is proposed to diagnose low intensity radiation of centimetre size. It can raise the radiation collection efficiency for about 1~3 orders comparing with traditional pinhole imaging [2]. Its aperture size can be smaller than the radiation source which results in source information involved in both central and penumbral area of the coded image. Coded image must be reconstructed to acquire source information. Although nonlinear reconstruction methods have been developed in recent years [4,8], they all have shortcomings such as long time cost, complicated and great amount of computing resources requirement. As a result, linear reconstruction methods are still widely used. This article simulates the large thick aperture imaging process of 60 Co-γ source and compares the reconstruction effectiveness of different linear reconstruction methods.

Coded imaging simulation
Monte Carlo method [9] can simulate the passage of particles through matter. Geant4 [10] is a toolkit developed by CERN and KEK in USA. It includes a complete range of functionality including tracking, geometry, physics models and hits. The physics processes offered cover the main interactions of γ-ray with different materials including photoelectric effect, Compton scattering, gamma conversion into an electronpositron pair. This paper adopts Geant4 to simulate the coded imaging process of 60 Co-γ source [11]. Main function calls several subroutines. Their function is introduced as below: (1) G4VUserDetectorConstruction: Define material and geometry of the aperture and detector.
(2) G4OpenInventor: Visualization of the geometry and particle track.
(4) G4VUserPrimaryGeneratorAction: Set parameters of emitting particles including type, position, energy and direction.
(6) G4VSteppingAction: Record particle information. Geometry of the program is shown in Fig 1.(mm) Image recorder Figure 1 Schematic of large thick aperture imaging Material of the aperture is tungsten. Shape of source is cylinder. Angle between γ emitting direction and zaxis is in the range of 0-1.3° so that simulation time is reduced significantly. Position information is recorded if γ-ray reaches the image recorder plane. If diameter of the source D, object distance s 1 , image distance s 2 and diameter of the aperture d satisfy Eq.(1), PSF(Point Spread Function) on image surface is larger than coded image. PSF is space invariant. In this case, linear reconstruction method can be used.

Principle and realization of linear reconstruction methods
There are four linear reconstruction methods in Matlab software. They are wiener filter, Lucy-Richardson, regularized filter and blind deconvolution. CT reconstruction method is also discussed because of its successful application in penumbral imaging. Generally, large aperture imaging process can be described by Eq. (2).
Taking the Fourier transform of both sides of Eq.(2) results in where O(u,v), and N(u,v) are the Fourier transform of coded image, source image and additive noise respectively.

Wiener filter
Wiener filter method searches for a reconstruction result ô (x,y) who has the least square error with o(x,y). Expression in frequency domain is shown as Eq.(4).
where S n (u,v)/S o (u,v) is the noise-to-signal power ratio of the additive noise(NSR). Wiener filter [12] can be realized using deconvwnr function in Matlab. In general, NSR is unknown and difficult to estimate accurately. So it is a variable input parameter in the function.

Lucy-Richardson
The algorithm is based on maximizing the likelihood of the reconstruction result's being an instance of the source image under Poisson statistics. Iteration equation [13] is shown as Eq.(5). Fig 1, L-R method is linear. Otherwise, it is nonlinear. L-R method can be realized using deconvlucy function in Matlab.

Regularized filter
The algorithm is a constrained optimum in the sense of least square error between ô (x,y)and o(x,y) under requirement of preserving image smoothness [12]. Expression in frequency domain is shown as Eq. (6).
where γ is a variable parameter, P(u,v) is the Fourier transform of the Laplacian. Regularized filter can be realized using deconvreg function in Matlab and NP (Noise Power) is a variable input parameter related to γ.

Blind deconvolution
In general, all of the reconstruction methods in which PSF is not a necessary input parameter are blind deconvolution. CT reconstruction method introduced below is one kind of blind deconvolution methods. The deconvblind function in Matlab deconvolves i(x,y) using the maximum likelihood algorithm like L-R method [12]. Both the o(x,y) and h(x,y) are returned.

CT reconstruction method
Knife edge imaging method [14] is widely used in the measuring of integral intensity distribution of radiation source in one direction. Large thick aperture can be seen as 360 knife edges. In this sense, intensity distribution of each angle in the coded image is a knife edge image. Differential of the knife edge images are the integral intensity distribution in 360 directions. If they are plotted in one photograph, we get sinogram. Then filter back projection method can be adopted to acquire the source information. There is no function in Matlab realizing CT reconstruction method directly. Diff function and iradon function are used to realize differential operation and filter back projection respectively [8].

Effectiveness evaluation
Indicators of the digital image quality evaluation include mean square deviation, peak signal to noise ratio, fidelity and correlation coefficient. Each indicator has limitation and can't characterize all the features of a digital image. Shape of the 60 Co-γ source in this article is cylinder. It turns out to be a circular uniform surface source observing along z-axis. Correlation coefficient is chosen as the indicator of the difference between ô (x,y) and o(x,y). Cor of image A and B is calculated using Eq.(7). , 2 2 , , x y x y x y The closer Cor is to 1, the smaller difference is between A and B in statistics [15]. Besides, direct perception is also important and need to be taken into consideration.

Ideal collimation
During the process of collimating the large thick aperture imaging system, centres of the source plane and image plane define a line which is reference optical axis.  Table 1 is maximized by adjusting the corresponding variable parameters in different reconstruction methods. It is easy to find that Wiener filter, L-R method, regularized filter and blind deconvolution can all restore the outline of the 60 Co-γ source. Comparing with L-R method and blind deconvolution, Cors of the reconstruction results using Wiener filter and regularized filter are closer to 1. However, obvious artifacts occur in (d) and (f) resulting in poor perception. Besides, NSR and NP are difficult to estimate precisely without a priori knowledge of the source. From Fig 2(h) we can see that CT reconstruction method gives out a source diameter 6 mm which is smaller than the actual diameter 20 mm. Reason for this is analysed as below. Raito of s 1 to s 2 in Fig 1 is 1:2. One knife edge image size is twice of the source size. That is Φ40 mm. Coded image size is about Φ80 mm in Fig 2(b). Coded image can be seen as 360 knife edge images fusing together. Since one knife edge image size is too large. They overlap when fusing together. Differential of one direction intensity distribution in coded image can't reflect the integral intensity distribution in source image, which leads to the ill-condition of CT reconstruction method. As a result, L-R method and blind deconvolution are preferred to be used in reconstructing the large thick aperture coded image.

Non-ideal collimation
The adjustment platform for putting the aperture on can perform lifting, translating, rotating and pitching operation. A laser theodolite is used to set up a reference optical axis which goes through both the source plane centre and image plane centre. Since the spot size of laser is not infinitesimal (~1 mm), errors may occur when adjusting the central axis of the aperture to coincide with the reference optical axis. Error is divided into two types: translation error and rotation error. Coded images are simulated under the condition of 2 mm, 4 mm, 6 mm translating distance and 0.95°, 1.91°, 2.86° rotating angle. Cors between the reconstruction results and the source image using L-R method and blind deconvolution are plotted in Fig 3 and      Reconstruction results with 2.86° rotating angle Reason for this is analysed as below: PSF in L-R method is a permanent input parameter. PSF has changed when the aperture rotates. But we still use the PSF shown in Fig 2(c) since we think the collimation process is ideal by observing. In the contrary, blind deconvolution will look for the real PSF automatically with the iteration process going on.

Conclusion
This paper discusses linear reconstruction methods for large aperture imaging. Principle of Wiener Filter, L-R method, Regularized Filter, blind deconvolution and CT reconstruction method is introduced. Realization method using Matlab software is also introduced. Result acquired by CT reconstruction method is ill because knife edge images overlap when fusing together. The other 4 methods can all restore the outline of the 60 Co-γ source. Wiener filter and regularized filter cause obvious artifacts resulting in poor perception. Under ideal collimation and translation error, both L-R method and blind deconvolution are effective. L-R method is better than blind deconvolution. However, if the aperture rotating angle is larger than 1.4°, blind deconvolution method should be chosen first because of its good faulttolerance ability for input parameter PSF. Experiment work will be done after finishing the aperture design in the future.