KR101283266B1 - Method and apparatus for monte-carlo simulation gamma-ray scattering estimation in positron emission tomography using graphics processing unit - Google Patents

Abstract

A method and apparatus are provided for generating a gamma ray scattering sinogram to reconstruct an image. In the reconstructed image without initial scattering correction, the photon generation position is determined. The GPU uses Monte Carlo simulation to compute the photon's movement. The three-dimensional measurement line is calculated from the position where the photon reaches the detector, and the calculated measurement line is converted into the three-dimensional sinogram form and stored in the prompt sinogram and gamma ray scattering sinogram. The final gamma ray scattering sinogram is generated by scaling the gamma ray scattering sinogram.

Description

METHOD AND APPARATUS FOR MONTE-CARLO SIMULATION GAMMA-RAY SCATTERING ESTIMATION IN POSITRON EMISSION TOMOGRAPHY USING GRAPHICS PROCESSING UNIT}

The examples below relate to Positron Emission Tomography (PET).

A method and apparatus are disclosed for processing Monte Carlo simulations in parallel using a GPU for gamma ray scattering estimation.

Tomography is an imaging technique in which a detector is rotated 180 ° or 360 ° to observe a particular section of an object non-invasively without overlapping another adjacent section.

A sinogram refers to a plurality of images obtained by tomography. The sinogram is reconstructed into a 3D image. An observer can observe a specific cross section of the photographed object through the reconstructed three-dimensional image.

In nuclear medicine tomography, radiopharmaceuticals administered to the subject's body are represented by images of the distribution of the radiopharmaceutical in the subject's body according to their biochemical properties.

Radiopharmaceuticals administered to the body are limited in radiation dose by reference values. Therefore, it is difficult to reconstruct a high resolution image from the projection value recorded by the very small amount of photons.

When imaging a tomography using nuclear medical imaging equipment, the noise is relatively strong compared to the weak intensity of the observed signal.

Therefore, scattering of gamma rays has a great influence on the reconstructed image, and an approximated gamma scattering estimation algorithm, such as a single scattering simulation, is widely used to correct image noise due to scattering of gamma rays.

Monte Carlo simulation is a method of generating gamma ray photons one by one and simulating the migration of gamma ray photons similarly to the real world.

Monte Carlo simulation is a way for researchers to predict and evaluate what can actually happen. However, Monte Carlo simulations typically take hours to days of computational time. Thus, Monte Carlo simulations are rarely used in actual clinical practice. Similarly, gamma ray scattering estimation and correction using Monte Carlo simulations are rarely used in research.

In the process of reconstructing a Positron Emission Tomography (PET) image into a 3D image, gamma ray scattering estimation process and the final 3D image reconstruction process using the same are the most time-consuming.

Parallel processing techniques can be used to reduce this computation time. Parallel processing is a method of dividing the same kind of work in a plurality of computing cores at the same time.

A multi-core central processing unit (CPU) or GPU can be used for parallel processing. In particular, GPUs typically have dozens or hundreds of cores to perform massive processing in parallel, and can handle mathematical operations as well as graphic operations.

One embodiment of the present invention, in the gamma ray scattering estimation using the Monte Carlo simulation known to be the most accurate, can provide a method and apparatus having a speed several hundred times more improved by using a GPU.

An embodiment of the present invention may provide a method and apparatus for estimating gamma ray scattering using Monte Carlo simulation having an improved speed to a level that can be used even in a clinic.

One embodiment of the present invention may provide a method and apparatus for generating a gamma ray scattering sinogram using a hybrid CPU-GPU structure.

According to one aspect of the invention, the step of receiving the actual measured sinogram generated by measuring the object, performing a Monte Carlo simulation for the initial occurrence of the gamma rays generated in the object-the Monte Carlo simulation is GPU Generating a gamma ray scattering sinogram with gamma ray scattering corrected based on the information about the scattered photons generated by the Monte Carlo simulation, and between the gamma ray scattering An image reconstruction method is provided, which includes generating a reconstructed image based on a nogram.

According to another aspect of the present invention, there is provided a method of determining a generation position of a photon in a reconstructed image, calculating a movement of two photons generated at the generation position by using Monte Carlo simulation, two detector positions, and a movement. Calculating the changed values of the detectors, wherein the two detector positions are detector coordinates when each of the two photons reaches the detector, and the changed values during the movement are determined by the movement of each of the two photons through and through the object. Are values changed during the step of: calculating a three dimensional measurement line connecting the two detector positions, generating a prompt sinogram and a gamma ray scattering sinogram based on the three dimensional measurement line and the gamma ray scattering sino A method for generating gamma ray scattering sinograms is provided, the method comprising scaling a gram. .

The generation position may be a plurality.

Each of the plurality of generation positions may be assigned to one core of the plurality of cores of the GPU.

Calculating the movement of the two photons generated at the generation position by using the Monte Carlo simulation and calculating the two detector positions and the values changed during the movement may be performed in parallel by the plurality of cores. .

Determining a generation position of a photon in the reconstructed image includes: recognizing positions having values in the reconstructed image as the plurality of occurrence positions and generating a table including the plurality of occurrence positions. It may include.

Each of the plurality of cores of the GPU may be assigned one or more occurrence positions determined based on the number of cores in the table.

Computing the movement of the two photons generated at the generation position by using the Monte Carlo simulation, wherein the core of the GPU is assigned the generation position, the core is to determine the initial movement direction at the assigned generation position of photons Randomly generating, the core randomly generating a moving distance of the photons, the core calculating a position at which the photons move by the moving distance in the initial moving direction, and the core is moved Calculating the attenuation of the photon as it moves into position and confirming that the photon has reached the detector.

Computing the attenuation of the photon according to the movement of the core to the moved position, the attenuation of the photon is a value corresponding to the three-dimensional coordinates of the moved position in the three-dimensional image of the attenuation coefficient of the core The core may be acquired by reading a value corresponding to the three-dimensional coordinates by using a three-dimensional linear interpolation method provided by the GPU.

Calculating the movement of two photons generated at the generation position by using the Monte Carlo simulation, the core determines whether scattering of the photons at the moved position occurs, if the scattering has occurred, The core may further include randomly generating a scattering angle of the scattering, and when the scattering occurs, the core calculates the Compton scattering value of the scattering.

Determining whether the scattering of the photons at the moved position of the core is determined, if the attenuation coefficient of the moved position of the core is less than a boundary value, the scattering does not occur, And determining that scattering has occurred if the attenuation coefficient of a location is greater than the boundary value.

Generating a prompt sinogram and a gamma ray scattering sinogram based on the 3D measurement line may include converting the 3D measurement line into a sinogram form, and converting the 3D measurement line into the sinogram form. Storing a line in a prompt sinogram and storing a three-dimensional measurement line converted to the sinogram form in the gamma ray scattering sinogram when the two photons are scattered photons. .

The converting of the 3D measurement line into a sinogram form may include: dividing the 3D measurement line into a vector in an XY direction and a Z direction, and converting the 3D measurement line into the sinogram form based on the vector. It may include the step of converting.

Adjusting the scale of the gamma ray scattering sinogram comprises: calculating a scale constant that minimizes the error between the actual measured sinogram and the prompt sinogram and adding the scale constant to the gamma ray scattering sinogram And scaling the gamma ray scattering sinogram by multiplying.

According to another aspect of the present invention, a receiver for receiving an actual measured sinogram generated by measuring an object, a CPU for generating a first reconstructed image using the actually measured sinogram and a gamma ray in the object A GPU that performs Monte Carlo simulations for these initial occurrence locations; Wherein the GPU performs the Monte Carlo simulation in parallel by initial generation positions of photons by using one or more cores, and wherein the CPU scatters gamma rays based on information about scattered photons generated by the Monte Carlo simulation An image reconstruction apparatus is provided, which generates this corrected sinogram and generates a reconstructed image based on the sinogram whose gamma ray scattering is corrected.

According to yet another aspect of the present invention, a generation position of a photon in a reconstructed image is determined, a three-dimensional measuring line connecting two detector positions, a prompt sinogram based on the three-dimensional measuring line and Compute the movement of the two photons generated at the generation position by using a Monte Carlo simulation and a CPU to generate a gamma ray scattering sinogram and scale the gamma ray scattering sinogram, the two detector positions and GPU for calculating changed values during movement-the two detector positions are detector coordinates when each of the two photons reaches the detector, and the changed values during the movement are each of the two photons passing through the object and moving Gamma-ray scattering sinogram generation apparatus is provided, wherein the values vary over time.

According to another aspect of the invention, the CPU determines the location of the photon's occurrence in the reconstructed image, the GPU three-dimensionally connecting two detector positions for the pair of photons generated at the generation position by using Monte Carlo simulation Generating information of a measurement line, wherein the detector positions are detector coordinates when each of the two photons of the photon pair reaches the detector, and the CPU rebinds the sinogram based on the information of the three-dimensional measurement line And generating a plurality of generation positions, and assigning each of the plurality of generation positions to each of a plurality of threads of the GPU, generating information of the 3D measurement line; And performing the sinogram rebinding is performed in parallel by the CPU and the GPU, respectively. This method is provided.

The plurality of generation positions may be divided into a plurality of groups.

Each of the plurality of groups may be processed by iterations of generating information of the 3D measurement line and performing the sinogram rebinding.

Information of the 3D measurement line may be transmitted from the GPU to the CPU through an asynchronous memory transfer.

According to another aspect of the invention, by using the CPU and Monte Carlo simulation to perform the sinogram revision based on the information of the three-dimensional measurement line, to determine the location of the photon in the reconstructed image at the location of the generation A GPU for generating information of the three-dimensional measurement line joining two detector positions for the generated photon pair, wherein the detector positions are detector coordinates when each of the two photons of the photon pair reaches the detector; There are a plurality of generation positions, and the GPU assigns each of the plurality of generation positions to each of a plurality of threads of the GPU, generating information of the 3D measurement line, and performing the sinogram rebinding. Steps are provided by a sinogram generating apparatus, each executed in parallel by the CPU and the GPU.

A gamma ray scattering estimation method and apparatus are provided in which computation time is reduced by using parallel processing.

A method and apparatus for generating a Monte Carlo simulation scattering sinogram using parallel processing of a GPU is provided.

A method and apparatus are provided for generating gamma ray scattering sinograms using a hybrid CPU-GPU architecture.

1 is an image reconstruction apparatus using a gamma ray scattering sinogram according to an embodiment of the present invention.
2 is a flowchart illustrating an image reconstruction method using a gamma ray scattering sinogram according to an embodiment of the present invention.
3 is a flowchart illustrating a Monte Carlo simulation line according to an embodiment of the present invention.
4 is a flowchart of a sinogram generation method according to an embodiment of the present invention.
5 is a flowchart illustrating a method of adjusting a scale of a gamma ray scattering sinogram according to an embodiment of the present invention.
6 is a flowchart illustrating a method of performing an operation according to movement of photons until reaching a detector according to an example of the present invention.

Hereinafter, an embodiment of the present invention will be described in detail with reference to the accompanying drawings. However, the present invention is not limited to or limited by the embodiments. Like reference symbols in the drawings denote like elements.

1 is an image reconstruction apparatus using a gamma ray scattering sinogram according to an embodiment of the present invention.

The image reconstruction apparatus 100 includes a receiver 110, a central processing unit (CPU) 120, and a graphics processing unit (GPU) 130. The image reconstruction apparatus 100 may be a sinogram generating apparatus that generates a gamma ray scattering sinogram.

The receiver 110 may receive data representing a sinogram actually measured through a wired or wireless network.

GPU 130 may include one or more cores 140. The cores 140 may perform graphics operations in parallel.

GPU 130 calculates the Monte Carlo simulation. Here, the Monte Carlo operation means calculating the photon movement.

Each of the cores 140 calculates the movement of the gamma ray photon pair (ie, the path of travel of the photon pair) that occurred at one initial generation position.

Movement of the photon pairs can include scattering. The photon's state value includes information about whether or not the photon is scattered.

When scattering occurs, the scattering angle and the moving distance of the photons are determined randomly.

When the photons reach the detector, the operation ends.

GPU 120 may be included in CPU 120.

CPU 120 may generate a prompt sinogram and gamma ray scattering sinogram by using the coordinates of the detector obtained by Monte Carlo simulation.

The prompt sinogram is a sinogram that contains all the data. Gamma ray scattering sinograms are sinograms that contain only gamma ray scattering.

In addition, the CPU 120 may calculate a scale constant having the least error by comparing the measured sinogram and the prompt sinogram, and multiply the calculated scale constant by the gamma ray scattering sinogram to obtain the final gamma ray scattering sinogram. Gram can be obtained.

The image reconstruction apparatus 100 may reconstruct an image using a gamma ray scattering sinogram.

In addition, a gamma ray scattering sinogram may be generated by the CPU 120 and the GPU 130. Therefore, the CPU 120 and the GPU 130 may configure a gamma ray scattering sinogram generating apparatus.

2 is a flowchart illustrating an image reconstruction method using a gamma ray scattering sinogram according to an embodiment of the present invention.

This embodiment may be performed by the image reconstruction apparatus 100. That is, the image reconstruction apparatus 100 including the GPU may reconstruct an image using a gamma ray scattering sinogram.

The reconstructed image may be a 3D image.

In step 210, for example, by the receiver 110, a sinogram for reconstructing an image is received.

The received sinogram is the actual measured sinogram generated by the real machine (ie, positron emission tomography) measuring the object.

The received data may be from a positron emission tomography camera. Alternatively, the received data may be transmitted from a server for providing data.

By using the actually measured sinogram, the CPU 120 may generate a reconstructed image without scattering correction, and an object may be recognized based on the reconstructed image.

In step 220, for example, by the GPU 130, a Monte Carlo simulation is performed on the initial occurrence locations where gamma rays occur in the object.

The initial occurrence locations may be locations within the reconstructed image where scattering is not corrected. That is, for example, the CPU 120 may determine initial generation positions at which gamma rays occur in the reconstructed image in which scattering is not corrected.

Monte Carlo simulation is a simulation in which virtual photons are generated to move the virtual photons generated randomly, similar to actual ones. In other words, Monte Carlo simulation is to model the actual photon movement.

Monte Carlo simulation can be implemented through parallel processing on GPU 130. That is, the Monte Carlo simulation may be performed in parallel for each of the initial generation positions of the photons by using the GPU 130.

Parallel processing is implemented by parallel processing code.

Even for parallel processing codes that implement the same algorithm, the operation speed of each parallel processing code is different from each other depending on the complexity inside the code.

Therefore, parallel processing code needs to be simply constructed using a small number of variables.

For example, a CPU (eg, the CPU 120) in the image reconstruction apparatus 100 performing Monte Carlo simulation is involved only in file input / output, etc., and the actual operation for the Monte Carlo simulation is performed in the parallel processing apparatus (eg, the GPU 130). If divided and processed in parallel, the parallel processing code executed in such a parallel processing apparatus needs to be composed of the most concise code.

In order to improve the operation speed of parallel processing code which is the core of parallel processing, it is very important to reduce the frequency of data transmission and reception to a minimum. In general, one operation can be completed within a few nanoseconds, but operations that read or write memory can take hundreds of nanoseconds or more.

If the received data is from a positron emission tomography, the sinogram required to reconstruct the image may only be received once before the operation begins, and the reconstructed image may only be received once after the operation is completed. Can be.

For step 220, it is described in detail below with reference to FIG.

In step 230, for example, a gamma ray scattering sinogram is generated by the CPU 120.

Gamma-ray scattering sinograms are sinograms with gamma-ray scattering corrected.

A gamma ray scattering sinogram can be generated based on the information about the scattered photons generated by the Monte Carlo simulation described above. That is, the CPU 120 may obtain information about scattered photons from the Monte Carlo simulation described above, and generate a gamma ray scattering sinogram based on the obtained information.

Scaling may be required in performing step 230. Scale adjustment is described in detail below with reference to FIG. 4.

In step 240, for example, the CPU 120 generates a reconstructed image with scattering corrected based on the generated gamma ray scattering sinogram.

3 is a flowchart illustrating a Monte Carlo simulation line according to an embodiment of the present invention.

Steps 310-340 may be included in step 220 described above.

In step 310, for example, by the CPU 120, the generation position of the photons in the reconstructed image is determined.

There may be a plurality of generation positions.

The reconstructed image may be a reconstructed image in which the scattering described above with reference to FIG. 2 is not corrected.

Step 310 may include steps 312 and 314.

In step 312, positions with values in the reconstructed image, for example, by the CPU 120, may be generated as the initial occurrence position of the photons.

The generation position may be a gamma ray generation position. Also, the photon may be a gamma ray photon.

In step 314, a table including a plurality of occurrence positions may be generated, for example by CPU 120, by collecting only these occurrence positions. The table may reflect the amount of occurrence of gamma rays (ie the number of occurrence positions).

Each of the plurality of cores 140 of the GPU 130 may read out different positions from each other in a table generated according to the number of the cores 140. This means that each of the plurality of cores 140 of the GPU 130 is assigned one or more occurrence positions among the plurality of occurrence positions in the table.

In other words, each of the plurality of generation positions may be assigned to one core 140 of the plurality of cores 140 of the GPU 130. Here, one or more occurrence locations assigned to the core 140 may be determined based on the number of the core 140.

Each of the plurality of generation positions may be sequentially read out alternately, regardless of the number of cores 140 included in the GPU 130. For example, if there are 60 generation positions and 30 cores 140 of the GPU 130, each of the cores 140 may be assigned two generation positions.

The plurality of cores 140 of the GPU 130 may perform operations (eg, steps 320 and 330) in parallel on one or more occurrence positions assigned to each of the plurality of cores 140. have.

In step 320, for example, by the core 140 of the GPU 130, the movement of the two photons generated at the location of generation by using Monte Carlo simulation is calculated.

For step 320, it is described in detail below with reference to FIG.

In step 330, for example, by the core 140 of the GPU 130, the two detector positions and the values changed during movement are calculated.

The two detector positions are the detector coordinates when each of the two photons reaches the detector.

The values changed during movement are the values changed while each of the two photons passes through and moves through the object.

When positron emission tomography is used, gamma rays occur in pairs. Detector position coordinates are required when each of the two photons reaches the detector as gamma rays occur in pairs.

The sinogram can be obtained by calculating the Line of Response (LoR) through the detector position coordinates.

Therefore, in one operation unit, each core 140 may be responsible for one gamma ray generating position, and output the final two detector positions through physical operation on the gamma ray photon pair generated at the gamma ray generating position. In this case, each core 140 may output values that are changed while the gamma ray photon pair passes and moves through the object (eg, the subject).

In step 340, for example, by the core 140 of the GPU 130, the two detector positions and the values changed during movement are stored.

The core 140 may send two detector positions and values changed during movement to the CPU 120.

4 is a flowchart of a sinogram generation method according to an embodiment of the present invention.

Steps 410-450 may be included in step 230 described above.

In step 410, for example, by the CPU 120, two detector positions each of the two photons has reached and values changed during the movement are received.

In step 420, for example, by the CPU 120, a three-dimensional measuring line (or equation of three-dimensional measuring line) connecting two detector coordinates is calculated.

In step 430, for example, by the CPU 120, a prompt sinogram and a gamma ray scattering sinogram are generated based on the calculated three-dimensional measurement line.

Step 430 may include the following steps 432, 434, 436, 438, and 440.

In steps 432 and 434, for example, by the CPU 120, the generated three-dimensional measurement line may be converted into sinogram format.

In step 432, for example, by the CPU 120, the generated three-dimensional measurement line is divided into vectors in the X-Y direction and the Z direction.

In step 434, for example, by the CPU 120, the three-dimensional measurement line is converted into sinogram form based on the vector.

In step 436, for example, by the CPU 120, the three-dimensional measurement line converted into the sinogram format is stored in the prompt sinogram. The three-dimensional measurement line may be stored in the prompt sinogram through the sinogram coordinates corresponding to the three-dimensional measurement line.

Prompt sinograms include all measurement lines, including scattering. Thus, the prompt sinogram may be considered to be data in the same format as the sinogram actually measured by the machine (ie, positron emission tomography). With this property, the scale adjustment constant can be calculated by comparing the prompt sinogram obtained by Monte Carlo simulation and the sinogram actually measured by the machine. A specific method of calculating the scale composition constant is described in detail below with reference to FIG. 5.

In step 438, for example, by the CPU 120, it is checked whether scattering has occurred while the two photons are moving.

Each of the two photons has a state value. The photon's state value includes information about whether or not the photon is scattered. Therefore, the CPU 120 may check whether scattering has occurred through the state value of the photons.

If scattering has occurred, in step 440, for example, by the CPU 120, the three-dimensional measurement line converted into the sinogram format is stored in the gamma ray scattering sinogram.

In order to estimate and correct gamma ray scattering, a gamma ray scattering sinogram is required. Therefore, scattering sinograms are stored separately in this embodiment.

In step 450, for example, by the CPU 120, the scale of the gamma ray scattering sinogram is adjusted.

For step 450, it is described in detail below with reference to FIG.

5 is a flowchart illustrating a method of adjusting a scale of a gamma ray scattering sinogram according to an embodiment of the present invention.

While photographing the actual positron emission tomography, hundreds of millions of photons generated by the machine (ie, positron emission tomography) cannot be accurately estimated. Thus, the scale of the gamma ray scattering sinogram can be adjusted by first calculating the simulation using a sufficient amount of photons and comparing the results of the simulation with the actual measured values.

In step 510, for example, by the CPU 110, the sinogram actually measured at the machine (ie, positron emission tomography) is received.

In step 520, for example, by the CPU 110, a prompt sinogram obtained through Monte Carlo simulation is received.

In step 530, for example, by the CPU 110, a scale constant is calculated that minimizes the error between both by comparing the measured sinogram and the prompt sinogram.

Since the distribution of the prompt sinogram closely matches the distribution of the actually measured sinogram, a constant is calculated.

Since the actual measured sinogram is the actual data obtained by the machine, it includes all the scattering. Thus, the actual measured sinogram can be compared with the prompt sinogram.

In step 540, for example, by the CPU 110, the scale of the gamma ray scattering sinogram is adjusted by multiplying the scale constant by the gamma ray scattering sinogram.

Through the steps 510 to 540 described above, the gamma ray scattering sinogram may have the same scale as the actually measured sinogram.

6 is a flowchart illustrating a method of performing an operation according to movement of photons until reaching a detector according to an example of the present invention.

Operation based on the movement of photons may be performed in the GPU 130. The core 140 of the GPU 130 may calculate a path of two photons (ie, photon pairs) generated at one generation position.

In step 610, each core 140 is assigned an initial generation location of photons.

In steps 620-690, each core 140 creates an initial path of photons. The initial path means a path from the initial occurrence position in any direction. In addition, each core 140 may generate a path including an arbitrary distance and scattering angle, and compute the moving photons according to the generated path.

Random values, such as a path in any direction, any distance and scattering angle, can be generated by a random number generator. Random values may be generated by using a random generation function provided by the GPU 130.

In step 620, each core 140 randomly generates an initial direction of movement at the assigned generation location of photons.

In step 630, each core 140 randomly generates a moving distance of the photons. That is, each core 140 randomly determines how much to move in a randomly determined movement direction.

In step 635, each core 140 calculates the position at which the photons moved by the moving distance in the initial moving direction.

In step 640, each core 140 determines whether scattering of photons occurs at the location where the photons have moved.

The core 140 may determine that scattering does not occur when the attenuation coefficient of the position where the photon is moved is less than or equal to the boundary value, and determine that scattering occurs when the attenuation coefficient of the moved position is larger than the boundary value.

The attenuation coefficient may be obtained through a three-dimensional image of the attenuation coefficient photographed by a tomography camera.

If scattering has occurred, steps 650 through 660 are executed. If scattering has not occurred, step 570 is executed.

In step 650, each core 140 randomly generates scattering angles.

The angle of scattering can affect the direction of movement of the photons in the next step.

In step 660, each core 140 calculates a Compton Scattering value of scattering. That is, each core 140 calculates the effect of Compton scattering of scattering.

Each core 140 can calculate Compton scattering values by using the Klein-Nishina formula.

The equation used to calculate the Compton scattering value reflects the energy reduction due to the scattering angle.

In step 670, each core 140 calculates the attenuation of the photons as they move to the location where the photons have moved. That is, each core 140 may sum the attenuation coefficients of the object corresponding to the path of the photons.

The attenuation of the photons can be calculated through a three-dimensional image of the attenuation coefficient photographed by the tomography camera described above. That is, the attenuation of the photons may be calculated by the core 140 acquiring a value corresponding to the three-dimensional coordinates of the position where the photons are moved in the three-dimensional image of the attenuation coefficient.

In this case, the core 140 may obtain a value corresponding to three-dimensional coordinates from the three-dimensional object by using the three-dimensional linear interpolation provided by the GPU 130.

In order to use the 3D linear interpolation method, the CPU 120 may connect the gamma ray emission reconstruction image and the attenuation reconstruction image to the 3D texture memory of the GPU 130. The GPU 130 may provide a 3D linear interpolation method in Monte Carlo simulation using 3D texture memory connection. Thus, the speed of the attenuation operation can be significantly increased.

In step 680, each core 140 checks whether photons have reached the detector.

In step 680, each core 140 may calculate the photon location when the photon arrives at the detector by considering the number of detector blocks and the diameter provided by the specifications of the modeled machine.

If the photon has not reached the detector, step 620 and subsequent steps 630-680 can be performed again.

If the photons have reached the detector, step 690 is performed.

In step 690, each core 140 transmits the position of the detector to CPU 120 when the photons arrive at the detector. Alternatively, each core 140 may store the location of the photons when they arrive at the detector.

7 is a flowchart illustrating an image reconstruction method using a hybrid CPU-GPU implementation according to an embodiment of the present invention.

Steps 710, 720, 730, and 740 may correspond to steps 210, 220, 230, and 240 described above with reference to FIG. 2, respectively.

Step 720 may be performed by the CPU 120 and the GPU 130. For example, at step 720, the CPU 120 may determine the location of the photon generation within the reconstructed image. GPU 130 may be used to perform photon movement. The GPU 130 may generate information of a three-dimensional measurement line connecting two detector positions for photon pairs generated at a generation position by using Monte Carlo simulation. The information on the three-dimensional measurement line is used to mean 1) information related to the three-dimensional measurement line or 2) information used for sinogram reconstruction with respect to the information on the three-dimensional measurement line and other photon pairs generated at the location of occurrence. can do.

Detector positions may be detector coordinates when each of the two photons of the photon pair reaches the detector. Here, the information of the 3D measurement line may mean the coordinates of the 3D measurement line.

In operation 730, the CPU 120 may perform sinogram rebinding based on the information of the 3D measurement line. Here, the sinogram rebinding may mean generating one or more sinograms of the prompt sinogram and the gamma ray scattering sinogram by using the information of the 3D measurement line obtained in operation 720. It can mean part of the process of generation

Generating the information of the three-dimensional measurement line 720 and performing the sinogram reconstruction 730 may be executed in parallel by the CPU 120 and the GPU 130, respectively. Steps 720 and 730 may be separately assigned to the CPU 120 and the GPU 130, respectively. Steps 720 and 730 may each be performed repeatedly for a particular calculation. The result of the calculation in the previous iteration may not affect the calculation in the next iteration.

In operation 720, there may be a plurality of generation positions. That is, step 720 and step 730 may be performed on a plurality of photon pairs having different generation positions. In operation 720, the GPU 130 may allocate each of the plurality of generation positions to each of the plurality of threads of the GPU 130. The plurality of threads may be processed by one or more cores 140. One core may process one thread, and a plurality of threads may be simultaneously processed by one or more cores 140.

By parallel execution as described above, the entire resources of the CPU 120 and the GPU 130 can be used efficiently, and the bottleneck of memory transfer between the CPU 120 and the GPU 130 is avoided. Can be. For parallel execution, Open Multi-Processing (OpenMP) may be used for the CPU 120 and Compute Unified Device Architecth (CUDA) may be used for the GPU 120. Can be.

8 is an operation flowchart of a hybrid CPU-GPU multi scattering simulation in 3D PET according to an embodiment of the present invention.

CPU 120 may perform steps 870 and 875. In addition, the CPU 120 may perform steps 810 and 880. GPU 130 may perform steps 820, 830, 835, 840, 845, 850 and 860. Steps 870, 875, and 880 may correspond to step 730 described above. Steps 810, 820, 830, 835, 840, 845, 850, and 860 may correspond to step 720 described above. Meanwhile, step 810 may correspond to step 710.

The kernel can quickly access constant memory. Thus, geometric parameters can be stored in constant memory. Here, the geometric parameters may include a diameter, a voxel and a detector size.

In step 810, the CPU 120 may generate a table. The table may be an emission coordinate table. The emission coordinate table may include the intensity of the emission position from the emission position and the radiotracer image.

Upon generating the table, the CPU 120 may generate an initial generation position of the photons.

At step 820, GPU 130 may launch a photon pair. Here, the photon pairs may be launched toward the directions indicated by the unit vectors in the directions. Step 820 may indicate initialization performed in the GPU 130.

The GPU 130 may allocate a thread i of the GPU 130 to each occurrence position q _i . Here, i may be an integer of 1 or more and n or less. n may be the maximum number of threads that the GPU 130 can process. That is, the GPU 130 may operate one or more threads, each of which is assigned a generation location.

The allocation of the above-mentioned thread to the origination location may correspond to step 610 described above with reference to FIG. 6.

GPU 130 uses a uniform random number generator (RNG) supported by CUDA, which uses a different seed for each thread, so that the unit vector in the initial direction u ₀ Can be defined. Here, the RNG may be an xorshift RNG having a period greater than 2 ¹⁹⁰ .

Two unit vectors opposite the photon propagation can be calculated as u. Here, u = ku ₀ and k = -1, 1 may be.

Here, the emission position p may be defined as in Equation 1 below.

The calculation of each of the two unit vectors described above may correspond to step 620 described above with reference to FIG. 6.

In step 830, GPU 130 may calculate the scattering length of each of the photon pairs. Step 830 may correspond to step 630 described above with reference to FIG. 6.

Here, the scattering length may mean a length l of a direction path. The length l of the directional path may be calculated based on Equation 2 below.

는 객체의 재료(material)의 감쇄 계수일 수 있다.

는 범위 [0, 1] 내의 균등분포(uniform distribution)일 수 있다. I_A가 텍스처 메모리를 나타낼 때, 수학식 2의

는 I_A(p)로 대체될 수 있다.here,

May be attenuation coefficient of the material of the object.

Can be a uniform distribution within the range [0, 1]. When I _A represents texture memory,

May be replaced by I _A (p).

In step 835, GPU 130 may calculate the movement and absorption of each of the two photons of the photon pair. Here, the absorption may be photoelectric absorption along the ray.

By utilizing the texture memory of the GPU 130 to provide fast trilinear interpolation, line integrals that define a defined ray in the attenuation image can be calculated. Can be. The GPU 130 may update the position of each of the photon pairs by calculating the shipment. Here, the decay map in the texture memory may represent a decay image.

When the photons are located outside of the Field of View (FoV), the GPU 130 may stop the integration.

At step 840, GPU 130 may calculate Compton scattering for each of the two photons of the photon pair. Step 840 may correspond to steps 640, 650, 660, and 670 that are specialized with reference to FIG. 6.

를 정의할 수 있다.If the photons are still within the object, the GPU 130 may calculate Compton scattering based on the Klein-Nissin formula. Here, the GPU 130 uses the uniform random number [0, 1] to scatter angle

Can be defined.

In addition, the GPU 130 may inspect the energy of photons after Compton scattering. If the energy of the two photons of the photon pair is lower than the Low Level Discriminator (LLD), the GPU 130 may reject the photon pair. It can mean to break the operation on, and not use the photon pair for sinogram rebinding.

If the energy of the two photons of the photon pair is greater than or equal to the LLD, the GPU 130 may continue with the steps 830, 835, 840, 845 and 850 described above by following the propagation path of each of the two photons.

At step 680, GPU 130 may check whether the photon pair has reached the detector. Step 680 may correspond to step 680 described above with reference to FIG. 6.

Here, that the photon pair has reached the detector may mean that two photons of the photon pair are located outside of the FoV. In addition, when the photon pair reaches the detector, the energy of the two photons that reach the detector may be greater than or equal to the LDD.

If the photon pair has not reached the detector, the GPU 130 may calculate a scatter vector at step 850.

Here, the calculation of the scattering vector may mean updating the direction of the scattering flag and the unit vector u. By using the scattering flag, the scattering event of the Compton scattering can be recorded.

After step 850 is performed, step 830 may be repeated.

If the photon pair has reached the detector, step 860 may be performed.

In step 860, after the photon pair reaches the detector, GPU 130 may convert the two detector positions into LoR. GPU 130 may calculate two detector locations and may calculate detector cross-sections. LoR can be calculated by converting two detector positions to LoR.

The information of the three-dimensional measurement line described above may include LoR, two detector positions and detector cross-sections.

In step 870, the CPU 120 may perform the conversion to the downsampled sinogram. Here, the conversion to the downsampled sinogram may include a sinogram liby corresponding to the movement of both. Can mean Ning. The CPU 120 may use the information calculated in step 860 for the sinogram revival.

In step 875, the CPU 120 may check whether it is finished. If the calculated photon pairs remain, step 820 may be performed again. If the movement of all photon pairs has been calculated, step 880 may be performed.

In step 880, the CPU 120 may upsample the scattering sinogram.

The CPU 120 may calculate a scaling factor for upsampling the scattering sinogram and obtain the scattering sinogram.

After quantum shift computed by GPU 130 and sinogram rebinding by CPU 120, CPU 120 may upsample the sinogram to the original sinogram size.

GPU 130 may be used for the upsampling. For example, for fast upscaling, fast linear interpolation using texture memory supported by GPU 130 may be used. Spatial smoothing using text memory can be computationally efficient thanks to large virtual detector arrays.

Finally, the CPU 120 can calculate the minimum square fitting by using a simulated prompt from the MSS, and measure the calculated factor to calculate the scaling factor. The nogram can be calculated. In addition, prompt sinograms can be directly compared to the scatter portion of the measurement.

Code 1 below represents a pseudo code executed in CPU 120 for the hybrid CPU-GPU multiscattering simulation described above.

[Code 1]

Code 2 below represents pseudo code executed on GPU 130 for the hybrid CPU-GPU multiscattering simulation described above.

[Code 2]

9 illustrates the sequential performance of multiple tasks in accordance with an example of the present invention.

The total number of threads in GPU 130 is limited. Thus, all tasks may not be assigned to a single execution.

Tasks for multiple scattering simulations can be divided into one or more tasks. Each of the one or more tasks may be continuously scheduled for hybrid CPU and GPU computations.

The GPU 130 and the CPU 120 may each sequentially process n tasks. Here, n may be an integer of 1 or more.

In order to make the most of the CPU 120 and the GPU 130, two jobs may be separated. The two jobs may be represented as a task to perform photon shift and a task to perform sinogram transformation, respectively.

In FIG. 9, n tasks for performing photon shifting processed by the GPU 130 are illustrated, and n tasks for performing sinogram transformation processed by the CPU 120 are illustrated. In one iteration, GPU 130 and CPU 120 may each process one task. In FIG. 9, the space between the dotted lines may represent one iteration. The result of calculating the task of the GPU 130 at a particular iteration may be used in the task of the CPU 120 at the next iteration.

The plurality of occurrence positions of step 720 may be divided into a plurality of groups, and each of the plurality of groups generates information of a three-dimensional measurement line (720) and performs sinogram rebinding. It can be processed by iterations of 730. The plurality of groups may be processed sequentially. Each of the plurality of tasks may process one group of the plurality of groups.

That is, a plurality of occurrence positions may be divided into n groups, and n groups may be processed by n + 1 repetitions. The information of the three-dimensional measurement lines for the generated positions processed at the k-th iteration can be used in sinogram rebinding at the k + 1 iterations. Here, k may be an integer of 1 or more and n + 1 or less.

For concurrent executions between CPU 120 and GPU 130, event query and asynchronous memory transfer can be used. GPU 130 may record the events at any point in the program and may query the events when the above events are completed. When a synchronous function is called, control can wait for the function to complete. On the other hand, when an asynchronous function is called, control can continue regardless of the completion of the function. Thus, asynchronous memory transfer can be used.

In particular, after execution of the task of the GPU 130, an asynchronous memory transfer function can be called, and the task of the next GPU 130 can be started regardless of the completion of the transfer. Information of the 3D measurement line may be transmitted from the GPU 130 to the CPU 120 through an asynchronous memory transfer. The GPU 130 may wait for completion of the transmission by constantly checking the event query.

While the GPU 130 is performing quantum movement, the memory transfer may be completed. The CPU 120 may calculate the sinogram rebinding corresponding to the quantum shift using OpenMP.

The above-described architecture can maximize the utilization of the CPU 120 and can greatly improve the overall performance in gamma ray scattering simulation.

Method according to an embodiment of the present invention is implemented in the form of program instructions that can be executed by various computer means may be recorded on a computer readable medium. The computer readable medium may include program instructions, data files, data structures, etc. alone or in combination. The program instructions recorded on the medium may be those specially designed and constructed for the present invention or may be available to those skilled in the art of computer software. Examples of computer-readable recording media include magnetic media such as hard disks, floppy disks, and magnetic tape, optical media such as CD-ROMs, DVDs, and magnetic disks, such as floppy disks. Magneto-optical media, and hardware devices specifically configured to store and execute program instructions, such as ROM, RAM, flash memory, and the like. Examples of program instructions include not only machine code generated by a compiler, but also high-level language code that can be executed by a computer using an interpreter or the like. The hardware device described above may be configured to operate as one or more software modules to perform the operations of the present invention, and vice versa.

As described above, the present invention has been described by way of limited embodiments and drawings, but the present invention is not limited to the above embodiments, and those skilled in the art to which the present invention pertains various modifications and variations from such descriptions. This is possible.

Therefore, the scope of the present invention should not be limited to the described embodiments, but should be determined by the equivalents of the claims, as well as the claims.

100: video reconstruction device
110:
120: CPU
130: GPU
140: core

Claims (18)

Receiving an actual measured sinogram generated by measuring an object;
Performing a Monte Carlo simulation on the initial occurrence locations of gamma rays in the object, wherein the Monte Carlo simulation is performed in parallel for each of the initial generation locations of photons by using a GPU;
Generating a gamma ray scattering sinogram with gamma ray scattering corrected based on information on scattered photons generated by the Monte Carlo simulation; And
Generating a reconstructed image based on the gamma ray scattering sinogram
Image reconstruction method comprising a. Determining a location of generation of photons in the reconstructed image;
Calculating the movement of the two photons generated at the location of occurrence by using a Monte Carlo simulation;
Calculating two detector positions and values changed during the movement, wherein the two detector positions are detector coordinates when each of the two photons reaches the detector, and the values changed during the movement are each of the two photons Values changed during the passage and movement of this object;
Calculating a three-dimensional measurement line connecting the two detector positions;
Generating a prompt sinogram and a gamma ray scattering sinogram based on the three-dimensional measurement line; And
Scaling the gamma ray scattering sinogram
Including, gamma-ray scattering sinogram generation method. The method of claim 2,
The generation position is a plurality,
Each of the plurality of occurrence positions is assigned to one core of the plurality of cores of the GPU,
Calculating the movement of the two photons generated at the generation position by using the Monte Carlo simulation and calculating the two detector positions and the values changed during the movement are performed in parallel by the plurality of cores Scattering sinogram generation method. The method of claim 3,
Determining a generation position of photons in the reconstructed image,
Considering positions having values within the reconstructed image as the plurality of occurrence positions; And
Generating a table including the plurality of occurrence locations
Lt; / RTI >
Wherein each of the plurality of cores of the GPU is assigned one or more occurrence locations determined based on the number of cores in the table. The method of claim 2,
Calculating the movement of two photons generated at the generation position by using the Monte Carlo simulation,
A core of a GPU being allocated the generation position;
Randomly generating, by the core, an initial direction of movement at the assigned occurrence position of photons;
The core randomly generating a moving distance of the photons;
Calculating a position at which the core moves the photon by the movement distance in the initial movement direction;
Calculating attenuation of the photons as the core moves to the moved position; And
Confirming that the photon has reached the detector;
Including, gamma-ray scattering sinogram generation method. The method of claim 5,
Computing the attenuation of the photon according to the movement of the core to the moved position,
The attenuation of the photons is calculated by obtaining a value corresponding to the three-dimensional coordinates of the moved position of the core in the three-dimensional image of the attenuation coefficient,
And the core is obtained by reading a value corresponding to the three-dimensional coordinates by using a three-dimensional linear interpolation method provided by the GPU. The method of claim 5,
Calculating the movement of two photons generated at the generation position by using the Monte Carlo simulation,
Determining whether scattering of the photons occurs at the moved position of the core;
When the scattering occurs, the core randomly generating scattering angles of the scattering; And
If the scattering occurs, the core calculating a Compton scattering value of the scattering
Further comprising, gamma scattering sinogram generation method. The method of claim 7, wherein
Determining whether or not scattering of the photons in the moved position of the core,
Determining that the scattering has not occurred if the attenuation coefficient of the moved position of the core is less than or equal to a boundary value, and determining that the scattering has occurred if the attenuation coefficient of the moved position is greater than the boundary value of the core;
Including, gamma-ray scattering sinogram generation method. The method of claim 2,
Generating a prompt sinogram and gamma ray scattering sinogram based on the three-dimensional measurement line,
Converting the three-dimensional measurement line into a sinogram format;
Storing the 3D measurement line converted into the sinogram format in a prompt sinogram; And
When the two photons are scattered photons, storing the three-dimensional measurement line converted into the sinogram format in the gamma ray scattering sinogram
Including, gamma-ray scattering sinogram generation method. 10. The method of claim 9,
Converting the three-dimensional measurement line in the sinogram format,
Dividing the three-dimensional measurement line into a vector in an XY direction and a Z direction; And
Converting the three-dimensional measurement line into the sinogram format based on the vector
Including, gamma-ray scattering sinogram generation method. The method of claim 2,
Adjusting the scale of the gamma ray scattering sinogram,
Calculating a scale constant that minimizes the error between the actual measured sinogram and the prompt sinogram; And
Scaling the gamma ray scattering sinogram by multiplying the scale constant by the gamma ray scattering sinogram
Including, gamma-ray scattering sinogram generation method. A receiver configured to receive an actual measured sinogram generated by measuring an object;
A central processing unit (CPU) for generating a first reconstructed image using the actually measured sinogram; And
A graphics processing unit (GPU) for performing a Monte Carlo simulation on initial generation positions at which gamma rays are generated in the object; The GPU performs the Monte Carlo simulation in parallel by initial occurrence positions of photons by using one or more cores.
/ RTI >
The CPU generates a gamma ray-corrected sinogram based on the information about scattered photons generated by the Monte Carlo simulation, and generates a reconstructed image based on the gamma ray-corrected sinogram. Image reconstruction device. Determine a generation location of photons in the reconstructed image, calculate a three dimensional measurement line connecting two detector positions, generate a prompt sinogram and a gamma ray scattering sinogram based on the three dimensional measurement line, A CPU for scaling the gamma ray scattering sinogram; And
A GPU that calculates the movement of the two photons generated at the generation position by using a Monte Carlo simulation line and calculates the two detector positions and the values changed during the movement-the two detector positions are each Detector coordinates upon reaching the detector, the values changed during the movement being the values changed while each of the two photons passes through the object and moves;
Including, gamma-ray scattering sinogram generating device. The method of claim 13,
The generation position is a plurality,
For two photons generated at each of the plurality of generation positions, the plurality of cores of the GPU calculates the movement of two photons generated at each of the generation positions in parallel by using Monte Carlo simulation, and the two detectors An apparatus for generating gamma ray scattering sinograms that calculates the positions and values changed during movement in parallel. Determining a generation position of photons in the reconstructed image by the CPU;
GPU generating information of a three-dimensional measurement line connecting two detector positions for the photon pair generated at the generation position by using Monte Carlo simulation, wherein the detector positions are such that each of the two photons of the photon pair will reach the detector. Detector coordinates when; And
The CPU performing sinogram revision based on the information of the 3D measurement line;
Lt; / RTI >
The generation position is plural,
The GPU allocates each of the plurality of occurrence positions to each of the plurality of threads of the GPU,
Generating the information of the three-dimensional measurement line and performing the sinogram rebinding are executed in parallel by the CPU and the GPU, respectively. 16. The method of claim 15,
The plurality of occurrence positions are divided into a plurality of groups,
And said plurality of groups are each processed by iterations of generating information of said three-dimensional measurement line and performing said sinogram rebinding. 16. The method of claim 15,
And information of the three-dimensional measurement line is transmitted from the GPU to the CPU through an asynchronous memory transfer. A CPU for performing sinogram revision based on information of a three-dimensional measurement line, which determines a location of generation of photons in the reconstructed image; And
GPU using Monte Carlo simulation to generate information of the three-dimensional measurement line connecting two detector positions for the photon pair generated at the generation position
Lt; / RTI >
The detector positions are detector coordinates when each of the two photons of the photon pair reaches the detector,
The generation position is plural,
The GPU allocates each of the plurality of occurrence positions to each of the plurality of threads of the GPU,
And generating the information of the three-dimensional measurement line and performing the sinogram rebinding are executed in parallel by the CPU and the GPU, respectively.

Images