Method and apparatus for ultra fast symmetry and SIMD based projectionback projection for 3D pet image reconstruction
Download PDFInfo
 Publication number
 USRE42034E1 USRE42034E1 US12498952 US49895209A USRE42034E1 US RE42034 E1 USRE42034 E1 US RE42034E1 US 12498952 US12498952 US 12498952 US 49895209 A US49895209 A US 49895209A US RE42034 E1 USRE42034 E1 US RE42034E1
 Authority
 US
 Grant status
 Grant
 Patent type
 Prior art keywords
 θ
 projection
 φ
 image
 data
 Prior art date
 Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
 Active
Links
Images
Classifications

 G—PHYSICS
 G01—MEASURING; TESTING
 G01T—MEASUREMENT OF NUCLEAR OR XRADIATION
 G01T1/00—Measuring Xradiation, gamma radiation, corpuscular radiation, or cosmic radiation
 G01T1/29—Measurement performed on radiation beams, e.g. position or section of the beam; Measurement of spatial distribution of radiation
 G01T1/2914—Measurement of spatial distribution of radiation
 G01T1/2985—In depth localisation, e.g. using positron emitters; Tomographic imaging (longitudinal and transverse section imaging; apparatus for radiation diagnosis sequentially in different planes, steroscopic radiation diagnosis)

 G—PHYSICS
 G06—COMPUTING; CALCULATING; COUNTING
 G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
 G06T11/00—2D [Two Dimensional] image generation
 G06T11/003—Reconstruction from projections, e.g. tomography

 G—PHYSICS
 G06—COMPUTING; CALCULATING; COUNTING
 G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
 G06T11/00—2D [Two Dimensional] image generation
 G06T11/003—Reconstruction from projections, e.g. tomography
 G06T11/006—Inverse problem, transformation from projectionspace into objectspace, e.g. transform methods, backprojection, algebraic methods

 A—HUMAN NECESSITIES
 A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
 A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
 A61B6/00—Apparatus for radiation diagnosis, e.g. combined with radiation therapy equipment
 A61B6/02—Devices for diagnosis sequentially in different planes; Stereoscopic radiation diagnosis
 A61B6/03—Computerised tomographs
 A61B6/037—Emission tomography

 G—PHYSICS
 G06—COMPUTING; CALCULATING; COUNTING
 G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
 G06T2211/00—Image generation
 G06T2211/40—Computed tomography
 G06T2211/424—Iterative
Abstract
Description
The present application claims priority from Korean Patent Application No. 1020060042155 filed on May 10, 2006, Korean Patent Application No. 1020070027305 filed on Mar. 20, 2007, and Korean Patent Application No. 1020070040515 filed on Apr. 25, 2007, the entire subject matter of which is incorporated herein by reference.
1. Field
The present invention relates to a method and apparatus for ultra fast symmetry and SIMD based projectionbackprojection for 3D PET image reconstruction. More specifically, the present invention relates to a method and apparatus based on the symmetry properties of the projection and backprojection processes, especially in the 3D OSEM algorithm that requires a plurality of projections and backprojections.
2. Background
There has been a remarkable progress in PET development over the recent years, especially in the areas of hardware, software and computer implementation of image reconstruction. Recent developments in PET scanners (e.g., HRRT (High Resolution Research Tomograph) developed by CTI (now Siemens)) allow greatly enhanced resolution and sensitivity. In such PET scanners, the amount of collected coincidence line data contains more than 4.5×10^{9 }coincidence lines of response generated by as many as 120,000 nuclear detectors. Such large amount of data and the reconstruction of this data set pose to be a real problem in HRRT. That is, they pose to be major problems in achieving further developments and applications of high resolution PET scanners. Thus, in such types of scanners, obtaining one set of reconstructed images often requires many hours of image reconstruction. For example, in HRRT with full data collection in normal brain scans (using SPAN 3), the image reconstruction time is almost eighty minutes. This makes it practically impossible to attempt any list mode based dynamic imaging since the image reconstruction time takes many days (as long as 43 hours or more for 32 frame dynamic image reconstruction).
In general, tomographic images can be reconstructed by two approaches, one being an analytic method and the other being an iterative approach. PET scanners of different types were developed in the mid 1970's and the application of various tomographic image reconstruction techniques was naturally introduced in the field. In case of an analytic approach such as backprojection and filtering or filtered backprojection (FB), an artifact known as a streak artifact is frequently generated. This is especially true when the detector arrangements are not uniform such as in the case of HRRT (e.g., Siemens' High Resolution Research Tomograph) where the detectors are arranged in a set of blocks in an octagonal shape. These types of detector arrangements often involve missing data due to the gaps between the blocks and result in a severe streak artifact in the case of the FB technique. Therefore, alternative approaches such as an EM (Expectation Maximum) algorithm have been sought. Generally, the EM approach requires several steps in the reconstruction process, the two major steps of which are: projection (forward projection) to create projection data from the image or object and backprojection into the image domain for the final image reconstruction. In the EM algorithms, these two processes are repeated until satisfactory images are obtained. Obviously, these repeated projection and backprojection processes are time consuming and have been the major drawback of the EM approach compared to straight filtered backprojection (FB) algorithm. In addition, in case of 3D image reconstruction, the computational burden increases out of proportion due to the astronomical increases in the coincidence lines or the line of responses (LOR). This is a major stumbling block in the daily operation of high resolution PET scanners. Thus, there is a strong need for improving the computational speed or the reconstruction time in EM approaches, especially with high end PET scanners such as HRRT.
Projection methods usually employ a system matrix, which is determined by the geometric factor of the scanner. As the resolution of the PET image improves and the number of slices increases, the size of the matrix is also increased drastically in proportion to the increases in LORs, thus resulting in not only the need for a large memory but also the total computation time. Current HRRT, for example, requires nearly eighty minutes of reconstruction time, in addition to the generation of sinograms and appropriate data streaming processes such as attenuation, random and scatter corrections, a set of precursors to the reconstruction processes. To remedy the computational burden of image reconstruction, a number of alternative proposals such as linear integration have been proposed, as well as the use of multiple CPUs or a cluster computer system approach. Most of the techniques, however, are not practically useful since such cluster computing requires a large data transfer time, although the overall computation is faster than a single unit.
Obtaining or generating projection data can be divided into two categories, namely, raydriven method and pixeldriven method. The raydriven method calculates the linear integration directly along the ray path connecting the centers of the two opposite detector cells, whereas the pixeldriven method calculates the linear integration along the ray path centered around an image pixel for the entire projection angles. The raydriven method is often used in projection, while pixeldriven method is used in backprojection.
In early reconstruction techniques, projection was obtained by weighing the ray passing through the areas of pixels with the assumption that the ray path is a strip with a finite width. It, however, involves a large amount of computation as well as the storage of a large number of matrix or data. Concurrently, Shepp and Logan proposed a simple and computationally efficient algorithm, which requires computing the length of the ray path intersecting with each pixel (instead of the areas).
To speed up the computation, there have been a number of attempts to reduce the reconstruction time. An incremental algorithm has also been developed in which the symmetric property between the neighboring pixels is considered to calculate the position of intersection of a ray. This idea was expanded to 3D reconstruction in cylindrical geometry using oblique rays. In 3D form with a multiring system such as HRRT, it became apparent that true 3D approaches will be required to fully utilize the oblique rays to thereby improve the statistics of the image.
Arrangements and embodiments may be described in detail with reference to the following drawings in which like reference numerals refer to like elements and wherein:
A detailed description may be provided with reference to the accompanying drawings. One of ordinary skill in the art may realize that the following description is illustrative only and is not in any way limiting. Other embodiments of the present invention may readily suggest themselves to such skilled persons having the benefit of this disclosure.
Overview of the Projection, Backprojection and Symmetry Properties
A. Overview of Projection and Backprojection in the Aligned (Reference) Frame with Rotated Projection Plane
In 3D tomographic image processing, the projection ray has 4 dimensions, namely, x_{r}, y_{r}, φ and θ, as shown in the FIG. 1. Projection is a process of converting a 3D object function or image information in 3D coordinates to 2D projection coordinates, and the projection plane as shown in FIG. 1.
With a fixed view angle φ and oblique angle θ, the projection plane is determined as shown in FIG. 1A. In this projection plane, the definition of the projection operation can be expressed in the form of linear integral given by:
where
{right arrow over (ι)}_{x} _{ r } _{,y} _{ r } _{,φ,θ }is the ray path,
x_{r}=x′=x cos φ+y sin φ,
y_{r}=z−(−x sin φ+y cos φ)tan θ.
Equation (1) describes that projection P_{φ,θ}(x_{r},y_{r}) is a sum of the pixel values in an image function I(x,y,z) along the path of projection ray, {right arrow over (ι)}_{x} _{ r } _{,y} _{ r } _{,φ,θ}, in the image domain. Here, δ(.) represents a sampling function.
It is known that a projection plane at angle θ can be rotated (or aligned) against reference axes, (x,y). Further, x′ coincides with x_{r }in FIG. 2. The well known relation between the rotated coordinate (x′,y′,z) and the image coordinate (x,y,z) in the cylindrical coordinate are provided by the following:
where
R_{φ }is the rotation matrix.
Equation (1) can then be simplified as the weighed sum along the ray path and can be written as follows:
where
I(x′,y′,z)=R_{φ}I(x,y,z)
x′=x_{r }
In (3),y′ is assumed to be the integration variable.
As a special case, when the oblique angle θ is equal to zero, the projection rays are in parallel with the y′ axis (see FIG. 1). In multilayer or 3D PET, oblique rays (θ≠0) are collected and used for true 3D reconstruction since the full utilization of the oblique rays will enhance the image quality due to the increased statistics. Extension of (3) to a 3D case can be represented by θ≠0. It should be noted that the coordinates in the transverse plane are now independent of θ and the ray projected to the transverse plane is parallel to the y′ axis. By the trigonometric relationship depicted in
z=y_{r}+y′tan θ (4)
By substituting (4) into (3), the following discrete form of projection data can be obtained:
where
z_{y} _{ r } _{,n,θ}=y_{r}+n tan θ=Z_{y} _{ r } _{,n,θ}+r_{y} _{ r } _{,n,θ}
Z_{y} _{ r } _{,n,θ}εINTEGER, r_{y} _{ r } _{,n,θ}εREAL, and 0≦r_{y} _{ r } _{,n,θ}<1.
Z_{y} _{ r } _{,n,θ }is integer value of Z_{y} _{ r } _{,n,θ}
r_{y} _{ r } _{,n,θ }is interpolation coefficient or remainder of Z_{y} _{ r } _{,n,θ}
r_{y} _{ r } _{,n,θ}+r_{y} _{ r } _{,−n,θ}=1 or r_{y} _{ r } _{,−n,θ}=1−r_{y} _{ r } _{,n,θ}
n is discrete value of y′

 is integration length of y′ at given y′ at given y_{r}, φ, and θ.
It should be noted that since projection rays and the coordinates (x,y,z) or (x′,y′,z) are not always coincident. Interpolation is always required. In case of SSP implementation, a 1D linear interpolation is needed along the z axis, as shown in (5). Note that r_{y} _{ r } _{,n,θ }and 1−r_{y} _{ r } _{,n,θ }are the interpolation coefficients along the g axis.
On the other hand, the backprojection process is the transpose of projection. Therefore, the same idea as in the projection (or projection data generation) process can be applied to backprojection or image data generation. In case of an iterative reconstruction such as the OSEM algorithm, a plurality of projections and backprojections are required. That is, algorithm requires many repeated image reconstructions as well as projection data generations.
Now the backprojection from projection data, P_{φ,θ}(x_{r},y_{r}), can be performed, as shown by the following:
By using rotation, the matrix defined in (2), (6) can also be written as:
In the aligned frame, x_{r}=x′ is defined.
To simplify the notation, let us define I_{φ }as the sum of the projection planes about θ for a fixed φ as given below as an intermediate stage image before the finally reconstructed image. The computation can then be further simplified, for example, by using the symmetry property in θ,
Using the rotation matrix given in (2) once again, (8) can be noted as follows:
I_{φ}(x′,y′,z)=R_{φ}I_{100 }(x,y,z)
or I_{φ}(x,y,z)=R_{−φ}I_{φ}(x′,y′,z) (9)
Using (9), (7) can further be noted as:
Equation (10) is equivalent to rotating the set of partially reconstructed image planes at a given angle φ and then integrating these partially reconstructed images over all φ's to form the finally reconstructed image I(x,y,z).
B. Symmetry Properties
As discussed above, the full utilization of symmetry properties will be the central part of the proposed SSP method. Hereafter, the symmetric points represent the similar coefficients (the same value but with different polarity and/or different coordinates) that will be shared in computation without additional calculations as well as the use of identical instructions. That is, the symmetric points have the same interpolation coefficient or complementary value. The applicant found that there exists sixteen symmetric points, which are practically usable, namely “Mirrorsymmetry,” “φsymmetry,” “y′symmetry” and “θsymmetry.”
(a) Mirrorsymmetry (+x′ and −x′, see
One of the interesting aspects of the image reconstruction is the use of the symmetric properties of each image point. That is, once one point is calculated with an appropriate interpolation coefficient, one can assign the same coefficient value in the opposite point(s) by using symmetry properties of various types. One of the simplest forms is the use of the y′axis symmetry, as shown in FIG. 4A. In this computation, a point, (x′_{0},y′_{0}) at +x_{r }is identical in value (x′_{1},y′_{1}) at −x_{r }with change in polarity only. An illustration of this “Mirrorsymmetry” is shown in FIG. 4A. In this illustration, an example is given for symmetric points due to “Mirrorsymmetry”, i.e., (x′_{0},y′_{0})→(x′_{1},y′_{1})=(−x′_{0},y′_{0}). This symmetry property reduces the computation time by half.
(b) φsymmetry (φ and φ+90° see
Similar to the case of the “Mirrorsymmetry”, one can also demonstrate that there is symmetry between the points at φ and φ+90° symmetry. In this case, (x′_{0},y′_{0})→(x′_{1},y′_{1})=(−y′_{0}, x′_{0}), i.e., the polarity as well as the coordinates are changed.
(c) y′symmetry (+n and −n, see
Unlike mirror and φsymmetry, y′ and θsymmetry are in the y′z plane. It should be noted that y′symmetry is intimately related to θsymmetry. Therefore, it is easy to illustrate the two together, as follows. The y′symmetry property has the effect of changing the interpolation coefficients to complementary values. In fact, y′symmetry has the same symmetry property as θsymmetry, as discussed below. First, the property of y′symmetry is illustrated in FIG. 5A. As shown in
z_{y} _{ r } _{,n,θ}=y_{r}+ntanθ=Z_{y} _{ r } _{,n,θ}+r_{y} _{ r } _{,n,θ} (11)
where
n, Z_{y} _{ r } _{,n,θ and y} _{r}εINTEGER
r_{y} _{ r } _{,n,θ}εREAL and 0≦r_{y} _{ r } _{,n,θ}<1
where
r_{y} _{ r } _{,n,θ }and 1−r_{y} _{ r } _{,n,θ }represent die interpolation coefficients.
Using this symmetry, two interpolation operations can be simultaneously performed, thereby reducing the number of loops in y′(=n) by half. Therefore, (5) can be rewritten as follows:
Equation (13) can be transformed into a factored form to reduce the number of instructions required to execute the computation. Furthermore, the interpolation coefficient r_{y} _{ r } _{,n,θ }can be reduced to simple functions of only n and θ and can be derived from (11). Because y_{r }is always an integer and r_{y} _{ r } _{,n,θ }satisfies the condition 0≦r_{y} _{ r } _{,n,θ}<1, r_{y} _{ r } _{,n,θ }is a function of only n and θ. Therefore, the notation r_{y} _{ r } _{,n,θ }can be simplified to r_{n,θ }(see and compare with FIGS. 3 and 5A). Further, (1−r_{y} _{ r } _{,n,θ}) is equal to r_{y} _{ r } _{,−n,θ}. Similarly, (1−r_{y} _{ r } _{,−n,θ}) is equal to r_{n,θ}. Using the simplified forms, (13) can be transformed to a factored form, as shown below:
(d) θsymmetry (+θ and −θ, see
Now, let us illustrate how the θsymmetry shares the calculation of interpolation coefficients with the y′symmetry discussed above. As shown in
Using (4), (5), (11) and (12), it can be shown that z and r are given as:
Using the relation given in (15), (14) can be written as follows:
Equation (14) is for P_{φ,θ}(x_{r},y_{r}), while (16) is for the P_{φ,−θ }(x_{r},y_{r}). For the computation of (16), the relations given in (15), i.e., r_{n,θ}, r_{n,−θ}, Z_{y} _{ r } _{,n,θ}, and Z_{y} _{ r } _{,n,−θ }can be used once again and further simultaneously calculate P_{φ,θ}(x_{r},y_{r}) and P_{φ,−θ(x} _{r}, y_{r}).
As shown above, y′ and θsymmetries share the common properties. The y′symmetry, therefore, is used to reduce the number of loop counts y′(=n), while the θsymmetry is used to reduce the number of loop counts θ.
In summary, a total of sixteen points of symmetry for the SSP method, namely, 2 (Mirrorsymmetry), 2 (φ+90° symmetry), 2 (y′symmetry) and 2 (θsymmetry), are used. Thus, multiplication of these leads to 16. In addition to these symmetry based reductions in computation time (16 times), as will be detailed, SIMD will speedup the computation further by another factor of 4. As such, a total speed gain of 64 (16×4) has been achieved with the new SSP method.
Parallelization of Computation using SIMD for SSP method
A. SIMD (Single Instruction Multiple Data) and Symmetry Properties
It is important to reduce the number of instructions per loop for practical purposes, while reducing the number of loops within the projector/backprojector, to realize a fast algorithm. The number of instructions per loop was reduced by using the SIMD technique. The use of SIMD reduces the number of instructions per loop, thereby reducing the overall computation time by a factor of four. The operand of SIMD should be determined carefully to improve the efficiency. As shown above, projection/backprojection of sixteen points were simultaneously performed. For this operation, sixteen symmetry points were divided by two groups, one for +θ and the other for −θ. Since a single SIMD instruction data set consists of four data sets, each group having the same (+θ or −θ), was coupled to one of the two SIMD instruction datasets.
These two groups having the same θ is based on (14) and (16).
As noted above, sixteen symmetry pairs are divided into two groups, each sharing the same θ. In addition, these two groups were once again divided into two subgroups. The members of subgroups a1 and b2 share the same interpolation coefficient (e.g., r_{y′,θ}), while the subgroup a2 and b1 have complementary values (e.g., 1−r_{y′,θ}) as interpolation coefficients. Each subgroup represents one packed SIMD data, as shown below:
where
SI means an SIMD image data set.
SP means an SIMD projection data set.
Equation (17) shows the arrangement of the SIMD data set for the upper group of FIG. 12. This process is referred to as packing and single SIMD instruction data set consists of four points. This process can be applied to the lower groups having −θ.
As noted above, the upper group and the lower group of
where
SI_{x} _{ r }(n,Z_{y} _{ r } _{,n,θ}), SI_{x} _{ r }(n,Z_{y} _{ r } _{,n,θ}+1) and SP_{φ,θ(x} _{r},y_{r}) are according to subgroup “a1” of FIG. 12.
SI_{x} _{ r }(−n,Z_{y} _{ r } _{,−n,θ}), SI_{x} _{ r }(−n,Z_{y} _{ r },−nθ+1) and SP_{φ,θ}(x_{r},y_{r}) are according to subgroup “a2” of FIG. 12.
SI_{x} _{ r }(n,Z_{y} _{ r } _{,−n,θ}), SI_{x} _{ r }(n,Z_{Y} _{ r } _{,−n,θ}+1) and SP_{φ,−θ}(x_{r},y_{r}) are according to subgroup “b1” of FIG. 12.
SI_{x} _{ r }(−n,Z_{r} _{ r } _{,n,θ}), SI_{x} _{ r }(−n,Z_{y} _{ r } _{,n,θ}+1) and SP_{φ,−θ}(x_{r},y_{r}) are according to subgroup “b2” of FIG. 12.
Z_{y} _{ r } _{,n,θ}=Z_{y} _{ r } _{,−n,−θ}
Z_{y} _{ r } _{m,n,θ}=Z_{y} _{ r } _{,n,−θ}
This SIMD packing concept can be applied equally to a backward projection.
As a precaution, it should be noted that there are several special cases within the above symmetric pairs. These special cases indicate that some symmetry pairs will share the same rotating image point, thereby resulting in errors. TABLE II shows these special cases, which require special attention and correction.
TABLE II  
THE SPECIAL CASES WITHOUT SYMMETRY COUNTERPART  
Symmetry pairs on image plane  Special case  
φsymmetry  I(x′,y′,z_{y} _{ r } _{,y′,θ}) and I(−y′,x′,z_{y} _{ r } _{,y′,θ})  The points of x′ = y′ 
Mirror  I(x′,y′,z_{y} _{ r } _{,y′,θ}) and I(−x′,y′,z_{y} _{ r } _{,y′,θ})  The points of x′ = 0 
symmetry  
φsymmetry  I(x′,y′,z_{y} _{ r } _{,y′,θ}) and I(x′,=y′,z_{y} _{ r } _{,y′,θ})  The points of y′ = 0 
B. Optimization of Job Distribution in MultiProcessor Environment
The recently available PC provides a multiprocessor environment. To utilize this parallel computation capability, an instruction called “thread programming technique” can be used. The proposed SSP method is capable of being expanded to multiprocessor systems. In this case, a job distribution with an equal amount of load is important for optimal performance.
As the method according to the present invention is based on the assumption of circular cylindrical FOV, the job load is designed to be distributed equally along the x′ or x_{r }axis. This is because there is uneven distribution of computational load if an equal range of x′ is assigned to each CPU, as illustrated in FIG. 6.
Description of the Overall Method
The schematics of the new SSP projector and backprojector are provided in
As shown in FIG. 7 and
To reduce processing time, the efficient use of cache was found to be an important factor in the overall computation. Memory allocation and loop orders should be optimized for the memory access pattern. Since x_{r}(=x′) is set to the variable of the outermost loop, the memory stride on read and write can be confined to a certain range within the L2 cache size. Furthermore, the variable of the innermost loop should be z or y_{r}, since most of the computation intensive operations in the SSP method are related to interpolation of z or y_{r}. Therefore, it is desirable to allocate the initial memory by the order [x′]:[y′]:[z] for backprojection while for projection, the order should be [φ]:[x_{r}]:[θ]:[y_{r}]. These memory allocations, therefore, will determine the loop orders. These loop orders of projector/backprojector are also depicted in FIG. 7 and FIG. 8. As a result of using these symmetries, the total number of loop counts can be reduced by half at each step.
For the SIMD, a memory packing (step) is performed prior to each projection/backprojection step. After the projection step, the packed SIMD data set in the projection domain is unpacked and returned to the original position. Similarly, the packed SIMD data set in the image domain is unpacked after each backprojection step.
Experimental Results
A. Methods
To evaluate the efficiency of the SSP method, the data or sonogram was applied to HRRT (High Resolution Research Tomograph developed by CPS/Siemens, Knoxville, Tenn., U.S.A). HRRT has 119,808 detectors and is designed for ultra high resolution brain scanning. Since HRRT has the largest number of detectors (with the smallest detectors size) among the human PET scanners that are currently available, it has the highest computational complexity for image reconstruction. The number of computations has been increased substantially (in proportion to the number of LORs) and is much larger compared to other existing PET scanners. The reconstruction software provided by the HRRT package supports iterative algorithms such as OPOSEM3D and is the preferred reconstruction algorithm in the HRRT system. One of the advantages of OPOSEM3D, among others, is preventing the occurrence of negative bias when the random rate is high. This algorithm was chosen to test the performance of our new SSP method since it requires the largest memory size and is the most computationintensive algorithm. Therefore, unless otherwise specified, the HRRT OPOSEM3D package S/W will be referred to as the existing method.
As can be seen, the SSP method is supported by commercially available computer systems such as a common PC. Three other platforms were chosen to compare the SSP method against the existing method. The first platform is PC 1, which has Intel dualcore 3.0 GHz, 4 GB RAM, 1 MB L2 cache per CPU. The second platform, PC2, is a highend PC with two dualcore AMD 2.4 GHz, 8 GB RAM, 1 MB L2 cache per CPU. The third platform is the current HRRT computer platform, which consists of an eight node cluster system (this is denoted as HRRT CPS system), wherein each node consists of two Xeon 3.06 GHz processors with 512 KB L2 cache and 2 GB RAM.
The concept of “span” in 3D means the mode of axial compression. The test was performed in span 3 and span 9, which are typically used in HRRT. The parameters for OPOSEM3D are as follows: 256×256×207 image pixels; six iterations; and 16 subsets. Thread programming is also used.
B. Results
Before discussing the computational speed, the accuracy of the reconstructed images (compared to the original HRRT package provided by CPS/Siemens) appears to have been emphasized. First, the results of the projector and backprojector of the SSP were compared with the existing HRRT package. As shown in
The performance aspect of the SSP method (i.e., execution time on PC 1 and 2) was also compared. The computation time was measured for the two existing compression modes, namely, span 3 and span 9. The results indicate that the SSP method, as expected, indeed enhanced the performance by almost 80 times, 2 (φ+90° symmetry)×2 (x_{r }or mirrorsymmetry)×2 (y′symmetry)×2 (θsymmetry)×4 (SIMD)+other optimizations, compared to the existing method (original HRRT package). It should be noted that this performance enhancement can be obtained more or less independent of the system architecture. Furthermore, as the method suggests, the relative improvement will be higher as more oblique rays or angles are used, as seen from TABLES III and IV.
TABLE III  
THE COMPARISONS OF RUNNING TIME  
OF THE VARIOUS OPERATIONS ON PC1.  
Existing  new  ratio  
Projection (span 9)  1539(s)  24(s)  64.1  
Backprojection (span 9)  1527(s)  25(s)  61.1  
Projection (span 3)  4276(s)  51(s)  83.8  
Backprojection (span 3)  4255(s)  53(s)  80.3  
Projection + backprojection  3066(s)  49(s)  62.6  
(span 9)  
Projection + backprojection  8531(s)  104(s)  82.0  
(span 3)  
PC1: Intel 3.0 GHz 1 dualcore CPU, 4 GB RAM, 1 MB L2 cache per CPU. 
TABLE IV  
THE COMPARISONS OF RUNNING TIME  
OF THE VARIOUS OPERATIONS ON PC2  
Existing  new  ratio  
Projection (span 9)  556(s)  9(s)  61.8  
Backprojection (span 9)  679(s)  11(s)  61.7  
Projection (span 3)  1600(s)  22(s)  72.2  
Backprojection (span 3)  1767(s)  22(s)  80.3  
Projection + backprojection  1235(s)  20(s)  61.7  
(span 9)  
Projection + backprojection  3608(s)  44(s)  82  
(span 3)  
PC2: AMD 2.4 Ghz 2 dualcore CPU, 8 GB RAM, 1 MB L2 cache per CPU. 
The use of more oblique rays or angles (e.g., span 3) results in better axial resolution and statistics. Since reconstruction using the EM algorithm consists of a large number of projections and backprojections, this improvement in speed is especially important. Finally, a comparative study of the HRRT CPS system was compared with the original HRRT CPS algorithm vs. our 2 dual core CPU (without cluster) configuration, PC2, with the SSP method and obtained a speed gain of a factor of 9 to 11. The reason for the differences between the calculation (64 times) and the actual system based computation (911 times) is probably due to the 16 XeonCPUs based original HRRT CPS system. The HRRT CPS system having 16 XeonCPUs is much more powerful than PC2 having only 2 dualcore CPUs. This H/W difference could have been the result of such reduced speed gain. In addition, SSP method can easily extend to the cluster version.
The proposed SSP method will improve the overall computation time, especially when a dynamic functional study is desired. Currently, the image reconstruction time for a normal HRRT PET scan with a commercially available PC (PC2 type) requires about seven to eight minutes for span 3 after completing the sinogram formation and the precorrection processes compared with eighty minutes with current HRRT reconstruction package (OPOSEM3D) utilizing an 8node cluster system.
The present invention provides a fast method of calculating the ray path, which leads to the projection data and reconstructed image in 3D by using the symmetry properties of the projection and image data. Exploiting these simple geometrical symmetry properties and with the help of SIMD, a new ultra fast reconstruction method was developed, which has a computational speed advantage of nearly two orders of magnitudes compared to the existing reconstruction package (specifically that of OPOSEM3D, the latest and most widely used PET method, especially for HRRT without compromising image quality). The key concepts introduced in the SSP method can be summarized as follows. First, interpolation operations are reduced from three dimension to one dimension by using rotation based projection, or the aligned frame concept. Second, the rotation based projection is combined with symmetric properties of the (θ,φ,y′,x_{r}) and coupled with SIMD technique with optimized L1 and L2 cache use. The sixteen symmetry points, which share the same interpolation coefficients (or their complementary values), were grouped, thus permitting them to be processed simultaneously in a projection/backprojection operation to thereby achieve a substantially improved overall reconstruction time. To simultaneously process these symmetric points in the projection/backprojection, SIMD operations have also been incorporated. In particular, the SIMD scheme allowed us to simultaneously access four data. In addition, the data size per loop suitable for an L2 cache size was optimized, as well as the data structures for minimizing the memory stride.
In summary, the projection and backprojection operations were performed in the aligned frame by using the symmetry concept and SIMD with cache optimization and reduced the number of instructions in the reconstruction of an image from the sinogram data provided by the HRRT PET. In other words, the proposed SSP method incorporates the symmetry properties of projection to maximum, thereby reducing the computation time by as much as 16 times. Together with the use of the SIMD operator and cache optimization, an overall computational speed gain by a factor close to 80 was achieved. Such an improvement in computation time will open the door for the dynamic functional study of high resolution PET such as HRRT, which suffered from an unacceptably long reconstruction time and has been a serious undesired obstacle in molecular imaging.
This improvement will provide numerous new opportunities for PET users, especially for HRRT users in the molecular imaging community. It is clear that the proposed method helps the PET researchers to contemplate a possible dynamic study, which was limited by the unacceptable image reconstruction time imposed by the reconstruction process. This newly developed SSP method can also be applied to precorrection processes such as attenuation and scatter correction. Finally, with the new SSP method, it will be possible to achieve true interactive scanning. The new SSP method can also be applied to large classes of existing PET scanners of various types as well as to the cluster systems for dynamic studies.
Any reference in this specification to “one embodiment,” “an embodiment,” “example embodiment,” etc., means that a particular feature, structure or characteristic described in connection with the embodiment is included in at least one embodiment of the invention. The appearances of such phrases in various places in the specification are not necessarily all referring to the same embodiment. Further, when a particular feature, structure or characteristic is described in connection with any embodiment, it is submitted that it is within the purview of one skilled in the art to effect such feature, structure or characteristic in connection with other ones of the embodiments.
Although embodiments have been described with reference to a number of illustrative embodiments thereof, it should be understood that numerous other modifications and embodiments can be devised by those skilled in the art that will fall within the spirit and scope of the principles of this disclosure. More particularly, various variations and modifications are possible in the component parts and/or arrangements of the subject combination arrangement within the scope of the disclosure, the drawings and the appended claims. In addition to variations and modifications in the component parts and/or arrangements, alternative uses will also be apparent to those skilled in the art.
Claims (21)
Priority Applications (8)
Application Number  Priority Date  Filing Date  Title 

KR20060042155  20060510  
KR1020060042155  20060510  
KR20070027305A KR20070110186A (en)  20060510  20070320  Method and apparatus for ultra fast symmetry and simd based projectionbackprojection for 3d pet image reconstruction 
KR1020070027305  20070320  
KR20070040515A KR20070110191A (en)  20060510  20070425  Method and apparatus for ultra fast symmetry and simd based projectionbackprojection for 3d pet image reconstruction 
KR1020070040515  20070425  
US11800657 US7456406B2 (en)  20060510  20070507  Method and apparatus for ultra fast symmetry and SIMD based projectionbackprojection for 3D pet image reconstruction 
US12498952 USRE42034E1 (en)  20060510  20090707  Method and apparatus for ultra fast symmetry and SIMD based projectionback projection for 3D pet image reconstruction 
Applications Claiming Priority (1)
Application Number  Priority Date  Filing Date  Title 

US12498952 USRE42034E1 (en)  20060510  20090707  Method and apparatus for ultra fast symmetry and SIMD based projectionback projection for 3D pet image reconstruction 
Related Parent Applications (1)
Application Number  Title  Priority Date  Filing Date  

US11800657 Reissue US7456406B2 (en)  20060510  20070507  Method and apparatus for ultra fast symmetry and SIMD based projectionbackprojection for 3D pet image reconstruction 
Publications (1)
Publication Number  Publication Date 

USRE42034E1 true USRE42034E1 (en)  20110118 
Family
ID=38789018
Family Applications (2)
Application Number  Title  Priority Date  Filing Date 

US11800657 Active US7456406B2 (en)  20060510  20070507  Method and apparatus for ultra fast symmetry and SIMD based projectionbackprojection for 3D pet image reconstruction 
US12498952 Active USRE42034E1 (en)  20060510  20090707  Method and apparatus for ultra fast symmetry and SIMD based projectionback projection for 3D pet image reconstruction 
Family Applications Before (1)
Application Number  Title  Priority Date  Filing Date 

US11800657 Active US7456406B2 (en)  20060510  20070507  Method and apparatus for ultra fast symmetry and SIMD based projectionbackprojection for 3D pet image reconstruction 
Country Status (2)
Country  Link 

US (2)  US7456406B2 (en) 
DE (1)  DE102007020879A1 (en) 
Families Citing this family (5)
Publication number  Priority date  Publication date  Assignee  Title 

WO2011027402A1 (en) *  20090904  20110310  株式会社島津製作所  Nuclear medicine data processing method and nuclear medicine diagnosis device 
FR2989785B1 (en) *  20120420  20140613  Imacisio  Method for Characterization of a toe system, acquisition method, image reconstruction method, device and system toe implementing a such a process 
KR20140042461A (en)  20120928  20140407  삼성전자주식회사  Method and apparatus to correct motion 
CN103593868B (en) *  20131012  20160427  沈阳东软医疗系统有限公司  Based on the projection image reconstruction method and apparatus of the pet 
CN104484857B (en) *  20141226  20170818  国网重庆市电力公司电力科学研究院  One kind of meter data reading method and system 
Citations (8)
Publication number  Priority date  Publication date  Assignee  Title 

US6282257B1 (en) *  19990623  20010828  The Board Of Trustees Of The University Of Illinois  Fast hierarchical backprojection method for imaging 
US20020054659A1 (en) *  20000814  20020509  Miwa Okumura  Radiation detector, radiation detecting system and Xray CT apparatus 
US20020106051A1 (en)  20001220  20020808  Wido Menhardt  Image reconstruction using multiple Xray projections 
US20020122528A1 (en)  20001229  20020905  Besson Guy M.  Computed tomography (ct) weighting for high quality image reconstruction 
US20020154727A1 (en)  20010216  20021024  Ruola Ning  System and method for fast parallel conebeam reconstruction using one or more microprocessors 
US20030161443A1 (en) *  20020228  20030828  The Board Of Trustees Of The University Of Illinois  Methods and apparatus for fast divergent beam tomography 
US20050123089A1 (en)  20031209  20050609  Man Bruno D.  Method and apparatus for reduction of artifacts in computed tomography images 
US7522695B2 (en) *  20050719  20090421  Ge Medical Systems Global Technology Company, Llc  Xray CT apparatus 
Family Cites Families (3)
Publication number  Priority date  Publication date  Assignee  Title 

JP2006062002A (en)  20040825  20060309  Oki Electric Ind Co Ltd  Method of segmenting substrate of semiconductor device 
KR20070027305A (en)  20050906  20070309  아성기계 주식회사  Steam nozzle of pattern processing apparatus for various designshape of fiber fabrics 
KR20070040515A (en)  20051012  20070417  삼성전자주식회사  Method for manufacturing a specimen 
Patent Citations (8)
Publication number  Priority date  Publication date  Assignee  Title 

US6282257B1 (en) *  19990623  20010828  The Board Of Trustees Of The University Of Illinois  Fast hierarchical backprojection method for imaging 
US20020054659A1 (en) *  20000814  20020509  Miwa Okumura  Radiation detector, radiation detecting system and Xray CT apparatus 
US20020106051A1 (en)  20001220  20020808  Wido Menhardt  Image reconstruction using multiple Xray projections 
US20020122528A1 (en)  20001229  20020905  Besson Guy M.  Computed tomography (ct) weighting for high quality image reconstruction 
US20020154727A1 (en)  20010216  20021024  Ruola Ning  System and method for fast parallel conebeam reconstruction using one or more microprocessors 
US20030161443A1 (en) *  20020228  20030828  The Board Of Trustees Of The University Of Illinois  Methods and apparatus for fast divergent beam tomography 
US20050123089A1 (en)  20031209  20050609  Man Bruno D.  Method and apparatus for reduction of artifacts in computed tomography images 
US7522695B2 (en) *  20050719  20090421  Ge Medical Systems Global Technology Company, Llc  Xray CT apparatus 
NonPatent Citations (2)
Title 

I.K. Hong, et al., "Fast Forward Projection and Backward Projection Algorithm Using SIMD," IEEE, San Francisco 2006, Conference Program, M14372. 
Official Korean Action, with Translation, dated Sep. 18, 2008, as cited in Korean Patent Application No. 1020070040515. 8 Pgs. 
Also Published As
Publication number  Publication date  Type 

DE102007020879A1 (en)  20090402  application 
US7456406B2 (en)  20081125  grant 
US20070278410A1 (en)  20071206  application 
Similar Documents
Publication  Publication Date  Title 

DaubeWitherspoon et al.  An iterative image space reconstruction algorthm suitable for volume ECT  
Gullberg et al.  An attenuated projectorbackprojector for iterative SPECT reconstruction  
Vandenberghe et al.  Iterative reconstruction algorithms in nuclear medicine  
US6014419A (en)  CT cone beam scanner with fast and complete data acquistion and accurate and efficient regional reconstruction  
Guan et al.  A projection access order for speedy convergence of ART (algebraic reconstruction technique): a multilevel scheme for computed tomography  
US5881123A (en)  Simplified cone beam image reconstruction using 3D backprojection  
US5404293A (en)  Cone beam reconstruction using helical data collection paths  
US7394923B2 (en)  Imaging system for generating a substantially exact reconstruction of a region of interest  
DaubeWitherspoon et al.  Treatment of axial data in threedimensional PET  
US6292525B1 (en)  Use of Hilbert transforms to simplify image reconstruction in a spiral scan cone beam CT imaging system  
US6078638A (en)  Pixel grouping for filtering cone beam detector data during 3D image reconstruction  
US6477221B1 (en)  System and method for fast parallel conebeam reconstruction using one or more microprocessors  
US5907594A (en)  Reconstruction of volumetric images by successive approximation in conebeam computed tomography systems  
Reader et al.  Fast accurate iterative reconstruction for lowstatistics positron volume imaging  
Kinahan et al.  Analytic 3D image reconstruction using all detected events  
Fessler  Statistical image reconstruction methods for transmission tomography  
US6574299B1 (en)  Exact filtered back projection (FBP) algorithm for spiral computer tomography  
US6724856B2 (en)  Reprojection and backprojection methods and algorithms for implementation thereof  
Kachelriess et al.  Hyperfast parallel‐beam and cone‐beam backprojection using the cell general purpose hardware  
De Man et al.  Distancedriven projection and backprojection  
US6332035B1 (en)  Fast hierarchical reprojection algorithms for 3D radon transforms  
Mueller et al.  Fast implementations of algebraic methods for threedimensional reconstruction from conebeam data  
Herraiz et al.  FIRST: Fast iterative reconstruction software for (PET) tomography  
Reader et al.  EM algorithm system modeling by imagespace techniques for PET reconstruction  
US6018561A (en)  Mask boundary correction in a cone beam imaging system using simplified filtered backprojection image reconstruction 
Legal Events
Date  Code  Title  Description 

FPAY  Fee payment 
Year of fee payment: 4 

FPAY  Fee payment 
Year of fee payment: 8 