CC BY
31
ITERATIVE RESTORATION OF X-RAY IMAGES TAKEN IN X-PINCH RAYS
A.B. Konovalov * V.V. Vlasov * AS. Uglov *, T.A. Shelkovenko **, S.A. Pikuz ** *Russian Federal Nuclear Center - Institute of Technical Physics, Snezhinsk, Chelyabinsk Region, 456770 Russia; ** Cornell University, Ithaca NY, 14853 USA
Imaging in X-pinch rays produces records of several superimposed components of an object image from several radiation focuses separated in space. The present paper studies the feasibility of restoring an arbitrary chosen true component through the removal of the non-true one by the example of two experimental images with two superimposed components. The generalized discrete blurring model and iterative restoration algorithm used are described. It is shown that standard boundary conditions give satisfactory results for «Ant» image but «Pinch» image requires a special image extrapolation approach based on the fitting of reflection direction relative to the boundary for each reproduced structure.
I. Introduction
Resent research in high energy density physics offers a fundamental possibility of developing relatively inexpensive and powerful sources for X-ray radiography. Indeed, X-pinch produced by the electric explosion of crossed metal wires emits soft X-rays of high intensity [1-3]. X-pinch is peculiar in the generation of not one but usually several radiation focuses separated in space. As a result, the X-ray photograph contains several (usually two or three) superimposed image components displaced relative to each other. So, it is necessary to restore an arbitrary selected true component (usually most bright) by removal of the rest, non-true components which make it difficult to recognize the details of reproduced structures. Given space invariance, this image blurring can be described with a convolution model which represents the blurred image as the convolution of the point spread function (PSF) and the true image. Here the PSF is a sum of Dirac delta functions one of which (in the origin of coordinates) reproduces the true component and the rest ones remove (maximally clear) the non-true components:
N
psf (x, y) = 8(x, y) + '£iS(x-xi,y-yi), (1)
i-1
where r, (jc, , y,) is a vector defining the displacement of the / -th non-true component relative to the true one, and N is the number of non-true components. The convolution is usually reversed with iterative deconvolution algorithms, for example, the expectation-maximization maximum likelihood (EMML) algorithm [4]. What makes the convolution model severely disadvantageous is the impossibility to consider boundary conditions which provide a priori information on the image extrapolated beyond the boundaries. Indeed, for correct restoration points near the boundary of a blurred image are likely to have been affected by information outside the field of view. If this information is not taken into account, the results of restoration of images with contrast low-frequency spatial structures near bounda-
ries may contain Gibbs bandpass artifacts which strongly distort the true image. In this work we used the blurring model by Nagy et al. [5] which considers boundary conditions. Restoration results obtained with the residual norm steepest descent (RNSD) algorithm implemented in [6] are provided for two experimental images, «Ant» and «Pinch», taken in X-pinch rays. It is shown that the use of standard boundary conditions allows satisfactory results to be obtained for «Ant» image, while a special approach to image extrapolation is needed to satisfactorily restore «Pinch» image.
2. Blurring model and restoration algorithm
The generalized discrete blurring model by Nagy et al. [5] is described by a system of linear algebraic equations
b = Ax,
(2)
where b and x are vectors representing, respectively, the blurred and the true image, and A is a large, ill-conditioned matrix describing blurring; its non-zero elements are defined by the PSF. The code package by Nagy et al. [6] implements three type of «standard» boundary conditions: zero, periodic and reflexive. The zero boundary conditions correspond to image extension by zeros (Fig. la). The periodic boundary conditions assume that the image is periodically repeated (extended) in all directions (Fig. lb). Finally, the reflexive boundary conditions mean that the image is specularly (i.e., normally) reflected at the boundary (Fig. lc).
0 0 0
0 / 0
0 0 0
/ / /
/ / /
/ / /
frc fr frc
Л f fc
frc fr frc
a)
b)
c)
Fig. 1. Standard boundary conditions: a) zero, b) periodic, c) reflexive (fc is obtained by the transposition of columns f, f, by the transposition of rows f, f,'n by the transposition of rows and columns)
The matrix A is banded block Toeplitz matrix with banded Toeplitz blocks [7] if the zero boundary conditions are used, or the banded block circulant matrix with banded circulant blocks [8] for the periodic boundary conditions, or the sum of banded block Toeplitz matrix with banded Toeplitz blocks and the banded block Hankel matrix with banded Hankel blocks [9] for the reflexive boundary conditions.
The extension of the image domain and the matrix A means that in the iterative inversion of system (2), the result of the matrix-vector multiplication Az, where z is a vector of b (or x) dimension, depends on theelements z corresponding to not only the restored region of the image but also to that one extrapolated in accord with the boundary conditions chosen. The dimension of the extrapolation region is defined by the PSF size. The matrix-vector multiplication is implemented in the package by Nagy et al. with the 2-D discrete fast Fourier transform; its detailed description can be found in [10].
To solve system (2), we have chosen the iterative RNSD algorithm [11] as it is easy to implement and converges much faster than, for example, the EMML algorithm. Moreover, the RNSD algorithm exhibits a semi-convergence behavior [12] with respect to the relative error ||xt -x||/||x|, where xk is approximation to x at the k -th iteration. This means that the depend-
ence of the error on the number of iterations has a global minimum. If iterations are terminated after the minimum has been reached, the solution will be regularized. This is extremely important as the matrix A is ill-conditioned, and the solution requires regularization. Since the blurred initial images are binary and do not take negative values, the algorithm was implemented so as to artificially enforce a nonnegativity constraint on the solution approximation at each iteration. The RNSD algorithm can be represented by the following sequence of steps: x = b
g = Ar(Ax-b)
X = diag(x) y = grXg for k = 1,2,... s = -Xg u = As
a = min(f/u7u,mins.<0(-x, /s,)) \=x+as X = diag(x) z
Г
ATu
-g + az = grXg end.
The function «diag(-)» generates a diagonal matrix containing the initial vector.
Figs. 2 and 3 show the images «Ant» and «Pinch»,
c) d)
Fig. 2. The initial «Ant» image (a) and its restorations with standard boundary conditions: b) zero, c) periodic and d) reflexive
Fig. 3. The initial «Pinch» image (a) and its restorations with standard boundary conditions: b) zero, c) periodic and d) reflexive
Konovalov A.M., Vlasov V.V., Uglov A.S., Shelkovenko T.A., Pikuz S.A.
Iterative restoration of X-ray images taken in X-pinch rays
respectively, obtained in Cornell University on XP generator and restored with the above algorithm. In both cases, the number of iterations is equal to 3. It is seen from Fig. 2 that for «Ant» image, the bandpass artifacts are present if only the periodic boundary conditions are used (Fig. 2c). The zero and reflexive boundary conditions give quite good results (Fig. 2b and 2d). As for «Pinch» image, it is seen that none of the standard boundary conditions gives an artifact-free restoration (Fig. 3). It is clear that some special boundary conditions are needed in this case.
3. Fitting of boundary conditions for «Pinch»
restoration
It is seen from Fig. 3 that the best result of «Pinch» restoration is get with the reflexive boundary conditions, i. e., in one of the three cases where contrast low-frequency structures in the extrapolation region are adjacent to the boundary. So, it is naturally to assume that the structures should be extended beyond the boundaries as it is done for the reflexive boundary conditions. However, the direction of structure reflection relative to the boundaries should differ from the normal direction and be selected separately for each structure.
Information on the reflection direction and on the «useful» size of the extrapolation region can be obtained from, for example, Fig. 3b or 3c. Indeed, using the restored image with distinct bandpass artifacts, for each the j -th boundary structure it is easy to create a vector that defines in which direction and for how many pixels the image should be extrapolated to neutralize the boundary effect (Fig. 4). Thereto it is convenient to produce a segmentation of the boundary structures limited by the artifact strips. Then we can approximate each segmented structure by a quadrangle. Desired vector Ar, coincides with the line segment that connects the midpoints of the quadrangle sides parallel to the boundary.
Let (xmin,xmax,jmin,ymax) be Cartesian coordinates of the boundaries of «Pinch» image. To extrapolate the upper structure on the top of the image, first reflect the initial image relative to the left boundary in the direction Ar, (Fig. 5a) and then relative to the upper boundary in the same direction (Fig. 5b). The corresponding affine coordinate transformations (x, y) -»(xb, yb) can be written as
4 =2*mm
Уь‘У*2(‘™\х)Лу'; . 2QW, -y)*x1 .
І Лд — ЛІ- ,
Луї
І Л — ~Ут<іХ ~ У-
and
(3-4)
To extrapolate the image beyond the right boundary, first reflect the upper structure relative to this boundary in the direction Ar2 (Fig. 5c) and then the lower structure in the direction Ar3 (Fig. 5d).
Fig. 4. Vectors defining the reflection direction and «useful» size of the extrapolation region
Fig. 5. Extrapolation steps for: a) the left boundary; b) the upper boundary; c) and d) - the right boundary
Coordinate transformations have the form as follows
Xfj — 2xmax — x,
2(xmax ~x)Ayj (5)
y* =y +----------A----------'
AXJ
Here j =2 in the first case and j =3 in the second one. The discrete values for pixels should be found through «nearest-neighbor» interpolation. The other values in the extrapolation region including values beyond the lower boundary may be taken zero.
Fig. 6 shows the result of «Pinch» image restoration with the use of the RNSD algorithm and the boundary conditions synthesized as described above. It is seen that the restored image of the true component looks quite satisfactory being free of the bandpass artifacts.
Fig. 6. «Pinch» image restored with the boundary conditions selected as described in the paper
4. Conclusion
Using two experimental X-ray images taken in X-pinch rays as an example, we have shown the possibility of restoring the true component of image blurred due to superposition of non-true components from «spurious» radiation focuses. It is shown that the generalized discrete model by Nagy et al. [5] and the residual norm steepest descent algorithm [11] are quite useable for restoration. However, in the general case, the boundary conditions implemented in Nagy’s code package [6] do not give images that are free of Gibbs bandpass artifacts. These artifacts can be eliminated if use the image extrapolation method described in the present paper. The method is based on the fitting of the direction of image reflection relative to boundaries separately for each reproduced structure. The restoration approach presented seems to be applicable to any image blurred due to spatially invariant superposition of its several components.
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 Grant 2151 and by the U. S. Department of Energy under Grant DE-FG03-98ER54496.
References
1. S.A. Pikuz, T.A. Shelkovenko, V.M. Romanova, D. B. Sinars, D. A. Hammer, S. N. Bland, and S. V. Lebedev, «Xpinch as a source for X-ray radiography», Nukleonika, vol. 46, pp. 21-25, 2001.
2. T.A. Shelkovenko, D.B. Sinars, S.A. Pikuz, K.M. Chandler, and D.A. Hammer, «Point-projection X-ray radiography using an X pinch as the radiation source», Rev. Sci. Instrum., vol. 72, pp. 667-670, 2001.
3. G. V. Ivanenkov, S.A. Pikuz, T.A. Shelkovenko, V.M. Romanova, I. V. Glazyrin, O. G. Kotova, and A. N. Slesareva, «Review of publications devoted to modeling of the electric explosion of thin metal wires», Preprint of Lebedev Physical Institute RAS #9, Moscow, 2004 (in Russian).
4. D.S. C. Biggs and M. Andrews, «Acceleration of iterative image restoration algorithms», Appl. Opt., vol. 36, pp. 1766-1775, 1997.
5. J.G. Nagy, K.Palmer, andL.Perrone, «Iterative methods for image deblurring: a Matlab object oriented approach», Numerical Algorithms, vol. 36, pp. 73-93, 2004.
6. K.P. Lee, J.G. Nagy, andL. Perrone, An Object Oriented Matlab Package «Restore Tools» for Image Restoration, 2002, Avaliable: http://www. mathcs. emory. edu/~nagy/RestoreTools/.
7. J. Kamm andJ. G. Nagy, «Kronecker product and SVD approximation in image restoration», Linear Algebra Appl., vol. 284, pp. 177-192, 1998.
8. H. Andrews and B. Hunt, Digital Image Restoration, Prentice-Hall, Englewood Cliffs, NJ, 1977.
9. M. K. Ng, R. H. Chan, and W.-C. Tang, «A fast algorithm for deblurring models with Neumann boundary conditions», SIAM J. Sci. Comput., vol. 21, pp. 851-866, 1999.
10. . J. G. Nagy and D. P. O'Leary, «Fast iterative image restoration with a space-varying PSF», Proc. SPIE, vol. 3162, pp. 388-399, 1997.
11. L. Kaufman, «Maximum likelihood, least squares, and penalized least squares for PET», IEEE Trans. Med. Imag., vol. 12, pp. 200-214, 1993.
12. J. G. Nagy and K. Palmer, «Steepest descent, CG, and iterative regularization of ill-posed problems», BIT, vol. 43, pp. 1003-1017, 2003.
