IMPROVEMENT OF THE QUALITY OF RESTORING IMAGE WITH RADON INVERSE TRANSFORM

K.Ya. Kudryavtsev

National Research Nuclear University MEPhI (Moscow Engineering Physics Institute)

KYKudryavtsev@mephi.ru

Contetns

1. Introduction

2. Basic relation

2.1. Algorithm description

2.2. Radon transform matrix value increasing algorithm.

3. Estimation of efficiency

4. Conclusion

References

Annex. MatLab 7.0 script for implementing the developed algorithm

Abstract

The creation of the image of internal contents of object without its physical corrupting applied in medicine, metallurgy and crystallography is based on X-ray radiation with the subsequent mathematical processing of received data. When the object is being scanned, the output ray intensity is equal to density distribution function along the ray path. Registration of emission calculated at different angles allows to obtain the original data used to reconstruct the image.  Data used for image reconstruction has the form of matrix, each element of which is the result of integral evaluation along the determined line.  The bigger matrix is, i.e. the more experimental data is available, the more exact image reconstruction will be possible. However, the number of experiments is often limited by medical, technical, economic and other reasons. As a mathematical tool for the reconstruction of the cross-sectional image of an object, a direct and inverse Radon transform is used. This paper proposes the algorithm of rising the quality of reconstructed image from Radon transform matrix by adding additional colums containing Radon transform in interjacent angle turn calculated with the help of developed approximation formulas. Doubling the number of Radon integral transformation matrix colums helps to significantly rise the quality of image being reconstructed with Radon inverse transformation.

1. Introduction

The construction of the image of the internal contents of the object without its physical destruction, applied in medicine, metallurgy and crystallography is based on X-ray radiography, followed by mathematical processing of received data. The mathematical theory of image reconstruction is based on the direct and inverse integral Radon transformations [1-3]. The methodology of obtaining a two-dimensional laminogram  includes two stages [4]. At the first stage, the projection data is formed, on the second stage - the image of the cross-section is reconstructed from the projection data. Modern computer tomographs mainly use systems with a fan-shaped or parallel beam. In the fan-shaped system, the rays diverge from the X-ray source in one subspace, but at different angles forming a fan. The radiation detectors are located on the arc on the other side of the studied object. In a parallel system (it will be analyzed hereinafter) a line-up of radiators is used, which forms parallel beams.

When the object is being scanned, the output intensity of the ray depends on the density of the studied object and is equal to the integral of the distribution function of the matter density along the ray trajectory. The registered radiation (Radon image or projection) measured for a set of parallel beams at different angles of rotation is represented as matrix R with the N × M dimension, each column of which corresponds to a certain rotation angle, and the number of lines is determined by the number of detector emitters. Actually, the reconstruction of the image is carried out using the inverse Radon transform operation (the inverse projection method [5]), the initial data for which are the values of the matrix R.

The larger the dimension of the matrix R is, i.e. the more detailed the studied object is scanned, the more accurate the image reconstruction will be. However, the high dimensionality of the R matrix is associated with a number of technical, medical and economic constraints. Therefore, it is desirable to have an algorithmic methodology for calculating the additional values of the matrix R on the basis of the available experimental data; that is where the issue of "missing data reconstruction” arises. In fact, the data is not missing. In order to improve the image recovery, it is desirable to increase the number of columns of the matrix R, thereby decreasing the discretization interval of the angle of rotation of the radiators line-up.

This paper [6] presents a comparison of various methodologies [7-13] for missing data recovery, such as:

·         filling with the general average,

·         filling with selection,

·         filling by regression,

·         spline interpolation methods,

·         methodology of maximum likelihood filling,

·         methods of factor analysis,

·         methods of cluster analysis,

·         algorithms of the family Zet (Wanga),

·         methods of using neural networks.

The results of a comparative analysis [6] allow us to conclude that the best recovery results among all methods are shown by using the method of the missing value by spline interpolation by the present elements.

Consequently, in order to calculate the values of elements of additional columns of the matrix R, it is rational to use the approach based on the approximation methods.

This paper proposes an algorithm for improving the quality of image reconstruction by adding additional columns containing the Radon transform in the intermediate rotation angles calculated with the derived approximation formulas. The twofold increase in the number of columns of the matrix of the Radon integral transformation makes it possible to significantly improve the quality of image reconstruction with the inverse Radon transform.

2. Basic relation

The main idea of the proposed approach was previously presented by the author in paper [14]. Here we propose a more detailed and rigorous justification of mathematical relations.

Let us assume that there is Radon transform  of some image I. Normally Radon  transform  is shown as  matrix.

– is determined by image I size and line discretization interval τ along which Radon integral transform.

– is determined by the number of angle rotation with  step (within 0 - 180 angle degrees) of rotating coordinates, inside of which lines for evaluating Radon transformation are drawn.

Matrix element  is equal to the integral, evaluated along the line of perpedicular axis of x-coordinate of rotating coordinate system and the beforementioned integral is situated on a distance   from the origin of coordinates. Angle rotation is equal to

,

where  – is a turn step of moving coordinate system. Line formula is evaluated with the following formula [1]

As shown in [5], according to Kotelnikov–Shannon formula in order to reconstruct the image in an acceptable quality the following correlation should be performed:

Hence increasing (doubling) of   value is of special interest, i.e. dicreasing of angle rotation step in the rotating coordinate system.

The simplest way to get  R-matrix is to add one column between each  -matrix’s columns and evaluate  elements with line approximation formulas:

However, as it will be shown below, the following methodology will not result in quality rising of the reconstructed image, so therefore it is required to develop a special algorithm, that will twice increase  value. In other words, if Radon source transform of dimention is evaluated with  angle rotation step, that means the algorithm should result in   Radon transform matrix of  dimention, that is to obtain a Radon transform evaluated with  angle rotation step.

2.1. Algorithm description

Let us analyze and  coordinate systems, whose center is connected with the center of the studied object and rotated relative to each other by an  angle. (see Pic.1).

Pic.1  and  coordinate systems turned through  angle relative to each other.

In  coordinate system the integral is being calculated along AA’ line on s=nτ distance from the origin of the coordinates. In , coordinate system turned through  angle about , the integral is evaluated along BB’ line also situated on s=nτ distance from the origin of the coordinates. Let us analyze  coordinate system turned through  angle about . Please note that  coordinate system is turned through  angle about .

It is required to estimate the value of the integral along CC’ line, which is perpendicular to  axis and which is situated on s=nτ distance from the origin of the coordintes, i.e. to evaluate Radon transform with  angle rotation step two times smaller, that the original .

Let us mark the integral value along AA’, BB’, CC’ lines with RAA’ ; RBB’ ; RCC’ accordingly. Integral RCC evaluation with formula

is incorrect, as CC’ line lies mostly out of BDA’ and ADB’ sectors. This averaging goes to RDD’ integral evaluation, i.e. to the integral along DD’ line, which goes between AA’ and BB’ lines and lies entirely within BDA’ and ADB’ sectors.

Hence let us assume that

The distance from the origin of the coordinates, where DD’ line lies is equal to

For the following transform let us introduce the following notation:

where . First argument of  function means the distance from the beginning of the coordinate origin, but the second angle of rotation. So, for example,   means, that the integral is evaluated along the line, which is perpendicular to  axis, and rotated through the angle equal to zero, and is situated on  distance from the coordinate origin.

For distance  the same formulas can be written, especially for EE’ line

If values  and  are evaluated, that means  value can be evaluated with the help of line approximation if CC’ line is between EE’ and DD’ lines.

For that the following formula should be performed

or

If the inequaltion on the right is always performed, fot the inequation on the left the following can be written

or

Maximum value for  is equal to . Hence

or

So the following restriction for  angle rotation connected to  value is obtained.

Assuming that the last ineqitation for  angle rotation is done, the following formula for  with the help of line approximation is obtained

By simplifying the formula the result is

And taking into consideration, that

The final result will be

The obtained formula is rather interesting. There is a summand in it

which can be estimated as a simple line approximation like  and some correction

That estimates line position, along which there should be an integration with reference to lines used for evaluating approximation.

It should be mentioned that, if

i.e. the correction is not used, because CC’ direct line is entirely situated within BDA’ and ADB’ sectors and divides them into halves.

2.2. Radon transform matrix value increasing algorithm.

Given data:  Radon transform  matrix. Angle rotation step is equal to , i.e. integral calculation will be done for the angles: , while .

Let us divide the following algorithm into several stages.

Stage 1. Create  -matrix and fill it with zeros.

Stage 2. Assign R-matrix odd columns to the corresponding columns of -matrix.

Stage 3. Evaluate R-matrix even columns values for  with the following formulas:

where

Stage 4. Evaluate  column value of R-matrix for  with the following formula:

where

Stage 5. Evaluate R-matrix even columns values for  with the following formula:

where

Stage 6. Evaluate R-matrix  column value for  with the following formula:

where

3. Estimation of efficiency

Estimation of efficiency of the proposed algorithm will be held in MatLab 7.0. 256 ő 256 px phantom image is used fot the test. (see Pic. 2)

Pic. 2. Original phantom image

Let us choose angle rotation step as  and  calculate Radon transform matrix with the help of radon function in MatLab 7.0. As a result . -matrix is obtained. Here, .

When performing Radon inverse transform with iradon function, the restored image is obtained (see  Pic.3)

Pic. 3. Restored phantom image (angle rotation step )

As we can see the restored image quality is not very high, as  angle rotation step is too large.

If angle rotation step is  and Radon transform   matrix is evaluated and than Radon inverse transform, the restored image has better quality. (see Pic.4)

Pic. 4. Restored phantom image (angle rotation step )

Choosing  -matrix as given source, let us make Radon transform for  R-matrix using the abovementioned algorithm with formulas (1)-(5). Using Radon inverse transform, the image like on Pic. 5 is obtained.

Pic. 5. Phantom image restored from approximation matrix.

Comparing Pictures 3, 4 and 5 we can conclude, that the proposed algorithm improves the quality of the reconstructed image. Pic. 5 is visually less sharp than Pic. 4, but its quality is better, than that of Pic. 3.

In order to emphasize the significance of the proposed approximmation correction  a Radon transform approximation matrix without correction should be built and than restored. (see Pic.6).

Pic. 6. Image restored from approximation matrix without approximation correction.

As it can be seen, ignoring the approximation correction do not lead to the increase of the image quality.

In addition to the subjunctive, visual estimation of the proposed algorithm effectivity, the mean absolute deviation from exact values of the evaluated approximation values should be found. As exact values, we take Radon transform matrix of dimention  with an angle rotation step . The estimation is calculated according to the following formula:

where  – Radom transform matrix elements with angle rotation step.

– approximation transform matrix elements, evaluated according to the previously discribed algorithm with formulas (1)-(5).

The evaluation result is . To put this in perspective, if approximation matrix is built without corrections, the mean absolute deviation is . According to Table 1 if lesser angle rotation steps are used ()  value becomes even less.

Table 1. Comparison of approximation matrix evaluation accuracy

 11.2847 0.3261 34.6 11.2554 0.2969 37.9 11.1721 0.2138 52.25 11.1182 0.1598 69.57 11.0718 0.113 97.98

The proposed approximation correction significantly raises approximation quality, that is a substancial result.

The implementation of the described algorithm is presented in the application as a MatLab 7.0 script.

4. Conclusion

The following paper proposes an algorithm for raising the quality of the reconstructed image from Radon transform matrix by adding additional columns with Radon transform in the intermediate rotation angles, calculated from the derived approximation formulas.

The proposed algorithm consists of six stages:

- At the first stage: an R matrix with double number of columns is created and filled with zeroes;

- In the second stage: odd R-matrix columns are assigned to the corresponding columns of original matrix;

- At subsequent stages: the values of the elements of the matrix R are calculated from formulas (2) - (5).

Experimental validation showed, that the proposed algorithm allows to significantly improve the reconstructed image quality with the help of Radon inverse transform.

All the formulas obtained were experimentally checked in MatLab 7.0 system, which confirmed their full correctness.

References

1. Natterer F. The Mathematics of Computerized Tomography (Classics in Applied Mathematics, 32), Philadelphia, PA: Society for Industrial and Applied Mathematics, 2001 ISBN 0-89871-493-1

2. The Physics of Medical Imaging. Edited by Steve Webb Adam Hilger, Bristol and Philadelphia, 1991. ISBN 0-885274-367-0

3. Kudryavtsev K.Ya., Mikhaylov D.M. Radon integral transform on sparse rectangular grid. International Journal of Tomography and Simulation ISSN 2319-33362015, V. 28, N. 1, p. 89-104

4. Gruzman I.S., Kirichuk V.S., Kosyh V.P., Peretjagin G.I., Spektor A.A. Digital image processing in information systems: Textbook. NSTU, 2002. 352 p.

5. Kak A.C, Slaney M. Principles of computerized tomographic imaginng (Classics in Applied Mathematics, 33), Philadelphia, PA 19104-2688: Society for Industrial and Applied Mathematics, ISBN 0-89871-494-X

6. Abramenkova I.V., Kruglov V.V. Methods for recovering gaps in data sets. - Software products and systems, 2005, ą2, p.18-22

7. Little R.J.A., Rubin D.B. Statistical analysis of data with omissions. - Moscow: Finance and Statistics, 1990.

8. Zagoruiko N.G. Applied methods of data and knowledge analysis. - Novosibirsk: Publishing house of the Institute of Mathematics, 1999.

9. Zagoruiko NG, Elkina VN, Timerkaev VS An algorithm for filling out blanks in empirical tables (Zet algorithm). Empirical prediction and pattern recognition. Novosibirsk, 1975. Issue. 61: Computing systems. p. 3-27.

10. Rossiev AA Modeling data using curves to restore gaps in tables. Methods of Neuroinformatics. Krasnoyarsk: KSTU, 1998.

11. Demidenko E.Z. Linear and nonlinear regression. Finance and Statistics, 1981.

12. Dvienko S.D. Non-hierarchical divisible clustering algorithm. Automation and telemechanics. 1999. ą 4. P. 117-124.

13. Kruglov VV, Borisov VV Artificial neural networks. Theory and practice. M .: Hot line - Telecom, 2002.

14. Kudryavtsev K.Ya.  Rising the Quality of Restor-ing Image with Radon Inverse Transform Basing on Spectral Approximation Formulas. International Journal of Tomography and Simulation ISSN 2319-3336 2016, V. 29, N. 2, p. 50-61.

Annex. MatLab 7.0 script for implementing the developed algorithm

P=phantom(256);

imshow(P);

dtheta=10;

theta1=0:dtheta:(180-dtheta);

theta2=0:dtheta/2:(180-dtheta/2);

Ntau=size(R1, 1);

NZero=(Ntau+1)/2;

Ntheta=size(R1, 2);

Ntheta2=2*Ntheta;

R1Appr = zeros(Ntau,Ntheta2);

R1ApprSimple = zeros(Ntau,Ntheta2);

Nel=size(R1Appr, 1) * size(R1Appr, 2);

Ktheta=(1-cos(pi*dtheta/360))/2;

for i=1:Ntheta

R1Appr(:,2*i-1)=R1(:,i);

end

for i=1:Ntheta

for n=1:NZero

j=NZero+n-1;

I0Nt=R1Appr(j,2*i-1);

I0Nt1=R1Appr(j-1,2*i-1);

if i < Ntheta

I2Nt=R1Appr(j,2*i+1);

I2Nt1=R1Appr(j-1,2*i+1);

else

I2Nt=R1Appr(j,1);

I2Nt1=R1Appr(j-1,1);

end

R1Appr(j,2*i)=(I0Nt+I2Nt)/2-(I0Nt-I0Nt1+I2Nt-I2Nt1)*(n-1)*Ktheta;

R1ApprSimple(j,2*i)=(I0Nt+I2Nt)/2;

end

end

for i=1:Ntheta

for n=1:NZero

j=NZero-n+1;

I0Nt=R1Appr(j,2*i-1);

I0Nt1=R1Appr(j+1,2*i-1);

if i < Ntheta

I2Nt=R1Appr(j,2*i+1);

I2Nt1=R1Appr(j+1,2*i+1);

else

I2Nt=R1Appr(j,1);

I2Nt1=R1Appr(j+1,1);

end

R1Appr(j,2*i)=(I0Nt+I2Nt)/2-(I0Nt-I0Nt1+I2Nt-I2Nt1)*(n-1)*Ktheta;

R1ApprSimple(j,2*i)=(I0Nt+I2Nt)/2;

end

end

Pogr=sum(sum(abs(R1Appr-R2)))/Nel

PogrSimple=sum(sum(abs(R1ApprSimple-R2)))/Nel