US RE42034 E1 Abstract A method and apparatus are provided for reconstructing 3D image. The method may include the steps of: detecting a plurality of line of response (LORs) emitted from an object; transforming the plurality of LORs into first sinogram data; back-projecting the first sinogram data with a plurality of projection angles to produce image data for the object; and projecting the produced image data with the plurality of projection angles to transform the image data into second sinogram data. The back-projecting may include filling pixels of image plane for each of the plurality of projection angles with the first sinogram data and rotating a coordinate axis of the image plane with corresponding projection angle to produce the image data. The projecting may include rotating the image data with corresponding projection angle in an opposite direction before projecting the image data with the plurality of projection angles. The projecting and the back-projecting may use symmetry properties in coordinate space.
Claims(21) 1. A method for reconstructing 3D image, comprising:
detecting a plurality of line of responses (LORs) emitted from an object;
transforming said plurality of LORs into first sinogram data;
back-projecting said first sinogram data with a plurality of projection angles to produce image data for said object; and
projecting said produced image data with said plurality of projection angles to transform said image data into second sinogram data,
wherein said back-projecting includes filling pixels of image plane for each of said plurality of projection angles with said first sinogram data and rotating a coordinate axis of said image plane with corresponding projection angle to produce said image data,
wherein said projecting includes rotating said image data with corresponding projection angle in an opposite direction before projecting said image data with said plurality of projection angles, and
wherein said projecting and said back-projecting use symmetry properties in coordinate space.
2. The method of
3. The method of
4. The method of
5. The method of
6. The method of
7. The method of
8. The method of
9. The method of
10. The method of
11. An apparatus for reconstructing 3D image, comprising:
means for detecting a plurality of line of responses (LORs) emitted from an object;
means for transforming said plurality of LORs into first sinogram data;
means for back-projecting said first sinogram data with a plurality of projection angles to produce image data for said object; and
means for projecting said produced image data with said plurality of projection angles to transform said image data into second sinogram data,
wherein said back-projecting includes filling pixels of image plane for each of said plurality of projection angles with said first sinogram data and rotating a coordinate axis of said image plane with corresponding projection angle to produce said image data,
wherein said projecting includes rotating said image data with corresponding projection angle in an opposite direction before projecting said image data with said plurality of projection angles, and
wherein said projecting and said back-projecting use symmetry properties in coordinate space.
12. The apparatus of
13. The apparatus of
14. The apparatus of
15. The apparatus of
16. The apparatus of
17. The apparatus of
18. The apparatus of
19. The apparatus of
20. The apparatus of
21. An apparatus for reconstructing a 3D image, the apparatus comprising:
apparatus for detecting a plurality of line of responses (LORs) emitted from an object; and a processor programmed with an algorithm for: transforming the plurality of LORs into first sinogram data; back-projecting the first sinogram data with a plurality of projection angles to produce image data for the object, the back-projecting including filling pixels of an image plane for each of said plurality of projection angles with said first sinogram data and rotating a coordinate axis of the image plane with a corresponding projection angle to produce the image data; and projecting the image data with the plurality of projection angles to transform the image data into second sinogram data, the projecting including rotating the image data with a corresponding projection angle in an opposite direction before projecting the image data with the plurality of projection angles, wherein the projecting and back-projecting use symmetry properties in coordinate space. Description The present application claims priority from Korean Patent Application No. 10-2006-0042155 filed on May 10, 2006, Korean Patent Application No. 10-2007-0027305 filed on Mar. 20, 2007, and Korean Patent Application No. 10-2007-0040515 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 projection-backprojection 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 back-projections. 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 back-projection (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, ray-driven method and pixel-driven method. The ray-driven method calculates the linear integration directly along the ray path connecting the centers of the two opposite detector cells, whereas the pixel-driven method calculates the linear integration along the ray path centered around an image pixel for the entire projection angles. The ray-driven method is often used in projection, while pixel-driven 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 multi-ring 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:
{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:
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:
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 multi-layer 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 By substituting (4) into (3), the following discrete form of projection data can be obtained:
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′
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 OS-EM 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:
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 “Mirror-symmetry,” “φ-symmetry,” “y′-symmetry” and “θ-symmetry.” (a) Mirror-symmetry (+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 “Mirror-symmetry” is shown in FIG. 4A. In this illustration, an example is given for symmetric points due to “Mirror-symmetry”, 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 “Mirror-symmetry”, 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 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 (Mirror-symmetry), 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 speed-up 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 data-sets. 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 a-1 and b-2 share the same interpolation coefficient (e.g., r_{y′,θ}), while the subgroup a-2 and b-1 have complementary values (e.g., 1−r_{y′,θ}) as interpolation coefficients. Each subgroup represents one packed SIMD data, as shown below:
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 “a-1” 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 “a-2” 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 “b-1” 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 “b-2” 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.
B. Optimization of Job Distribution in Multi-Processor Environment The recently available PC provides a multi-processor 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 multi-processor 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 OP-OSEM3D and is the preferred reconstruction algorithm in the HRRT system. One of the advantages of OP-OSEM3D, 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 computation-intensive algorithm. Therefore, unless otherwise specified, the HRRT OP-OSEM3D 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 dual-core 3.0 GHz, 4 GB RAM, 1 MB L2 cache per CPU. The second platform, PC2, is a high-end PC with two dual-core 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 OP-OSEM3D 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 mirror-symmetry)×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.
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 (9-11 times) is probably due to the 16 Xeon-CPUs based original HRRT CPS system. The HRRT CPS system having 16 Xeon-CPUs is much more powerful than PC2 having only 2 dual-core 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 (OP-OSEM3D) utilizing an 8-node 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 OP-OSEM3D, 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. Patent Citations
Non-Patent Citations
Classifications
Rotate |