PSNR BASED OPTIMIZATION APPLIED TO ALGEBRAIC RECONSTRUCTION TECHNIQUE FOR IMAGE RECONSTRUCTION ON A MULTI-CORE SYSTEM

The present work attempts to reveal a parallel Algebraic Reconstruction Technique (pART) to reduce the computational speed of reconstructing artifact-free images from projections. ART is an iterative algorithm well known to reconstruct artifact-free images with limited number of projections. In this work, a novel idea has been focused on to optimize the number of iterations mandatory based on Peak to Signal Noise Ratio (PSNR) to reconstruct an image. However, it suffers of worst computation speed. Hence, an attempt is made to reduce the computation time by running iterative algorithm on a multi-core parallel environment. The execution times are computed for both serial and parallel implementations of ART using different projection data, and, tabulated for comparison. The experimental results demonstrate that the parallel computing environment provides a source of high computational power leading to obtain reconstructed image instantaneously.


Introduction
Image reconstruction methods are central to many of the new applications of medical imaging such as Positron Emission Tomography (PET), Computed Tomography (CT), Magnetic Resonan-ce Imaging (MRI) and Electron Magnetic Reso-nance Imaging (EMRI).They are most commonly used to visualize detailed internal structure and limited function of the object of interest.
Image reconstruction is a mathematical process that generates images from projection data acquired at many different angles around the object of interest.The projections are collected by sweeping the magnetic field at projection angles defined by the magnetic field gradient directions [1,2].To perform image reconstruction, the projections   (), collected along a set of fieldgradient orientations in polar coordinates, are used to obtain the image f(x, y) [3] as given in the equation (1).

𝑓(𝑥
Here r is taken on the x-y plane such that  =  + , and   * () is the projection   () A. Bharathi Lakshmi and D. Christopher Durairaj, PSNR Based Optimization 87 filter according to the expression inside the square brackets [3].Image reconstruction has been carried out using different types of reconstruction algorithms [4,1].Reconstruction methods utilize projection data as input and generate the estimate that resembles the internal structure as output [5,6].Data sets with 36 projections measured from 0 0 to 180 0 around the phantom object were considered in the present study.The same data set was used for testing the capability of the algorithms from restricted number of projections, by skipping projections at uniform angular distribution.
Reconstruction of images is usually done in two ways: Analytical and Iterative.Analytical method such as Back Projection (BP) or Filtered Back Projection (FBP) is used for different imaging modalities such as CT and PET in clinical settings because of its speed and easy implementation [3].For noisy projection data as well as for limited number of projections, the FBP method of image reconstruction shows very poor performance.Hence currently there is considerable interest to evaluate the use of other reconstruction methods for medical imaging techniques [6].FBP algorithm produces high-quality images with excellent computational efficiency.However, FBP produces low Signal-to-Noise Ratio (SNR) images when limited number of projections is used [12].
An Iterative method using a non-linear fit to the projection data has shown to give ripple free images [7].Iterative Methods are based on optimization strategies incorporating specific constraints about the object and the reconstruction process.The iterative reconstruction techniques perform better than the FBP method when reconstruction is attempted with limited number of projection data [3].Some of the accepted iterative algorithms are Additive Algebraic Reconstruction Technique (AART) and Multiplicative Algebraic Techniques (MART) [12].
However, the quality of the reconstructed images obtained from AART algorithm depends on number of iterations.Based on the number of available number of projections and the size of the phantom, the number of iterations differs.It is therefore necessary to find the best iteration in order to exploit correctly the promising iteration based on a better Peak to Signal Noise Ratio (PSNR) of the reconstructed images.Based on the equation(1) an optimization program has been developed for the given data set.The best PSNR value is obtained and verified whether the same PSNR value is achieved even after the selected iteration.
Parallel computing is emerging as a principle theory in high performance computing [14].In recent years, parallel computing with massive data has emerged as a key technology in imaging techniques also.Shared memory parallelization has been proved to be a best way to attain better runtime performance recently for image reconstruction [15].A shared-memory multiprocessor (SMP) consists of a number of processors accessing one or more shared memory modules.The penalty of using inter-processor communication is not up to the mark on SMP compared to distributed memory architectures [15].For a relatively large data size, it is advantageous to use SMP architecture.It has also been shown that shared memory parallelization is more suitable than distributed memory parallelization for image processing tasks and leads to better throughput as most of the computers now have two or more processors which share the memory [16].These features have motivated us to perform the parallelization of Algebraic Reconstruction Technique (ART) on a SMP parallel architecture.The present study focuses on reducing the computational complexity of ART using parallel programming techniques.Section 2 describes about ART briefly.The design and implementation ART algorithm in both parallel and sequential version are given in section 3 Section 4 discusses the results.

Radon Transformation
The main application of image reconstruction from projection technique is mostly related to medical image processing.The Procedure to implement Image Reconstruction from Projection (IRP) technique in the practical applications are "scanning" or "data acquisition" is considered to be the first and the very important step [17].Such data acquisition is done by means of PET, CT, MRI or EMRI in a procedure by passing rays in specific intervals of angles.
The Radon Transformation is a fundamental tool that computes projections of an image matrix along specified directions [18].The 2D Radon transformation is the projection of the image density along a radial line oriented at a specific angle.The value of a 2-D function at an arbitrary point is uniquely obtained by the integrals along the lines of all directions passing the point.The Radon transformation shows the relationship between the 2-D object and its projections.Figure 1 shows a 2-D function (, ).Integrating along the line, whose normal vector is in θ direction on s axis results in the (, ) projection represented in equation (2).The points on the line whose normal vector is in θ direction and passes the origin of (, )-coordinate satisfy the equa-tion  +  = 0.The general equation of the radon transform is acquired as where  is zero for every argument except to 0 and its integral is one [19].
The projection data obtained thus from Radon Transformation is utilized as input by the Reconstruction algorithm that produce estimates of the original internal structure as output [5,20].The size of data sets acquired by the different imaging modalities are usually huge because of the complex data type of the raw collection data, multiple gradients in the experiments, high dimen-sions of the resultant 3-D images, higher k-space requirement of whole body imaging and the number of points collected from the imager.The iterative methods, hence, suffers more reconstruction time.

Algebraic Reconstruction Technique (ART)
Image reconstructions based on Iterative methods create two-dimensional images from scattered or incomplete projections such as the radiation readings acquired during a medical imaging study.Algebraic Reconstruction Technique (ART) falls under the category of Iterative methods.
ART is one of the methods used for solving the linear system which appears in image reconstruction.ART can be broadly classified as either sequential or simultaneous or block itera-tive [21].ART is a fully sequential method and has a long history and literature.Originally it was proposed by Kaczmarz [22], and independently, for use in image reconstruction by Gordan, Bender and Herman [23].The vector of unknowns is updated at each equation of the system, after which the next equation is addressed.If system of equation is (0.1) consistent, ART converges to a solution of this system.If the system is incon-sistent, every subsequence of cycles through the system converges, but not necessarily to a least square solution [24].
ART perform corrections during iterations, without increasing the computation time.The image (, ) is a continuous two dimensional function and an infinite number of projections are mandatory for reconstruction [12].In practice (, ) is calculated using a finite number of points   ( = 1, 2, 3, … , ) where N represents the total number of cells, from a finite number of projections as shown in figure .2.
In figure 2 a ray is a fat line running through the (, )-plane where each ray is of width r.A line integral is called a ray-sum represented as   measured with i th ray as shown in figure .2.
The relationship between the   ′  and   ′ may be expressed as where M is the total number of rays (in all the projections) and wij is the weighting factor that represents the contribution of the j th cell to the i th ray integral and   represents a set of matrix equation for the data point   .The expanded form for equation(3) for the j th sample is given by Equation( 4) can also be expressed in the form of algebraic equations as Here,   is the weighting factor that represents the contribution of the  ℎ cell to the  ℎ sample sum and   represents a set of matrix equations for the data point   .Most of the   in Eqn. 4 is zero since only a small number of cells contribute to any given ray-sum.The density values   are iteratively adjusted until the calculated projections agree with the measured projections [12].Each projected density is thrown back across the reconstruction space in which the densities are iteratively modified to bring each reconstructed projection into concur with the measured projection [25].The projection data set is sustained in a vector and a weight sparse matrix   is constructed.Every row in   sparse matrix may contain  +  − 1 (where  x  is the resolution).As every row stands for the length of the segments obtained by the intersection of ray with the grid, and all reconstruction algorithms use rows of sparse matrix, the best method to store this matrix is in compressed row storage [17].
For each sample, the correction coefficient is computed as: . The average value of the correction coefficient is calculated.Correction is applied for each cell j as given: where λ is the relaxation parameter.This procedure is iteratively performed for all of the projection angles.As the size of the data set increases, the computation time increases.

OpenMP Architecture and Directives
Parallel computing is a form of computation in which many calculations are carried out simultaneously; large problems are divided into smaller ones, solved concurrently.The parallelism can be applied in image processing applications by three main ways: 1) Data Parallel 2) Task Parallel and 3) Pipeline Parallel.In Data Parallel approach, the data is divided and distributed among the computing units.The data parallelism to image data can be applied using one of three basic ways: i) Pixel Parallel ii) Row or Column parallel and iii) Block Parallel [27].This algorithm is parallelized in row/column parallel.In task parallel, image processing instructions/low level operations are grouped into tasks and each task is assigned to a different computational unit.If image processing application requires multiple images to be processed, then pipeline processing of images can be done [28].pART is implemented using OpenMP parallel computing in C language.OpenMP is a programing model for SMP computer systems.Data in memory can either be shared between all threads or private for one thread.Data transfer between threads is transparent to the programmer.OpenMP uses fork-and-join model of parallel execution.The program written with OpenMP begins execution as a single-process, called the master thread.The master thread executes the current program sequentially until it bump into parallel directives such as #pragma omp.The master thread forks a number of worker threads when it enters a parallel region.A parallel region is a block of code that is executed by all threads concurrently.
The "parallel for" or "for" is a wok sharing directive that distributes the workload of a "for" loop among all the threads.Data sharing of variables is mentioned at the beginning of the parallel region or work sharing construct using the SHARED or PRIVATE Clauses.

Data Set
The reconstruction system uses Shepp Logan phantom data of different sizes such as 64, 128 and 256.The figure.The projections of the Shepp Logan phantom in various angles are plotted in the figure 5.This is obtained by using radon function in matlab passing at which specific angles the object should be rotated.This system uses five different angles, such as 6 0 , 9 0 , 12 0 , 15 0 , 18 0 obtaining 30, 20, 15, 12, 10 numbers of projections respectively.Rows 1, 2 and 3 of figure.5 refer to the projections taken from the images sizes 64x64, 128x128 and 256x256 respectively.Columns A, B, C, D, E refer to the 10, 12, 15, 20 and 30 projections taken in 18, 15, 12, 9 and 6 angles respectively.

Pseudo Code
The pseudo code for ART algorithm implemented  for(all the rows) { correct the error by multiplying the difference with the calculated data.
apply the correction.} } recursively call part function for remaining projections } in MEX function executed sequentially is given in figure 6.
The pseudo code for pART algorithm implemented in MEX function executed parallel shows in figure 7.

UML Diagram
The operation of ART in sequential and parallel is symbolised in the figure 8 The data is read from the corresponding number of projections.This data is supplied into the MEX function to execute under single and multiple processors.For each projection the error value is calculated and the correction is applied Initially, the program starts with initialization at each core.Then the calculated and measured projection data is co-distributed between the workers.After that, each worker calculates the correction and applies the correction as per the algorithm till all the projections are completed.Then the processed data is collected in a vector  when the parallelism ends.

Results and Analysis
The results of constructing Shepp Logan Phantom image using ART in both sequential and parallel is given in figure 10.In this work, the time complexity of the phantom image of different size (64, 128 and 256) is compared in 2, 4 and 8 cores.
Peak-signal-to-noise ratio (PSNR) is used as a metric to check perceptual similarity between the and reconstructed images.The PSNR value measured in db is tabulated in table 1.According to Chen et al (1998), PSNR above 40 db indicates a good perceptual fidelity.It can be observed that PSNR for the different size of images using various angles is above 60 db which indicates the excellent perceptual fidelity.
In figure 11 the PSNR value of the reconstructed image using ART for various sizes in different number of projections is graphed.
Reconstruction time taken by the Algebraic Reconstruction Technique for different size of phantom image in sequential and parallel using 2, 4 and 8 cores in an AMD Processor under LINUX platform.The time complexity of the reconstructed image of various sizes under 2, 4 and 8 cores is given with respect to the number of projections.
Reconstruction time taken by the Algebraic Reconstruction Technique for different size of phantom image in sequential and parallel using 2, 4 and 8 cores in an AMD Processor under LINUX platform.The time complexity of the reconstructed image of various sizes under 2, 4 and 8 cores is given with respect to the number of projections.The optimized number of iteration to reconstruct an image in the three represented sizes at 30, 20, 15, 12 and 10 number of projections is tabulated in Table 2 and plotted in figure 12.
In figure 13, 14 and 15 the time complexity of phantom image of size 64, 128 and 256 reconstructed using 2, 4 and 8 cores with respect to 30, 20, 15, 12 and 10 is plotted respectively.Table 3, 4 and 5 tabulates time complexity for 64, 128 and 256 size images respectively.Table 3, 4 and 5 shows the reconstruction time taken by 1 Core (row1), 2 core (row2), 4 core (row 3), and 8 core (row 4) when using 30, 20, 15, 12 and 10 projections in the ART for image size 64, 128 and 256 respectively.It is observed that the time gradually reduces as the number of cores increases, for a given sets of projections.A graph is plotted to show the performance of the parallel system.It elucidates the time complexity of the system for a given number of projections using different cores of the parallel processor.
The time complexity of the system implementing parallel processor has got a considerable reduction of time consumptions which is certainly a high degree of utility to the user.A minor change in the time consumption will have a revolutionary impact while it is employed.The specific value of this finding is that the maximum number of core reconstruct the image is fast even for minimum number of projections.

Conclusion
The number of iteration mandatory to reconstruct an image is optimized.The images are reconstructed sequentially as well as in parallel environment using different projection data sets.In this study, of Shepp Logan Phantom data is reconstructed using ART and pART.The results have shown encouraging indication of the efficiency of the parallelization of ART algorithm.In general, the pART algorithm gives a paramount computational efficiency better than ART.The computational efficiency of both ART and pART is reported in this article.

Figure. 6 .
Figure. 6. Art algorithm art() { if ( not yet reached all the projections) { for (all elements in the projection) { calculate the value by multiplying the vector and the calculated data in the corresponding projection } calculate the error by subtracting the measured data from the calculated value for(all the rows) { correct the error by multiplying the difference with the calculated data.apply the corrected value to the vector.} } recursively call art function for remaining projections } (a) and figure 8(b) respectively.Parallel activity is pictured as Fork and Join.

Figure. 10 .
Figure.10.The Reconstructed Shepp Logan Phantom.Rows 1, 2 and 3 refer to the 64x64, 128x128, 256x256 size of the Image respectively.Columns A, B, C, D and E refer to the reconstructed image from the 10, 12, 15, 20 and 30 projections of an image taken in 18, 15, 12, 9 and 6 angles in Sequential and parallel.

TABLE 4 TIME
Figure. 15.A graph showing the Time Complexity of reconstructing Phantom image of size 256x256 parallel in 2, 4 and 8 core with respect to projections.