KR20100109826A - Method and apparatus for ultra fast symmetry and simd based projection-backprojection for 3d pet image reconstruction - Google Patents

Abstract

PURPOSE: A symmetry and SIMD based high speed forward-backward projection method and a device thereof are provided to reduce the number of calculation times by using a symmetric property on a coordinate space which is obtained through the rotation of a coordinate axis. CONSTITUTION: A backward projection phase changes first sinogram data including data about a plurality of projection angles into image data. A forward projection phase changes the image data into second sinogram data. The backward projection phase comprises a step of generating first intermediate image data, second intermediate image data, and image data. The forward projection phase comprises a step of generating forward projection intermediate image data and second sinogram data.

Description

METHOD AND APPARATUS FOR ULTRA FAST SYMMETRY AND SIMD BASED PROJECTION-BACKPROJECTION FOR 3D PET IMAGE RECONSTRUCTION}

BACKGROUND OF THE INVENTION 1. Field of the Invention [0001] The present invention relates to a tomography apparatus for reconstructing and providing a three-dimensional image of a body or a living tissue in real time. More particularly, the present invention relates to an iteration scheme used in the process of reconstructing a three-dimensional image. More particularly, the present invention relates to a method and apparatus for optimizing the iterative operation of forward projection and back-projection using symmetry in projection or image coordinate space.

Significant advances have recently been made in the development of Positron Emission Tomograph (PET) in the computer implementation of hardware, software and image reconstruction. The High Resolution Research Tomograph (HRRT) developed by CTI (now Siemens) is a representative example of recent developments, and has a very improved performance in addition to sensitivity. For these modern PET scanners, the number of coincidence lines of response generated by as many as 120,000 nuclear detectors is about 4.5 x 10 ⁹ . Has This vast amount of data and its reconstruction, for image reconstruction and processing of related image data, poses a real problem for HRRT and has become a major neck for subsequent development and application of high resolution PET scanners. . Thus, for this kind of PET scanner, it often requires several hours of image reconstruction time to obtain a set of reconstructed images. For example, for a typical brain scan (using Span 3), for an HRRT with an entire data population, the image reconstruction time is nearly 80 minutes, which takes several days (32 frames of 43 hours or more on dynamic image reconstruction) makes it virtually impossible to attempt list mode based dynamic imaging.

In general, tomographic images can be reconstructed by two techniques, an analytical method and an iterative approach. In the mid-1970s, various types of PET scanners were developed, resulting in various types of tomography image reconstruction techniques. Analytical methods, such as backprojecton and filtered or filtered backprojection (FB), particularly in the case of the High Resolution Research Tomograph (HRRT), where detectors are placed in a set of blocks in octagonal form. If the detector arrangement is not uniform, artificial distortion, also known as streak artifact, is often generated. PET scanners using such block detectors often cannot avoid data loss due to gaps between the blocks, resulting in severe linear artificial distortion in the FB technique. Because of this, iterative methods such as EM (Expectation Maximum) algorithm, rather than analytical methods, have recently been used in image reconstruction. In general, the EM technique requires several steps in the reconstruction process, the largest of which in terms of computational amount is a forward projection for generating projection data from an image or a scan object. FP) and a backprojection (BP) step for reconstructing the final image data on the image domain using the projection data. In the EM technique, these two processes are repeated many times until a satisfactory image is obtained. Due to this repeated FP and BP process, the EM technique has a large amount of calculation compared to the straight filtered backprojection method, which is representative of the analytical method, and this calculation has been pointed out as a major drawback of the EM technique. Moreover, in three-dimensional image reconstruction, the amount of computation increases even more due to the huge increase in the number of coincidence lines or lines of responses (LOR). This indicates that the use of EM techniques in high resolution PET scanners with a larger number of lines of responses (LOR) is inevitable. As a result, in order to use EM techniques in ultra-high resolution PET scanners such as HRRT, there is a need to dramatically increase the computation speed.

The projection method usually uses a system matrix determined by the scanner's geometric factors. As the resolution of PET images increases and the number of slices increases, the size of the matrix also increases significantly in proportion to the increase in the number of LORs, which not only requires a large amount of memory, but also increases the overall computation time. Current HRRTs are, for example, a series of precursors to the reconstruction process, in addition to appropriate data streaming processes such as generation of a sinogram, attenuation, randomness and scatter correction. This requires approximately 80 minutes of reconstruction time. In order to reduce the computation time required for image reconstruction, cluster computer system techniques have been proposed and used in practice. The overall computation time is somewhat reduced by using cluster systems, but there are problems with cluster systems as well. As data grows, compute nodes on the cluster system must increase, and data communication between these compute nodes increases, which is a problem that this data communication is a very slow process. In other words, if there is an enormous increase in the data to be processed, such as an ultra-high resolution 3D PET scanner, the cluster system also requires significant computational time.

Most EM techniques based on conventional projection computations are extremely large and have inefficient computational complexity and require extremely long reconstruction times. In addition, even if the parallel processing technique using a cluster system using dozens of CPUs is adopted, it is faster than it is not, but it does not bring a satisfactory improvement. As mentioned earlier, in the case of the parallel processing technique, the data transfer between the compute nodes is increased, thereby substantially reducing the overall gain in parallel processing. In the end, we need an algorithm that can reduce parallelism and make parallelism more efficient.

An object of the present invention is to provide a three-dimensional image reconstruction method capable of real-time processing by reducing the number of operations by using the symmetry in the coordinate space obtained through the rotation of the coordinate axis in the iterative operation of the forward and backward projection method will be.

It is another object of the present invention to provide a three-dimensional image reconstruction method optimized for a single instruction plural data method (SIMD) which executes operations on a plurality of different data simultaneously in one calculation process.

Still another object of the present invention is to uniformly divide the width of the domain region to be computed, so that the plurality of divided regions are most efficiently performed in parallel by each of the plurality of processors, thereby reducing the execution time as much as possible. It is to provide a reconstruction method.

According to one aspect of the present invention for achieving the above object, there is provided a method for reconstruction of a three-dimensional image, the method comprising the steps of: detecting a plurality of reaction lines emitted from a predetermined object, a plurality of Converting a response line into first sinogram data, a back-projection step of projecting the first cosine curve data at a plurality of projection angles to generate image data for the predetermined object, and the And a forward projection step of projecting the generated image data to the plurality of projection angles and converting the generated image data into second cosine curve data, wherein the backward projection step includes the first cosine curve data for each of the plurality of projection angles. Filling the pixels of the image plane, rotating the coordinate axis of the image plane by a corresponding projection angle, and generating image data; The step includes rotating the image data in the opposite direction by the projection angle before projecting at the plurality of projection angles, and the forward projection and the backward projection steps use symmetry in coordinate space.

According to one embodiment of the invention, the symmetry is mirror symmetry (± x'symmetry) and symmetry in the viewing angle ( phi -symmetry), symmetry in the inclination angle (θ-symmetry) and interpolation coefficient along the projection line. At least one of symmetry (y'-symmetry) in.

According to another embodiment of the present invention, the forward projection step and the backward projection step include executing a plurality of symmetrical projection operations in parallel using a single instruction multiple data method (SIMD).

According to another embodiment of the present invention, the backward projection step and the forward projection step divide the image data into a plurality of even subsets, and perform projection operations on each of the plurality of subsets by each of the plurality of processors. .

According to another aspect of the present invention, there is provided an apparatus for reconstructing a three-dimensional image, the apparatus comprising: means for detecting a plurality of reaction lines emitted from a predetermined object, the plurality of reaction lines using a first cosine curve ( means for converting into sinogram data, means for performing a back-projection of projecting first cosine curve data at a plurality of projection angles to generate image data for the predetermined object, and the generated Means for performing a forward projection that projects the image data at a plurality of projection angles and converts the second cosine curve data into the second cosine curve data, wherein the backward projection comprises the first cosine curve data at an image plane for each of the plurality of projection angles. Filling a pixel of the image plane and rotating the coordinate axis of the image plane by the corresponding projection angle to generate the image data. Before the projection, and the angle of attack comprises rotated in the direction opposite to each projection by the image data, the forward and reverse projection projecting utilizes the symmetry on the coordinate space.

BRIEF DESCRIPTION OF THE DRAWINGS The above and other objects, advantages and features of the present invention will become apparent from the following detailed description, which is disclosed with reference to the accompanying drawings which illustrate certain embodiments of the present invention.

Overview of Forward / Reverse Projection and Symmetry Characteristics

A. Overview of forward and reverse projection in alignment (reference) frames with a rotating projection plane

와의 관계에 주목하라. 각도 +θ는 이미지 평면에서 상위측으로의 경사각을 가리킨다. x_r과 y_r는 3D 객체의 2D 투영 평면의 좌표이다. 도 1c는 θ= 0°에서의 수평면과 x_r 상에 투영된 y′선을 따르는 선형 적분의 예를 도시한다. 여기서, φ는 좌표축 (x,y)를 기준으로 한 투영 광선 세트의 회전을 나타낸다. 점선 화살표는 투영 광선을 나타낸다. 각 투영 광선은 네 개의 변수, 즉 (x_r,y_r,φ,θ)로 결정된다. 도 1a에 도시된 것처럼 투영 평면은 θ와 Φ에 의해 결정되며 좌표축 (x_r,y_r)의 2D 투영 데이터(또는 bed of nails)로 구성된다.FIG. 1A shows the projection of a 3D object and the object onto a 2D projection plane, and FIG. 1B is a y'-z top view of FIG. 1A. Here, the path of the projection light beam to be performed with the projection planes and linear integration, i.e.

Note the relationship with The angle + θ indicates the inclination angle from the image plane to the upper side. x _r and y _r are the coordinates of the 2D projection plane of the 3D object. 1C shows an example of linear integration along the y 'line projected on the horizontal plane and x _r at θ = 0 °. Here, φ represents the rotation of the set of projection light beams with respect to the coordinate axis (x, y). Dotted arrows indicate projection rays. Each projection ray is determined by four variables: (x _r , y _r , φ, θ). As shown in FIG. 1A, the projection plane is determined by θ and Φ and consists of 2D projection data (or bed of nails) of the coordinate axes (x _r , y _r ).

As shown in Fig. 1, in 3D tomographic image processing, the projection light ray has four dimensions (i.e., x _r , y _r , φ and θ ). Projection (or forward projection) is the process of converting 3D object functions or image information that is 3D coordinates into a projection plane that is 2D projection coordinates. The projection plane is as shown in FIG. Given a fixed viewing angle φ and tilt angle θ , the projection plane is determined as shown in FIG. 1A. In this projection plane, the definition of the forward projection operation can be expressed in the form of linear integration as in equation (1).

here,

는 광선경로이고,

Is the ray path,

이 이미지 도메인에서 투영 광선 경로

를 따르는 이미지 함수 I(x, y, z) 에서의 화소 값들의 합임을 나타낸다.

는 샘플링 함수이다.Equation 1 is a projection

The projection ray path in this image domain

Denotes the sum of pixel values in the image function I ( x, y, z ).

Is a sampling function.

Fig. 2 is an alignment in which the position of the projection beam and the projection plane is fixed while having the fixed image plane on (x, y) and having the image plane rotating as opposed to the concept of the projection beam rotation frame in which the projection beam (ray path) rotates. The concept of the projection frame is shown. 2A shows that the projection light beam (ray path) and the projection plane are rotated over the image plane on a fixed reference coordinate (x, y). This is a common technique applied to most image reconstructions. On the contrary, the alignment projection frame shown in FIG. 2B is an image plane rotated instead of the projection ray and the projection plane, and the coordinate system (x ', y') in which the x and y axes of the image coordinate system (x, y, z) are rotated by the corresponding projection angle. image space with z, and projection plane with cosine curve coordinate system (x _r , y _{r ,} φ, θ ). The technique using the alignment projection frame illustrated in FIG. 2B is called rotation based projection, and is the basis of the SSP (Symmetry and SIMD based Projection-backprojection) method of the present invention. This technique makes it very easy to use the geometric symmetry that exists on the image plane by rotating the image plane, and the use of the symmetry feature reduces the overall forward and backward projection time. Note that for backward projection, the image function I ( x ' , y' , z) represents the intermediate stage of the image to be reconstructed and is temporary image data.

The projection plane at angle φ can be rotated (or aligned) about the reference axis (x, y). Note that in Figure 2 x 'matches x _r . The relationship between image coordinates (x, y, z) and rotated coordinates (x ', y', z) in a cylindrical coordinate system is given by Equation 2.

는 회전 매트릭스.here,

Is the rotation matrix.

Equation 1 may simplify the sum of the image pixels along the ray path as a form of weight sum as shown in Equation 3.

In Equation 3, y 'is used as an integral variable.

As described above, the alignment projection frame shown in FIG. 2B used for rotation-based projection has a coordinate system (x ', y', which rotates the x and y axes of the image coordinate system (x, y, z) by a corresponding projection angle. It is configured to include an image space having z) and a projection plane having a cosine curve coordinate system (x _r , y _{r ,} Φ, θ), where the x _r axis of the projection plane is the x of the coordinate system (x ', y', z). The axis 'parallel to the axis and the y' axis of the coordinate system (x ', y', z) of the alignment frame is configured to be parallel to the projection rays.

As a special case, when the inclination angle θ is 0, the projection light beam becomes parallel to the y 'axis (see Fig. 1). In 2D PET, which uses only 0 degrees of θ in PET consisting of multiple planes, more data is available for image reconstruction for 3D PET using multiple oblique rays (if θ is not zero). Through this, the reconstructed image quality can be improved. For this 3D PET, Equation 3 can be extended to the 3D case (that is, θ is not 0), in which case the coordinates in the horizontal plane are independent of θ and for the horizontal plane of the projection beam. Remember that the projectiles are parallel to the y 'axis. According to the triangular relationship shown in FIG. 3, z may be described as in Equation 4 below.

및

사이의 관계를 도시한다. 이러한 관계는 모든 y'-z 평면 상의 어떤 x'에 대하여도 유효하다는 점에 주목하여야 한다. 또한 y'는 이산값 n으로 표시된다는 점도 유의해야 한다.3 illustrates z, y ', x', θ, x _r , y _r , using a y'-z plan view.

And

Shows the relationship between. It should be noted that this relationship is valid for any x 'on all y'-z planes. It should also be noted that y 'is represented by a discrete value n.

When equation (4) is applied to equation (3), discrete projection data, such as equation (5) below, can be obtained.

here,

및

And

는

의 정수값

Is

Integer value of

는

의 나머지 또는 보간계수

Is

Remainder or interpolation factor

n is the discrete value of y '

은 주어진 y _r ,Φ 및 θ에 대한 y' 의 적분 길이.

Is the integral length of y ' for a given y _r , Φ and θ.

혹은 그의 보수

를 사용한다.Interpolation is necessary because the coordinates of the projection pixels and the image pixels of actual values in the (x, y, z) or (x ', y', z) coordinate system do not always coincide. Unlike conventional methods that normally require 3D interpolation or 2D interpolation, the SSP only needs 1D linear interpolation along the z axis as shown in Equation 5. One-dimensional linear interpolation is the interpolation coefficient along the z axis, as shown in Equation 5.

Or his reward

Use

On the other hand, the backward projection process is to reconstruct a set of images from the sinogram generated by the projection data described above. Thus, the same manner as in the projection (or forward projection or generation of projection data) process can be applied to backward projection or image reconstruction. In the case of iterative reconstruction, such as the OS-EM algorithm, since a number of forward and backward projections are required, the algorithm requires not only projection data generation but also multiple repeated image reconstructions. Projection Data P _{Φ, θ} ( x _r , y _r ) may be performed as given in Equation 6 below.

Using the rotation matrix defined in Equation 2, Equation 6 may be described as Equation 7 below.

Note that in the alignment frame, x _r = x .

To simplify the notation, define the sum of projection planes for all θ at a particular Φ as I _Φ as an intermediate image before the final reconstructed image, as given in Equation 8 below. Then, for example, the symmetry property of θ can be used to further simplify the calculation.

Using the rotation matrix given in Equation 2 again, Equation 8 can be written as Equation 9.

Using Equation 9, Equation 7 may be expressed as Equation 10 below.

들을 각각의 Φ에 해당하는 -Φ의 회전변환 후 이들 각각을 모든 Φ에 대해 적분하여 최종 재구성 이미지 I(x,y,z)를 생성하는 것을 나타낸다.Equation 10 is an image reconstructed as an intermediate step that is rotated at a given angle Φ .

Corresponding to each of the Φ-integrated with respect to each of which after rotational transformation of the Φ Φ to indicate that all generate the final reconstructed image I (x, y, z) .

B. Symmetry Properties

The SSP method according to an embodiment of the present invention utilizes a symmetric property. The symmetric points in the interpolation operation herein hereinafter represent similar coefficients (having the same value, but having different polarities and / or different coordinates) that are commonly available without additional calculations as well as the use of the same instruction. That is, the symmetric points have the same interpolation coefficient or complementary value. Applicant has found that there are actually four symmetry relationships available, i.e. 16 symmetry points, each of which is "mirror symmetry", " φ -symmetry", "y'-symmetry" and " θ- Symmetrical ".

(a) mirror symmetry (+ x '  And- x ' Wow y ' shaft; See FIG. 4A)

A feature of the above invention in image reconstruction is to use the symmetric properties of each image point. Once a point has been calculated using the appropriate interpolation coefficients, various types of symmetry attributes can be used to reuse the same coefficient value at opposing point (s). One of the simplest forms is to use mirror image symmetry as shown in FIG. 4A. In this case, the point (x ' ₀ , y' ₀ ) at + x _r is the same as the value (x ' ₁ , y' ₁ ) at -x _r with only the polarity being different. This “mirror symmetry” is illustrated in FIG. 4A. 4A shows exemplary symmetry points due to “mirror symmetry” (ie, (x ' ₀ , y' ₀ ) → (x ' ₁ , y' ₁ ) = (-x ' ₀ , y' ₀ )). It is. This symmetry reduces the computation time by half.

(b) φ -Symmetry ( φ  And φ + 90 °; 4b)

As in the case of "mirror symmetry", there is symmetry between the points at φ and φ + 90 °. In this case, since (x ' ₀ , y' ₀ ) → (x ' ₁ , y' ₁ ) = (-y ' ₀ , x' ₀ ), the polarities as well as the coordinates are changed. 4B shows a numerical example of such a case of φ -symmetry. As shown in FIG. 4B, once the point is calculated at φ and the interpolation coefficients in the forward projection (or backward projection) calculation, for example for (x ' ₀ , y' ₀ ) are obtained, it is rotated 90 ° (ie The same coefficient values can be reused in the forward projection (or backward projection) for the point φ + 90 °, that is, (-y ' ₀ , x' ₀ ). This symmetry property cuts the computation time back in half.

4 shows examples of "mirror symmetry" and "φ-symmetry" used in the proposed SSP method. Figure 4a shows an example of "mirror symmetry" for y ', knowing the coordinate points (x' ₀ , y ' ₀ ), the symmetry points (x' ₁ , y ' ₁ ) for y' without additional calculation This can double or increase the computational speed or halve the computation time Figure 4b shows the symmetry characteristic of the projection data with respect to Φ , i.e. coordinate values (x ' ₀ , y'). _If the interpolation coefficient of ₀ ) is known at Φ , the interpolation coefficient at the symmetry point (x ' ₁ , y' ₁ ) at Φ + 90 ° can also be obtained without further calculation, ie (x, y) at Φ a coordinate transformation with a simple change to the (-y, x) in the Φ + 90 °, this has the effect of replacing the interpolation coefficient operation of the forward projection (or reverse projection) corresponding to the Φ + 90 °. this is the same symmetry characteristics Cut the calculation time in half.

(c) y ' Symmetry (+ n and -n; see FIG. 5A)

Unlike "mirror symmetry" and φ -symmetry, y ' -symmetry and θ -symmetry are symmetry in the y' - z plane. Note that y' -symmetry is closely related to θ -symmetry. Therefore, it is convenient to show the two together as follows. y' -symmetry has the effect of changing the interpolation coefficient to a complementary value. y' -symmetry has the same symmetry characteristic as θ -symmetry, as described below. First, the y' -symmetry characteristic is shown in FIG. 5A. As shown in FIG. 5A, the variable z is represented by Equation 11 below (see FIG. 3).

here,

및 y _r ∈ INTEGER

And y _r ∈ INTEGER

∈ REAL 및 0 ≤

< 1

AL REAL and 0 ≤

<1

here,

및

은 보간 계수.

And

Is the interpolation factor.

Using this symmetry, two interpolation operations can be performed simultaneously, thereby reducing the number of loops in y ' (= n) by half. Therefore, Equation 5 may be written as Equation 13 below.

는 n과 θ만의 단순한 함수로 줄일 수 있으며, 이는 수학식 11로부터 도출될 수 있다. y _r 은 항상 정수이고

는 0 ≤

< 1 의 조건을 만족하므로,

는 n과 θ만의 함수이다. 따라서,

라는 표기는 r_{n , θ} 로 단순화될 수 있다(도 3 및 도 5a 참조 및 비교). 또한, (1-

)는

와 동일하다. 마찬가지로, (1-

) 는 r_{n ,θ} 와 동일하다. 이를 수학식 13에 적용시 수학식 13은 이하의 수학식 14와 같이 분해된 형태로 변환할 수 있으며 이 경우 필요한 연산의 수가 줄어들게 된다.Where interpolation factor

Can be reduced to a simple function of only n and θ , which can be derived from Equation (11). y _r is always an integer

Is 0 ≤

Since the condition of <1 is satisfied,

Is a function of only n and θ . therefore,

May be simplified to r _{n , θ} (see and compare FIGS. 3 and 5A). In addition, (1-

)

. Similarly, (1-

) Is r _{n , θ} . When this is applied to Equation 13, Equation 13 may be converted into an exploded form as shown in Equation 14, in which case, the number of required operations is reduced.

(d) θ Symmetry (+ θ Wow - θ ; See FIG. 5B)

We now look at how θ -symmetry shares the calculation of interpolation coefficients with y' -symmetry described above. As shown in Fig. 5B, when the sign of θ changes, the difference in the interpolation coefficient is merely changed by the complementary value.

Using equations (4), (5), (11) and (12), z and r can be expressed as in the following equation (15).

therefore,

Using the relation given in Equation 15, Equation 14 may be expressed as Equation 16 below.

및

를 다시 이용할 수 있으며, P _Φ,θ (x _r ,y _r )와 P _Φ,-θ (x _r ,y _r )를 동시에 계산할 수 있다.Equation 14 is for P _{Φ, θ} ( x _r , y _r ), but Equation 16 is for P _{Φ, -θ} ( x _r , y _r ). For the calculation of equation (16), the relationship given in equation (15), i.e. r _{n , θ} , r _{n , -θ} ,

And

Can be used again, and P _{Φ, θ} ( x _r , y _r ) and P _{Φ, -θ} ( x _r , y _r ) can be calculated simultaneously.

As above, y' -symmetry and θ -symmetry share common characteristics. 5c briefly illustrates y' -symmetry and θ -symmetry. y ' -Symmetry is loop count θ -symmetry is used to reduce the number of loop counts θ , while reducing the number of y ' (= n).

In summary, the SSP method uses a total of 16 symmetry points (or pairs). In other words, when using mirror symmetry, Φ + 90 ° -symmetry, y' -symmetry, and θ -symmetry, two forward projections (or backward projections) are possible in one operation for each symmetry. Forward (or backward) projections with two symmetries are computed simultaneously. In addition to this symmetry-based computation time reduction (16 times), the SIMD as detailed below further speeds up the computation by another 4 times. This results in a total speed gain of 64 (16x4) for the new SSP method.

5 shows the y'-symmetry and θ-symmetry of the proposed SSP method. 5A shows y'-symmetry and FIG. 5B shows θ-symmetry. 5C shows the relationship between y'-symmetry and θ-symmetry or the overall symmetry obtainable by y'-symmetry and θ-symmetry. The left side of each figure shows the z-y 'plane cut-view, and the right side enlarges the left side. Each quadrant is a detailed illustration of the interpolation coefficient r. This y'-symmetry or θ-symmetry characteristic indicates that once one of the interpolation coefficients is calculated for one point, the remaining interpolation coefficients for the other three points can be obtained without further calculation. Thus, by using y'-symmetry and θ-symmetry, the amount of calculation can be reduced to 1/4.

SSP  Way of SIMD Parallel computation

A. SIMD ( Single Instruction Multiple Data ) And symmetrical characteristics

In order to realize a faster forward-backward projection method, it is also important to reduce the number of loops within the forward-backward projection while reducing the number of instructions per unit loop actually called. According to an embodiment of the present invention, the number of instructions per unit loop can be reduced by using a SIMD technique. The use of SIMD reduces the number of instructions per unit loop, thus ideally reducing the overall computation time to 1/4. Operands of the SIMD must be carefully determined to increase efficiency. As described above, according to one embodiment of the present invention, forward / backward projection can be performed simultaneously at 16 points. In this process, 16 symmetric points are divided into two groups, one divided by + θ and the other divided by- θ . The data set of a single SIMD instruction currently supported on a commercial CPU consists of four data. That is, one SIMD instruction can access and process four data at the same time. The set of four data that can be processed at the same time is called a SIMD data set. The task of grouping them is called packing. The task of unpacking them in one set after an operation is called unpacking. That is, each group (+ θ or − θ ) having the same θ has a total of eight forward (or backward) pairs, and thus each consists of two SIMD instruction data sets.

These two groups with the same θ are based on equations (14) and (16). Table 1 shows sixteen pairs of symmetrical relationships between the rotated image data set and the projection data set. As can be seen from the contents described in relation to the above four symmetry properties and Table 1, the symmetric points (i.e., symmetry pairs) that form a symmetry relationship on an alignment frame are coordinates and projection data on an image coordinate system that is image data. (I.e., cosine curve data), the coordinates on the projection coordinate system (i.e., the cosine curve coordinate system).

As shown in Table 1, the sixteen pairs of symmetry determined through the mirror symmetry relationship, φ-symmetry relationship, θ-symmetry relationship and y'-symmetry relationship are divided into two groups each sharing the same θ . In addition, these two groups are further divided into two subgroups. Members of subgroups a-1 and b-2 share exactly the same interpolation coefficients (eg r _{y ' , θ} ), and subgroups a-2 and b-1 complement (eg 1-r _{y' , θ)} ) As the interpolation factor. Each subgroup represents one packed SIMD data as shown in Equation 17.

Where SI stands for SIMD image data set and SP stands for SIMD projection data set.

Equation 17 describes the construction of a SIMD data set for the upper group of Table 1. It suggests how to do packing in SIMD. The packing scheme presented in this formula is centered on θ , but can be equally applied to two subgroups with -θ (one group is packed into one SIMD data set).

As described above, the upper and lower groups in Table 1 are based on Equations (14) and (16). Equations 14 and 16 may be modified to the SIMD version disclosed in Equation 17. That is, this is as shown in Equation 18 below.

here,

및

는 표 1의 하위 그룹 "a-1" 에 따른 것

And

Is according to subgroup "a-1" in Table 1

및

는 표 1의 하위 그룹 "a-2"에 따른 것

And

Will be according to subgroup "a-2" in Table 1

및

는 표 1의 하위 그룹 "b-1"에 따른 것

And

Will be according to subgroup "b-1" in Table 1

및

는 표 1의 하위 그룹 "b-2"에 따른 것

And

Will be according to the subgroup "b-2" in Table 1

This SIMD packing concept applies equally to backward projection.

It should be noted that there are some special cases within this symmetric pair. This special case means that in each symmetry some forward projection (or backward projection) is the same image point (where the image point is the point of the rotated image). Note that x ' , y' , and z coordinates are shared), resulting in overlapping errors. Table 2 discloses special cases requiring special care and modification. As shown in Table 2, if the point corresponding to the symmetry is itself, care must be taken not to duplicate the point in projections.

B. Optimizing Job Distribution in a Multiprocessor Environment

Modern commercial PCs provide a multiprocessor environment. To take advantage of this parallel computing power, instructions called "thread programming techniques" can be used. The proposed SSP method can be extended to a multiprocessor system. In this case, it is important to evenly distribute the work for optimal performance.

Since the method according to an embodiment of the present invention is based on the assumption of a circular or cylindrical FOV, the same x 'is shown in FIG. 6. If the range is assigned to each CPU, then the distribution of computation is uneven, so distribute the workload equally along the x ' , ie, the x _r axis.

6 illustrates the concept of uniform work distribution based on the sum of the beam path lengths. The amount of computation is equally distributed for each group S1, S2, S3 and S4 on x ' . In the case shown in Figure 6, the workload was divided into four identical groups of S1, S2, S3 and S4. This redistribution shows that each group has almost the same area, that is, the same amount of computation.

Technology of the whole way

A schematic of the new SSP forward and backward projections is shown in FIGS. 7 and 8, respectively. In the case of an SSP, according to an embodiment of the present invention, in order to reduce the number of loop counts, that is, to simultaneously process several data simultaneously through one calculation loop, on an aligned frame and an aligned frame through a rotation operation. Use several geometric symmetry relationships that are easier to use. In addition, by combining symmetry features and SIMD techniques, the total number of instructions required for computation can be reduced.

7 is a flowchart of a forward projection according to an embodiment of the present invention. The innermost block 706 in the loop includes z loops, interpolation operations, integration 708 along the ray path and scanning 710 along the x 'axis. By performing the packing process, projection operation using the SIMD instruction is possible. However, after projection over all ranges of θ, the intermediate result (x ', y', z coordinate system) of the rotated projection must be unpacked since it must be stored back in the actual image coordinate system (x, y, z). To maximize cache hit rate and data input / output speed, the size of the I-pack and SIMD data sets (I-Pack + P-Pack) are designed to be smaller than the cache sizes of L1 and L2, respectively.

8 is a flow chart of backward projection in accordance with an embodiment of the present invention. First, eight projection data having a symmetric relationship are packed for use in SIMD. This packed projection data (P-Pack) is linearly interpolated along the y _r direction at a particular θ to produce a packed backward projection image necessary to produce intermediate resultant image data I _Φ ( x ', y ', z ). Create data (I-Pack) After this I-Pack has been processed for all θ, the unpacking process allocates the data of the I-Pack to its original position in the intermediate result image data I _Φ ( x ', y ', z ). This process is repeated by looping the y _r and x _r loops until the intermediate result image I _Φ ( x ', y ', z ) is completely reconstructed. Intermediate result image I _Φ ( x ', y ', z ) is reconstructed for each view angle, Φ , and each I _Φ ( x ', y ', z ) has its own view angle Φ The final image I, which is then rotated and transformed together and added together to produce an image in the original Cartesian coordinate system. Reconstruct ( x , y , z ) (see Equation 10). As with the forward projection, in order to maximize the data input / output speed during the calculation, the size of the SIMD data set (I-Pack + P-Pack) and I-Pack is designed to be smaller than the cache size of L1 and L2, respectively.

As shown in Figures 7 and 8, in accordance with one embodiment of the present invention, a rotating operation may be used to use an aligned frame. Prior to the projection step, the projection plane is aligned to the aligned frame, and the incoming image plane to be used in the forward or backward projection is aligned to the aligned frame by rotating according to each viewing angle (which is You can say that it converts from x, y, z coordinate system to x ', y', z coordinate system). After the backward projection step, the original image is recovered from the aligned frame by rotating in the reverse direction from the projection process by the original viewing angle. The image data in the rotated state before being restored to the original image coordinate system is called an intermediate result image.

To reduce process time, effective use of cache is an important factor in the overall calculation. Memory allocation and loop order should be optimized according to memory access patterns. By selecting x _r ( = x ' ) as the outermost variable, the memory stride at the data I / O during the operation can be reduced to within the L2 cache size. Moreover, since most computationally intensive tasks in the SSP method are related to the interpolation of z or y _r , the variables in the innermost loop are selected z or y _r to minimize memory stride during the calculation process. For such memory input / output optimization, data structures used in image reconstruction or backward projection operations are most preferably assigned in the order of [ x ' ]: [ y' ]: [ z ], and in the case of forward projection, [ Φ ]: [ x _r ]: [θ]: [ y _r ] in order. This allocation of memory structure also affects the loop order. In other words, for backward projection operations, if the data structure used in backward projection operations is assigned in the order [ x ' ]: [ y' ]: [ z ], then the backward projection operation is where the z loop is placed at the innermost and x ' The loop will contain the outermost placed execution loop, and for forward projection operations, if the data structure is assigned in the order [ Φ ]: [ x _r ]: [θ]: [ y _r ], then y _r The innermost loop and y _r in the order of θ loop, x _r loop and Φ loop It will contain the execution loop placed outward from the loop. As the most optimized case, the loop order of the forward / backward projection is shown in FIGS. 7 and 8. Also, as a result of applying the previously mentioned symmetries, the total number of loop counts can be reduced by half at each step.

For SIMD, the memory packing step is performed before each forward / reverse projection step. After the projection step, the compressed SIMD data set on the projection area undergoes an unpacking step and then returns to its original position. Similarly, the compressed SIMD data set on the image area is unpacked after each backward projection step.

Experiment result

A. Method

To evaluate the effectiveness of the SSP method, the SSP method was applied to HRRT (High Resolution Research Tomograph) data developed by CPS / Siemens, Knoxville, Tennessee, USA. The HRRT has 119,808 sensors and is designed for ultra high resolution brain scanning. HRRT has the highest computational complexity in image reconstruction because it has the largest number of (smallest) detectors used in human PET scanners in recent years. The number of calculations has increased significantly (in proportion to the number of LORs) and is much larger than other existing PET scanners. The reconstruction software provided in the HRRT package supports an iterative algorithm such as OP-OSEM3D, which is the preferred reconstruction algorithm in HRRT systems. One of the advantages of OP-OSEM3D is that it prevents the generation of negative bias when the random rate is high. Since this algorithm requires the largest memory size and is the most computationally intensive algorithm, we adopt this algorithm to test the performance of the new SSP method. Hereinafter, unless otherwise stated in the present specification, the HRRT OP-OSEM3D package S / W is represented by a conventional method.

As disclosed, the SSP method is supported by a commercial computer system such as a general PC. In this specification, three different platforms are adopted to compare the SSP method with the existing method. The first platform is PC1 with Intel dual-core 3.0GHz, 4GB RAM and 1MB L2 cache per CPU. The second platform PC2 is a high-performance PC with two dual-core AMD 2.4GHz, 8GB RAM and 1MB L2 cache per CPU. The third platform consists of an eight node cluster system (denoted HRRT CPS system), each node being a current HRRT computer platform consisting of two Xeon 3.06GHz processors with 512KB L2 cache and 2GB RAM.

The concept of "span" in 3D means axial compression mode. Applicants performed the tests in Span 3 and Span 9 which are commonly used in HRRT. The parameters for OP-OSEM3D are as follows. 256x256x207 image pixels, 6 iterations and 16 subsets. Thread programming is also used.

B. Results

Before discussing the computational speed, the accuracy of the reconstructed image must be guaranteed (compared to the original HRRT package provided by CPS / Siemens). First, Applicants compared the results of SSP's forward and backprojectors with existing HRRT packages. As shown in FIGS. 9 and 10, there was little difference. 9b shows a profile with projection data at 45 ° and 90 °. The same comparison was made for the backward projection data as well as the final reconstructed image, the results of which are shown in FIGS. 10 and 11. As shown, the projected data and image quality of the SSP are equivalent to the existing method, but have a gain of almost two orders of magnitude in computational speed.

Applicants compared the performance aspects of the SSP method, that is, the execution times of PC1 and PC2. Computation time was measured for two conventional compression modes, Span 3 and Span 9. The results, as predicted, the conventional method compared to the (original HRRT package) it is almost 80 times the SSP method (that is, 2 (Φ + 90 ° symmetry) x 2 (x _r or mirror symmetry) x 2 (y '- Symmetric ) x 2 (θ - symmetry) x 4 (SIMD) + other optimizations). It should be noted that this performance improvement can be obtained somewhat independently of the system architecture. In addition, as disclosed in Tables 3 and 4, the method according to one embodiment of the present invention will have more related improvements as more oblique rays or angles are used. The more oblique rays or angles are used, the better the axial resolution and statistics. Since reconstruction using the EM technique consists of a large number of forward and backward projections, this speed improvement is particularly important. Finally, we conducted a comparative study with a HRRT CPS system with the original HRRT CPS algorithm and PC2 in a 2 dual core CPU (no cluster) configuration performing the SSP method, with a speed gain factor of 9 to 11. Obtained. The reason for the difference between our calculation (64 times) and the actual system based calculation (9-11 times) is probably due to the original HRRT CPS system based on 16 Xeon-CPU. The HRRT CPS system with 16 Xeon-CPUs is much more powerful than the PC2 with only two dual-core CPUs. This difference in H / W may have resulted in a reduction in speed gain. In addition, the SSP method can be easily extended to cluster versions.

The Applicant expects that the SSP method will improve the overall computation time, particularly when dynamic functional studies are needed. Recently, image reconstruction time for a HRRT PET scan using a common commercial PC (PC2 type) needs about 7 to 8 minutes after the completion of cosine curve formation and precorrection in the case of span3. On the other hand, the current HRRT reconstruction package (OP_OSEM3D) using an eight node cluster system requires 80 minutes.

9A shows a comparison of projection data between the existing method and the proposed SSP method. The upper part shows the existing method, the center is the proposed SSP method, and the lower part shows the difference between the SSP method and the existing method. 9B is a cross-sectional view of the cosine curve at a particular field of view. The top shows the comparison at Φ = 45 °, the middle at Φ = 90 °, and the bottom at Φ = 45 ° and Φ = 90 °, respectively.

10A shows a brief comparison of the backward projection image between the existing method and the proposed SSP method. The upper part shows the existing method, the center is the proposed SSP method, and the lower part shows the difference between the SSP method and the existing method. 10B is a cross-sectional view (contour) on the x-axis (y = 154, z = 103). The top shows the comparison between the existing method and the proposed SSP method, and the bottom shows the difference between them.

Figure 11 shows the comparison of the reconstructed images repeated six times and the differences. 11A shows an image set reconstructed by the existing method and the proposed SSP method and their differences. 11B is a cross-sectional view (contour) on the x-axis (y = 154, z = 103).

Table 3 shows a comparison of the execution time of various operations on PC1.

Table 4 shows a comparison of the execution time of various operations on PC2.

According to the present invention, there is provided a fast method for calculating the ray path, forward projection and reconstructed three-dimensional images using the symmetrical properties of the projection and image data. With this simple geometrical symmetry and with the help of SIMD, the present invention achieves computational speed gains for existing reconstruction packages, especially HRRT, without compromising image quality, compared to the recently most widely used PET method, OP-OSED3D. It provides a faster reconstruction method by digit size. The main concepts introduced in the SSP method can be summarized as follows. First, using the concept of rotation-based projection or alignment frames, interpolation operations are reduced from three dimensions to one dimension. Second, the rotation-based projection is combined with the symmetry characteristic of ( θ , Φ , y ' , x _r ) and with the SIMD technique using the optimized L1, L2 cache. By grouping 16 symmetry points that share the same interpolation coefficients (or their complements), they are processed simultaneously in forward / reverse projection operations, which significantly improves the overall reconstruction time. To simultaneously process these symmetry points in forward / reverse projections, SIMD operations have been adopted. Specifically, the SIMD technique makes it possible to access four data simultaneously. In addition, to minimize memory stride, we optimized the data size per loop suitable for the L2 cache size as well as the data structure.

In summary, the present invention enables forward and backward projection operations using a symmetric concept and cache optimized SIMD on an alignment frame, thereby dramatically speeding up the reconstruction of images from cosine curve data measured from scanners of PET or CT. Improved. In other words, the SSP method proposed in the present invention incorporates the symmetry characteristic of the projection to reduce the computation time as much as 16 times, and achieves an additional computation time reduction of about 4 to 5 times by using the SIMD operator and cache optimization. As a result, the overall calculation speed gain in image reconstruction is almost 80 times higher than that of the conventional method. This improvement in computation speed enables the use of high resolution PET, such as HRRT, to perform dynamic functional studies. Until now, in the case of dynamic PET research, there have been many problems in using ultra-high resolution PET due to the excessively long reconstitution time, which has been a major obstacle in research using molecular imaging.

The invention will provide new opportunities for PET users, especially researchers using molecular imaging, or anyone who wants to use ultra-high resolution PET for diagnosis and treatment. The dynamic PET research presented above is one of them. The proposed SSP method is also applicable to all processes that require forward and backward projection. For example, it can be applied immediately to precorrection processes such as attenuation and scattering correction in PET. These advantages mean that the SSP method proposed by the present invention can be used to realize a truly interactive scanner with real-time image reconstruction. In addition, the SSP method of the present invention has an advantage that the existing PET can be easily applied to a cluster system and the like adopted for image reconstruction.

1A shows a projection of a 3D object and the object onto a 2D projection plane.

1B is a y'-z top view of FIG. 1A;

1C shows an example of linear integration along the horizontal plane at θ = 0 ° and the y ′ line projected on x _r ;

2A shows the rotation of a projection ray (ray path) or frame over an image plane on a fixed reference coordinate (x, y).

FIG. 2B is a diagram showing a case where the projection ray (ray path) frame coincides with fixed coordinates (x ', y').

및

사이의 관계를 도시한 도면.3 illustrates z, y ', x', θ, x _r , y _r , using a y'-z plan view.

And

Figure showing the relationship between.

4A shows an example of "mirror symmetry" used in the proposed SSP method.

4b shows an example of "φ-symmetry" used in the proposed SSP method.

5A shows the y'-symmetry of the proposed SSP method.

5b shows θ-symmetry of the proposed SSP method.

5C shows the relationship between y'-symmetry and θ-symmetry or the overall symmetry obtainable by y'-symmetry and θ-symmetry.

6 illustrates the concept of uniform work distribution based on the sum of the beam path lengths.

7 is a flow chart of a forward projection in accordance with an embodiment of the present invention.

8 is a flow chart of backward projection in accordance with an embodiment of the present invention.

9a shows a comparison of projection data between an existing method and a proposed SSP method.

9B is a cross-sectional view of the cosine curve at a particular field of view.

10a shows a comparison of a brief backward projection image between an existing method and a proposed SSP method.

10B is a cross-sectional view (contour) of the x-axis (y = 154, z = 103).

FIG. 11A illustrates an image set reconstructed by the existing method and the proposed SSP method and differences thereof. FIG.

11B is a sectional view (contour) of the x-axis (y = 154, z = 103).

Claims (16)

As a method of reconstructing an image of 3D PET, A backward projection step of converting first cosine curve data including data for a plurality of projection angles into image data, and A forward projection step of converting the image data into second cosine curve data, The reverse projection step, For each of the plurality of projection angles, an alignment frame with respect to the projection angle using a projection ray, wherein the alignment frame with respect to the projection angle has the axes x, y of the image coordinate system (x, y, z) by the corresponding projection angle. A projection plane having an image space having a rotated coordinate system (x ', y', z) and a cosine curve coordinate system (x _r , y _r, Φ, θ) having an axis x _r parallel to x 'of the rotated coordinate system. Wherein the axis y 'of the coordinate system of the alignment frame is parallel to the projection rays-an alignment frame for the corresponding projection angle by performing a backward projection operation on the first cosine curve data having the cosine curve coordinate system. Generating first intermediate image data having a coordinate system (x ', y', z) of-, Generating second intermediate image data having the image coordinate system (x, y, z) by rotating the first intermediate image data for each of the plurality of projection angles by the corresponding projection angle, and Adding the second intermediate image data for the plurality of projection angles to generate the image data, The forward projection step, Forward projection intermediate image data for the plurality of projection angles by placing the image data having the image coordinate system (x, y, z) in each of the alignment frames for the plurality of projection angles-the forward projection intermediate image Data each having a coordinate system (x ', y', z) of the alignment frame for a corresponding projection angle; and The second cosine curve having the cosine curve coordinate system by performing a forward projection operation on the forward projection intermediate image data using the projection light beam on the alignment frame for the corresponding projection angle for each of the plurality of projection angles. Generating data; The backward projection operation and the forward projection operation each use backward projection interpolation coefficients and forward projection interpolation coefficients, and a plurality of symmetric coordinate pairs having a symmetric relationship on the alignment frame, wherein the symmetric coordinate pairs are coordinates on the image coordinate system and the The backward projection interpolation coefficient and the forward projection interpolation coefficient computed for one symmetric coordinate pair of-as the backward projection interpolation coefficient and the forward projection interpolation coefficient for the plurality of symmetric coordinate pairs. Use, The backward projection operation uses a single instruction multiple data (SIMD) operation and uses the plurality of symmetric coordinate pairs to reduce the number of times a SIMD backward projection instruction is called in the backward projection operation. One SIMD backward projection operation instruction, which is divided and packed into at least one SIMD projection data set and called for the backward projection operation, processes at least one of the SIMD projection data sets, The forward projection operation uses a SIMD operation and packs the plurality of symmetric coordinate pairs into one or more SIMD image data sets to reduce the number of times a SIMD forward projection instruction is called in the forward projection operation, and packs the forward projection. And one SIMD forward projection calculation instruction called for operation processes at least one of said SIMD image data sets. The method of claim 1, Wherein the symmetry relationship comprises one or more of a mirror symmetry relationship, φ-symmetry relationship, θ-symmetry relationship, and y'-symmetry relationship. The method of claim 1, The plurality of symmetrical coordinate pairs symmetrical on the alignment frame include sixteen symmetrical coordinate pairs determined through mirror symmetry, φ-symmetry, θ-symmetry, and y'-symmetry, In the backward projection step, at least one of a representative backward projection interpolation coefficient and a complement of the representative backward projection interpolation coefficient computed for one symmetric coordinate pair of the sixteen symmetric coordinate pairs is determined for the sixteen symmetric coordinate pairs. Used as the backward projection interpolation coefficient, In the forward projecting step, at least one of a representative forward projection interpolation coefficient computed for one symmetric coordinate pair of the sixteen symmetric coordinate pairs and the complement of the representative forward projection interpolation coefficient for the sixteen symmetric coordinate pairs And reconstructing an image of 3D PET using the forward projection interpolation coefficient and the forward projection interpolation coefficient. The method of claim 3, The generating of the first intermediate image data may include packing the sixteen symmetric coordinate images into four SIMD projection data sets such that four pairs of symmetric coordinates are included in one SIMD projection data set, and performing the backward projection operation. And one SIMD backward projection operation instruction per one SIMD projection data set in order to be called. The method of claim 3, Generating the second cosine curve data may include packing the sixteen symmetrical coordinate pairs into four SIMD image data sets such that four pairs of symmetric coordinates are included in one SIMD image data set and the forward projection operation. And one SIMD backward projection operation instruction invoked for processing at least one of said SIMD image data sets. The method of claim 1, In the backward projection step, the first cosine curve data is allocated to the memory in the order of x ', y', z in the alignment frame, and the backward projection operation includes a z loop disposed at the innermost and an x 'loop positioned at the most. A method of reconstructing an image of 3D PET, including an outer loop disposed outward. The method of claim 1, In the forward projection step, the forward projection intermediate image data is allocated to the memory in the order of Φ , x _r , θ y _r on the alignment frame, and the forward projection operation is y _r. The loop is placed at the innermost and y _r in the order of θ loop, x _r loop, and Φ loop A method of reconstructing an image of 3D PET, comprising a performance loop disposed outwardly from the loop. The method of claim 1, The backward projection step and the forward projection step are performed in a multiprocessor system, and based on the sum of the beam path lengths of the projection beams, the processor is arranged along the x 'axis of the alignment frame such that an equal amount of work is allocated to each processor. A method of reconstructing an image of 3D PET, which allocates a workload of a backward projection operation and the forward projection operation. Device for reconstructing the image of 3D PET, Means for performing backward projection to convert first cosine curve data including data for a plurality of projection angles into image data, and Means for performing forward projection to convert the image data into second cosine curve data, The means for performing the reverse projection, For a plurality of projection angles, an alignment frame using projection rays, the alignment frame rotates the axes x, y of the image coordinate system (x, y, z) by a corresponding projection angle (x ', y', z) And a projection plane having a cosine curve coordinate system (x _r , y _{r ,} Φ, θ) having an image space having an axis x and an axis x _r parallel to x 'of the rotational coordinate system, wherein the axis y of the coordinate system of the alignment frame 'Is parallel to the projection beams-coordinate system (x', y ', z) of the alignment frame for each projection angle by performing a backward projection operation on the first cosine curve data having the cosine curve coordinates. A first intermediate image data generation unit for generating first intermediate image data having- A second intermediate image data generator for generating second intermediate image data having the image coordinate system (x, y, z) by rotating the first intermediate image data for the plurality of projection angles by the corresponding projection angle, respectively; An image combining unit generating the image data by adding the second intermediate image data for the plurality of projection angles; Means for performing the forward projection, The image data having the image coordinate system (x, y, z) is disposed on the alignment frame for each of the plurality of projection angles so that the coordinate system (x ', y', on the alignment frame for each of the plurality of projection angles) a forward projection intermediate image data generator for generating forward projection intermediate image data having z), and A second cosine curve data having the cosine curve coordinate system by performing a forward projection operation on the forward projection intermediate image data using the projection light beam on the alignment frame, for the plurality of projection angles; Including cosine curve data generation unit, The means for executing the backward projection and the means for executing the forward projection respectively use backward projection interpolation coefficients and forward projection interpolation coefficients in executing the backward projection operation and the forward projection operation, respectively, and are symmetrical on the alignment frame. The backward projection interpolation coefficient and the forward projection computed on one symmetric coordinate pair of a plurality of relational symmetric coordinate pairs, wherein the symmetric coordinate pair consists of the image data and the first or second cosine curve data. An interpolation coefficient is used as the backward projection interpolation coefficient and the forward projection interpolation coefficient for the plurality of symmetric coordinate pairs, The means for performing the backward projection uses a single instruction multiple data (SIMD) operation and the plurality of symmetric coordinate pairs to reduce the number of times a SIMD backward projection instruction is called in the backward projection operation. Dividing into one or more SIMD projection data sets and packing, causing one SIMD backward projection operation instruction to be called per one of the SIMD projection data sets for the backward projection operation, The means for performing the forward projection uses a SIMD operation and packs the plurality of symmetric coordinate pairs into one or more SIMD image data sets to reduce the number of times a SIMD forward projection instruction is called in the forward projection operation, And reconstruct an image of 3D PET such that one SIMD forward projection operation instruction is called per one SIMD image data set for the forward projection operation. 10. The method of claim 9, Wherein the symmetry relationship comprises one or more of a mirror symmetry relationship, a φ-symmetry relationship, a θ-symmetry relationship, and a y'-symmetry relationship. 10. The method of claim 9, The plurality of symmetric coordinate pairs symmetrical on the alignment frame include sixteen symmetric coordinate pairs determined through mirror symmetry, φ-symmetry, θ-symmetry, and y'-symmetry, The means for performing the backward projection is characterized in that at least one of the complement of the representative backward projection interpolation coefficient and the representative backward projection interpolation coefficient computed for one symmetric coordinate pair of the sixteen symmetric coordinate pairs Used as the backward projection interpolation coefficient for the pair, The means for performing the forward projection comprises: performing at least one of the complement of the representative forward projection interpolation coefficient and the representative forward projection interpolation coefficient computed for one symmetric coordinate pair of the sixteen symmetric coordinate pairs and the sixteen symmetric coordinates. And reconstruct an image of 3D PET using the forward projection interpolation coefficient and the forward projection interpolation coefficient for a pair. The method of claim 11, The means for performing the reverse projection packs the sixteen symmetric coordinate images into four SIMD projection data sets such that four pairs of symmetric coordinates are included in one SIMD projection data set, one for the backward projection operation. And reconstruct an image of a 3D PET that invokes one SIMD backward projection operation instruction per SIMD projection data set. The method of claim 11, The means for performing the forward projection packs the sixteen symmetric coordinate pairs into four SIMD image data sets such that the four symmetric coordinate pairs are included in one SIMD image data set, one for the forward projection operation. And reconstruct an image of 3D PET invoking one SIMD forward projection computation instruction per SIMD image data set. 10. The method of claim 9, The means for performing the backward projection allocates the first cosine curve data to the memory in the order of x ', y', z in the alignment frame, with a z loop placed at the innermost and an x 'loop positioned at the outermost. And reconstruct an image of a 3D PET that executes the backward projection operation including a placed execution loop. 10. The method of claim 9, The means for performing the forward projection allocates the forward projection intermediate image data to memory in the order of Φ , x _r , θ y _r on the alignment frame, and y _r The loop is placed at the innermost and y _r in the order of θ loop, x _r loop and Φ loop And reconstruct an image of 3D PET that executes the forward projection operation including a performance loop disposed outwardly from the loop. 10. The method of claim 9, The means for performing the backward projection and the means for performing the forward projection operate in a multiprocessor system and x of the alignment frame such that an equal amount of work is assigned to each processor based on the sum of the beam path lengths of the projection beams. An apparatus for reconstructing an image of 3D PET, which allocates the workload of the backward projection operation and the forward projection operation along an axis.

Images