WO2011027402A1 - 核医学用データ処理方法および核医学診断装置 - Google Patents

核医学用データ処理方法および核医学診断装置 Download PDF

Info

Publication number
WO2011027402A1
WO2011027402A1 PCT/JP2009/004381 JP2009004381W WO2011027402A1 WO 2011027402 A1 WO2011027402 A1 WO 2011027402A1 JP 2009004381 W JP2009004381 W JP 2009004381W WO 2011027402 A1 WO2011027402 A1 WO 2011027402A1
Authority
WO
WIPO (PCT)
Prior art keywords
nuclear medicine
parallel
cross
lor
section
Prior art date
Application number
PCT/JP2009/004381
Other languages
English (en)
French (fr)
Inventor
赤澤礼子
北村圭司
山田賢志
Original Assignee
株式会社島津製作所
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 株式会社島津製作所 filed Critical 株式会社島津製作所
Priority to JP2011529700A priority Critical patent/JP5263402B2/ja
Priority to CN200980161257.2A priority patent/CN102483459B/zh
Priority to US13/388,792 priority patent/US8987674B2/en
Priority to PCT/JP2009/004381 priority patent/WO2011027402A1/ja
Publication of WO2011027402A1 publication Critical patent/WO2011027402A1/ja

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01T—MEASUREMENT OF NUCLEAR OR X-RADIATION
    • G01T1/00—Measuring X-radiation, gamma radiation, corpuscular radiation, or cosmic radiation
    • G01T1/29—Measurement performed on radiation beams, e.g. position or section of the beam; Measurement of spatial distribution of radiation
    • G01T1/2914—Measurement of spatial distribution of radiation
    • G01T1/2985—In depth localisation, e.g. using positron emitters; Tomographic imaging (longitudinal and transverse section imaging; apparatus for radiation diagnosis sequentially in different planes, steroscopic radiation diagnosis)
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus for radiation diagnosis, e.g. combined with radiation therapy equipment
    • A61B6/02—Devices for diagnosis sequentially in different planes; Stereoscopic radiation diagnosis
    • A61B6/03—Computerised tomographs
    • A61B6/037—Emission tomography

Definitions

  • the present invention relates to a nuclear medicine data processing method and a nuclear medicine diagnostic apparatus for performing forward projection and backprojection operations on list data generated from event data obtained by detecting radiation.
  • the above-described nuclear medicine diagnosis apparatus that is, an ECT (Emission Computed Tomography) apparatus will be described by taking a PET (Positron Emission Tomography) apparatus as an example.
  • the PET device detects a plurality of photons generated by the annihilation of positrons (Positrons), and only detects ⁇ rays simultaneously with a plurality of detectors (that is, only when they are counted simultaneously). It is configured to reconstruct an image.
  • the process of drug accumulation in the target tissue is measured over time, whereby quantitative measurement of various biological functions is possible. Therefore, the image obtained by the PET apparatus has function information.
  • a positron (positron) radioactive isotope for example, 15 O, 18 F, 11 C, etc.
  • a detector array comprising a number of ⁇ -ray detectors arranged in a ring shape so as to surround the body axis that is the longitudinal axis of the subject.
  • calculation is performed by a computer in the same manner as in ordinary X-ray CT (Computed Tomography), and the image is specified in the plane, and an image of the subject is created.
  • Non-Patent Document 1 is an indispensable technique when an image is reconstructed in a nuclear medicine diagnostic apparatus such as a PET apparatus.
  • an event for detecting ⁇ rays is referred to as an “event”
  • data that is simultaneously counted is referred to as “coincidence data”
  • single event data data that is not simultaneously counted
  • Data for each event transferred from the gantry of the PET apparatus is called “list data”.
  • This list data includes position information such as X coordinate, Y coordinate, and Z coordinate for each event, time information, and the like. These pieces of information are arranged in time series.
  • the list data is converted into a histogram to obtain histogram data.
  • This histogram data is data integrated at a predetermined time.
  • this histogram data there is a sinogram with the vertical axis representing the projection direction and the horizontal axis representing the pixel.
  • Non-Patent Document 2 a mode using list data
  • list mode the data increases by the number of events (events), and the large amount of calculation becomes a problem.
  • the amount of calculation is further enormous.
  • MIMD Multiple Instruction-Multiple Data
  • MPI Message Passing Interface
  • SIMD type operations can simultaneously process a plurality of data with a single instruction, they are suitable for performing the same operation processing on a large number of operations.
  • parallel operations such as forward projection processing and back projection processing using a GPU
  • writing to the same memory may occur (so-called “memory contention”).
  • memory contention Ingenuity is necessary.
  • Non-Patent Document 3 uses a sinogram. In the sinogram, data of the same pixel is duplicated, resulting in memory contention. Then, the parallel study of the successive approximation method using list data is also performed as a parallel operation using the same GPU (for example, refer nonpatent literature 4).
  • Non-Patent Document 4 uses a mechanism peculiar to GPU, and has a problem that it is difficult to apply to other SIMD type mechanisms other than GPU. Since the GPU is originally for image processing, the use of a mechanism unique to the GPU means that data is treated as a texture image and a pixel (pixel) shader or a vertex shader is used.
  • the present invention has been made in view of such circumstances, and an object of the present invention is to provide a nuclear medicine data processing method and a nuclear medicine diagnostic apparatus that are highly versatile and capable of speeding up operations. To do.
  • the nuclear medicine data processing method performs forward projection to obtain post-projection-processed data by performing forward projection operation in parallel for each LOR on list data generated from event data obtained by detecting radiation.
  • the forward projection processing step performs forward projection operation in parallel for each LOR (Line Of ⁇ ⁇ ⁇ Response) on list data generated from event data obtained by detecting radiation.
  • LOR Line Of ⁇ ⁇ ⁇ Response
  • backprojection calculation is performed in parallel for each area obtained by dividing the cross section of the image space. Since list processing is performed in parallel for each LOR for each LOR, for example, forward projection and backprojection calculations can be realized without using a GPU-specific mechanism, and can be applied to other calculation mechanisms. Is. Further, parallel processing is performed on the list data, and the parallel processing is performed by forward projection calculation and back projection calculation, respectively, so that memory contention can be prevented and high speed can be realized.
  • the back projection processing step is performed for each region obtained by dividing the cross section of the image space. Back projection operations are performed in parallel. As a result, the versatility is high and the calculation speed can be increased.
  • the main direction determination step determines the main direction of the LOR based on the LOR vector and the arrangement direction of the detector for detecting radiation.
  • the cross section determining step determines a cross section of the image space through which the LOR passes based on the main direction determined in the main direction determining step.
  • the calculation area determination step determines the calculation area where the LOR in the cross section traverses. More specifically, the forward projection processing step and the back projection processing step perform the following arithmetic processing.
  • the forward projection processing step calculates in parallel for each calculation region
  • the back projection processing step calculates the post-projection processing data or list data in parallel for each region obtained by dividing the calculation region into a plurality of regions. . In this way, determining the main direction makes it easy to classify the main direction and perform each calculation process.
  • the back projection process step divides the cross section into a plurality of regions, and the divided and dispersed regions belong to one region. It is preferable to perform operations in parallel for each area to which the user belongs. By classifying the regions in this way, it is possible to prevent localization in a specific thread and improve the calculation efficiency.
  • “Thread” indicates a unit of parallel data processing.
  • a more preferable example is to associate a plurality of row regions separated by a predetermined interval as one region described above. In this way, by classifying the regions evenly distributed, it is possible to further prevent localization in a specific thread.
  • An example of the nuclear medicine data processing method of these inventions described above is to use a sequential approximation method in which the forward projection process step and the backprojection process step are executed a plurality of times to sequentially approximate and update the image.
  • image reconstruction is performed using the successive approximation method as described above. Therefore, the present invention may also be applied to the successive approximation method.
  • filtered back projection FBP: “Filtered” Back Projection
  • FBP Filtered” Back Projection
  • the nuclear medicine diagnosis apparatus provides forward projection processing means for obtaining forward-projected data by performing forward projection operation on list data generated from event data obtained by detecting radiation in parallel for each LOR. And backprojection processing means for performing a backprojection operation in parallel for each area obtained by dividing the cross section of the image space.
  • the forward projection processing means performs forward projection operation on the list data generated from the event data obtained by detecting radiation in parallel for each LOR, and performs forward projection processed data. Since the backprojection processing means performs the backprojection operation on the data or list data after the forward projection processing in parallel for each region obtained by dividing the cross section of the image space, the versatility is similar to the nuclear medicine data processing method. It is high and the calculation speed can be increased.
  • the forward projection processing is performed by performing the forward projection operation on the list data generated from the event data obtained by detecting the radiation in parallel for each LOR. Since post-projection data is obtained and backprojection calculation is performed in parallel for each area obtained by dividing the cross section of the image space, the post-projection processing data or list data is highly versatile and the calculation speed can be increased.
  • PET Positron Emission Tomography
  • FIG. 6 is a schematic diagram showing a flow of list data in a series of steps S1 to S5.
  • FIG. 1 is a side view and a block diagram of a PET (Positron Emission Tomography) apparatus according to an embodiment
  • FIG. 2 is a schematic perspective view of a ⁇ -ray detector.
  • the PET apparatus includes a top plate 1 on which a subject M is placed as shown in FIG.
  • the top plate 1 is configured to move up and down and translate along the body axis Z of the subject M.
  • the subject M placed on the top 1 is scanned from the head to the abdomen and foot sequentially through the opening 2a of the gantry 2, which will be described later. Get the image. Note that there is no particular limitation on the scanned part and the scanning order of each part.
  • the PET apparatus includes a gantry 2 having an opening 2a and a ⁇ -ray detector 3.
  • the ⁇ -ray detector 3 is arranged in a ring shape so as to surround the body axis Z of the subject M, and is embedded in the gantry 2.
  • the ⁇ -ray detector 3 corresponds to the detector in the present invention.
  • the PET apparatus includes a top plate driving unit 4, a controller 5, an input unit 6, an output unit 7, a memory unit 8, a coincidence circuit 9, and an arithmetic processing unit 10.
  • the top plate driving unit 6 is a mechanism for driving the top plate 1 so as to perform the above-described movement, and is configured by a motor or the like not shown.
  • the arithmetic processing unit 10 corresponds to the forward projection processing means and back projection processing means in the present invention.
  • the controller 5 comprehensively controls each part constituting the PET apparatus according to the present embodiment.
  • the controller 5 and the arithmetic processing unit 10 are configured by a central processing unit (CPU) or the like.
  • the arithmetic processing unit 10 is configured by a SIMD type mechanism represented by a GPU or the like.
  • the arithmetic processing unit 10 may be a SIMD type mechanism other than the GPU, and is not particularly limited as long as parallel arithmetic is possible.
  • the input unit 6 sends data and commands input by the operator to the controller 5.
  • the input unit 6 includes a pointing device represented by a mouse, a keyboard, a joystick, a trackball, a touch panel, and the like.
  • the output unit 7 includes a display unit represented by a monitor, a printer, and the like.
  • the memory unit 8 includes a storage medium represented by ROM (Read-only Memory), RAM (Random-Access Memory), and the like.
  • the count value (count) simultaneously counted by the coincidence circuit 9 the data relating to the coincidence counting such as the detector pair consisting of the two ⁇ -ray detectors 3 and the LOR, and the arithmetic processing unit 10 perform the calculation.
  • Various processed data are written and stored in the RAM, and are read from the RAM as necessary.
  • the ROM stores in advance a program for performing imaging including various types of nuclear medicine diagnosis, and the controller 5 and the arithmetic processing unit 10 execute the program to perform nuclear medicine diagnosis according to the program. Do each.
  • a program related to CUDA is stored in advance in a ROM so that a parallel computing architecture called “CUDA (provided by NVida)” can execute a parallel operation using a GPU.
  • the arithmetic processing unit 10 executes steps S1 to S5 described later.
  • the ⁇ -rays generated from the subject M to which the radiopharmaceutical is administered are converted into light by the scintillator block 31 (see FIG. 2) of the ⁇ -ray detector 3, and the converted light is photoelectron of the ⁇ -ray detector 3.
  • a multiplier tube (PMT: Photo Multiplier Tube) 33 (see FIG. 2) multiplies and converts it into an electrical signal. The electric signal is sent to the coincidence circuit 9 as an event.
  • the coincidence circuit 9 checks the position of the scintillator block 31 (see FIG. 2) and the incident timing of the ⁇ rays, and only when the ⁇ rays are simultaneously incident on the two scintillator blocks 31 on both sides of the subject M. The sent event is determined to be appropriate data.
  • the coincidence counting circuit 9 rejects. That is, the coincidence counting circuit 9 detects that ⁇ rays are simultaneously observed in the two ⁇ ray detectors 3 based on the above-described electrical signal.
  • the event sent to the coincidence circuit 9 is sent to the arithmetic processing unit 10.
  • the arithmetic processing unit 10 performs image reconstruction by forward projection processing or back projection processing to obtain an image of the subject M.
  • the image is sent to the output unit 7 via the controller 5. In this manner, nuclear medicine diagnosis is performed based on the image obtained by the arithmetic processing unit 10. Specific functions of the arithmetic processing unit 10 will be described later.
  • the ⁇ -ray detector 3 includes a scintillator block 31, a light guide 32 optically coupled to the scintillator block 31, and photoelectrons optically coupled to the light guide 32.
  • a multiplier (hereinafter simply abbreviated as “PMT”) 33 is provided.
  • Each scintillator element constituting the scintillator block 31 converts ⁇ rays into light by emitting light with the incidence of ⁇ rays. By this conversion, the scintillator element detects ⁇ rays.
  • Light emitted from the scintillator element is sufficiently diffused by the scintillator block 31 and input to the PMT 33 via the light guide 32.
  • the PMT 33 multiplies the light converted by the scintillator block 31 and converts it into an electric signal. The electric signal is sent to the coincidence circuit 9 (see FIG. 1) as an event as described above.
  • FIG. 3 is a flowchart showing a flow of a series of list data processing methods
  • FIG. 4 is a schematic diagram showing coincidence with a ⁇ -ray detector for explaining detection probabilities
  • FIG. 6 is a schematic diagram that defines the relationship between the arrangement of ⁇ -ray detectors arranged in a three-dimensional manner and three-dimensional coordinates.
  • FIG. 6 is a schematic diagram in which the LOR direction vector is shown together with the conceptual diagram of projection
  • FIG. FIG. 8 is a schematic diagram of a calculation region in a cross section of the image space through which the LOR passes, FIG.
  • FIG. 8 is a schematic diagram for explaining the parallel calculation of the LOR and projection values for one subset
  • FIG. 9 is a diagram through which the LOR passes.
  • FIG. 10 is a schematic diagram showing division of the calculation area in the cross section of the image space
  • FIG. 10 is a schematic diagram showing the flow of list data in a series of steps S1 to S5. 4 and 5, only the scintillator block 31 is illustrated as the ⁇ -ray detector 3, and the light guide 32 and the PMT 33 are not illustrated.
  • a description will be given by taking as an example a pixel composed of three-dimensional voxels.
  • the list-mode 3D-DRAMA method (Dynamic-Row-Action-Maximum-Likelihood-Algorithm) using list data is adopted as the successive approximation algorithm, and corrections such as absorption correction, sensitivity correction, and random correction are incorporated.
  • parallel operation using a GPU is implemented by the above-mentioned CUDA.
  • the update formula of the list-mode 3D-DRAMA method is expressed by the following formula (1).
  • i (t) in the above equation (1) indicates the t-th coincidence event.
  • l (lowercase L) is the number of subsets
  • a i (t) j is a probability that ⁇ -ray photons emitted from voxel j are detected by LOR i (t) , and is also referred to as a “system matrix”.
  • the detector response function is also taken into consideration.
  • C j represents a normalization matrix
  • p lj represents a blocking factor.
  • Non-Patent Document 2 Two types of blocking factors have been introduced in Non-Patent Document 2 described above, but in this study, convergence is high and a method of calculating using all list data is employed.
  • ⁇ (k, l) is a relaxation coefficient and is related to the convergence speed.
  • ⁇ and ⁇ 0 are parameters for controlling ⁇ (k, l) .
  • a ′ i (t) j is a i (t) j after absorption correction, and the correction-related definition is shown in the following equation (2).
  • a i (t) j represents an absorption correction coefficient
  • ⁇ j represents a distribution of absorption coefficients.
  • r i is a value obtained by predicting the random count rate from the count rates s i0 and s i1 of the single count and the time window 2 ⁇ of the coincidence count.
  • FP in the above formula (1) represents Forward Projection, that is, forward projection, and calculation in the forward projection processing is performed in the portion in the above formula (1).
  • BP in the above equation (1) represents Back Projection, that is, back projection, and the calculation in the back projection process is performed in the portion in the above equation (1).
  • the initial image x j (0, 0) is set appropriately.
  • the initial image x j (0,0) may be an image having a uniform pixel value, for example, and x j (0,0) > 0.
  • x j ( 0,0) and a ′ i (t) j corrected by the above equation (2) are used, and repeatedly substituting it into the above equation (1), x j ( 0,0), ..., x j ( 0, L-1) is determined sequentially, x finally the obtained x j of (0, L-1) by the x j (1, 0) Raise to j (1,0) .
  • x j is sequentially incremented (x j (0,0) , x j (1,0) ..., X j (k, 0) ).
  • the number of times k representing repetition is not particularly limited, and may be set as appropriate.
  • Step S1 Determination of Main Direction
  • the horizontal axis is x
  • the vertical axis is y
  • the body axis of the subject M is z.
  • the body axis z is an axis in the depth direction.
  • the direction vector parallel to the LOR is defined as d (d x , d y , d z ) (where
  • 1, unit vector), and the LOR vector d and ⁇ -ray detection Based on the arrangement direction of the vessel 3, the main direction of the LOR is determined.
  • FIG. 6 is a schematic diagram when the x direction is the main direction. This step S1 corresponds to the main direction determining step in the present invention.
  • Step S2 Cross Section Determination
  • the cross section of the image space through which the LOR passes is determined based on the main direction determined in step S1.
  • the plane in which the main directions are orthogonal is determined as the cross section of the image space through which the LOR passes.
  • the cross section is the yz plane as shown in FIGS.
  • This step S2 corresponds to the cross section determining step in the present invention.
  • Step S3 Calculation Area Determination Based on the cross section determined in step S2 (the yz plane in FIG. 6), the calculation area crossed by the LOR in the cross section is determined.
  • the cross section that is the start point of the LOR (denoted by “slice_start” in FIG. 6) and the cross section that is the end point of the LOR (denoted by “slice_end” in FIG. 6) are determined.
  • the coordinates on the cross section In FIGS. 6 and 7, since the cross section is the yz plane, the coordinates on the cross section that is the start point are y and z coordinates, and the coordinates on the cross section that is also the end point are y and z coordinates.
  • one of the y-coordinates of the start and end points is the minimum coordinate (indicated as “Min_0” in FIG. 7), and the other is the maximum coordinate (in FIG. 7, “Max_0”). ).
  • one of the z coordinates of the start point and the end point is the minimum coordinate (indicated as “Min_1” in FIG. 7), and the other is the maximum coordinate (indicated as “Max_1” in FIG. 7). That is, on the cross section, the LOR crosses the maximum coordinate position from the minimum coordinate position on each cross section. Therefore, the calculation area traversed by the LOR is an area delimited by Min_0, Max_0, Min_1, and Max_1 as shown in FIG.
  • This step S3 corresponds to the calculation area determination step in the present invention.
  • Step S4 Forward Projection Processing
  • the list data is subjected to calculation processing in parallel in the calculation area determined in step S3.
  • the system matrix is obtained for each voxel j, and the projection value is obtained by multiplying and adding the pixel value in the portion indicated by FP in the above expression (1) within the range of the determined calculation area.
  • the parallel operation in the above equation (1) is performed for each LOR.
  • the memory into which the projection value is substituted is configured to write one element per event per subset, for example, as shown in FIG. In FIG. 8, “Thread_0”, “Thread_1”, “Thread_2”, “Thread_3”,... Indicate processing units of parallel data to be written.
  • Each LOR has a one-to-one correspondence with each thread per subset.
  • This step S4 corresponds to the forward projection processing step in this invention.
  • Step S5 Backprojection processing (region division) Back projection processing is performed on the post-forward projection processing data or list data on which forward projection processing has been performed in step S4.
  • each cross-section of the image space through which the LOR passes is divided into a plurality of regions, and the divided and dispersed regions belong to one region, and are parallel to each of the belonging regions.
  • the backprojection value is updated by adding in the part indicated by BP in the above-described equation (1) in each of the areas to which it belongs.
  • row regions that are separated by a predetermined interval belong as one region.
  • the calculation area is a cross section (slice) of 256 ⁇ 256 pixels
  • FIG. 9B is an enlarged view of FIG.
  • it is defined as one area of every 7th row, it is divided into 8 regions per slice. That is, it is divided by “Thread_0”, “Thread_1”,... “Thread_7” shown in FIG.
  • This step S5 corresponds to the back projection process step in this invention.
  • list data for each event (ie, for each LOR) (indicated by “Event_0”, “Event_1”, “Event_2”,... In FIG. 10) is classified in the main direction by determining the main direction in step S1.
  • the main direction (denoted by in FIG. 10, "the main direction x list data") list data in the x-direction
  • the main direction is classified into list data in the y direction (represented as “list data in the main direction y” in FIG. 10), and in each step S5, the back projection process is performed in parallel.
  • the successive approximation method described above is not limited to the above-described DRAMA method, and may be a static (that is, static) RAMLA method (Row-Action Maximum Likelihood Algorithm) or an ML-EM method (Maximum Likelihood Expectation Maximization). ) Or OSEM method (Ordered Subset ML-EM).
  • step S4 forward projection process
  • step S5 back projection processing
  • the data after the forward projection processing or the list data is back-projected in parallel for each region obtained by dividing the cross section of the image space. Since list processing is performed in parallel for each LOR for each LOR, for example, forward projection and backprojection calculations can be realized without using a GPU-specific mechanism, and can be applied to other calculation mechanisms. Is.
  • step S5 back projection processing divides the cross section of the image space. Back projection operations are performed in parallel for each region. As a result, the versatility is high and the calculation speed can be increased.
  • step S1 main direction determination
  • step S2 determines the cross section of the image space through which the LOR passes based on the main direction determined in step S1 (determining the main direction).
  • step S3 determines the calculation area where the LOR in the cross section traverses.
  • step S4 forward projection process
  • step S5 back projection process
  • step S4 forward projection processing
  • step S5 back projection processing
  • step S4 forward projection processing
  • step S5 back projection processing
  • the above-described step S5 (back projection process) is performed by dividing the above-described cross section into a plurality of regions, and the divided and dispersed regions belong to one region.
  • the calculation is preferably performed in parallel for each area. By classifying the regions in this way, it is possible to prevent localization in a specific thread and improve the calculation efficiency.
  • each row region separated by a predetermined interval belongs as one region (see FIG. 9). In this way, by classifying the regions evenly distributed, it is possible to further prevent localization in a specific thread.
  • step S4 forward projection process
  • step S5 back projection process
  • step S4 forward projection process
  • step S5 back projection process
  • the PET apparatus obtains post-projection-processed data by performing the forward projection operation in parallel for each LOR on the list data generated from the event data obtained by detecting ⁇ , and the forward projection.
  • an arithmetic processing unit 10 that performs the arithmetic processing of the above-described steps S1 to S5 that performs backprojection calculation in parallel for each region obtained by dividing the cross section of the image space.
  • the present invention is not limited to the above embodiment, and can be modified as follows.
  • step S4 forward projection processing
  • step S5 back projection processing
  • the Feldkamp method using filtered back projection may be applied.
  • the main direction is classified by determining the main direction, but it is not always necessary to determine the main direction. However, in order to facilitate each calculation process, it is more preferable to classify the main direction by determining the main direction as in the embodiment.
  • the main direction is determined after setting the horizontal axis to x, the vertical axis to y, and the body axis of the subject M to z, but the main direction is not necessarily limited to these directions.
  • the main direction is determined after arranging the ⁇ -ray detectors 3 arranged in a ring shape on the xy plane as shown in FIG. 5, but depending on the arrangement state of the ⁇ -ray detectors.
  • the main direction may be determined as appropriate.
  • the main direction may be determined based on the planes on the two detectors to be counted simultaneously.
  • the larger direction is set as the main direction when determining the main direction, but the present invention is not necessarily limited to this. If the main direction is determined based on a certain rule, for example, the smaller direction may be set as the main direction.
  • the plane in which the main direction is orthogonal is determined as the cross section of the image space through which the LOR passes.
  • the plane in which the main direction is orthogonal is not necessarily determined as the cross section.
  • a plane crossing at an angle other than 90 ° may be determined as the cross section.
  • each row area separated by a predetermined interval is assigned as one area, but is divided into block areas other than the row area, and some of them are collectively assigned as one area. You may let them.
  • the nuclear medicine diagnosis apparatus has been described.
  • the present invention can be applied to a gamma camera having a collimator mechanism and SPECT.

Landscapes

  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • High Energy & Nuclear Physics (AREA)
  • Molecular Biology (AREA)
  • Physics & Mathematics (AREA)
  • Medical Informatics (AREA)
  • Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Biophysics (AREA)
  • Optics & Photonics (AREA)
  • Pathology (AREA)
  • Radiology & Medical Imaging (AREA)
  • Biomedical Technology (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Surgery (AREA)
  • Animal Behavior & Ethology (AREA)
  • General Health & Medical Sciences (AREA)
  • Public Health (AREA)
  • Veterinary Medicine (AREA)
  • Nuclear Medicine (AREA)

Abstract

　この発明のＰＥＴ装置および処理方法では、ＬＯＲ毎に（γ線を検出することにより得られた事象データから生成された）リストデータについて並列にステップＳ４（順投影），Ｓ５（逆投影）の演算処理を行っているので、他の演算機構に適用することができ汎用的である。また、リストデータについて並列処理を行い、その並列処理をステップＳ４（順投影）およびステップＳ５（逆投影）でそれぞれ行っているので、メモリの競合を防止することができ、高速化が実現できる。その結果、汎用性が高く、演算の高速化を図ることができる。

Description

核医学用データ処理方法および核医学診断装置
　この発明は、放射線を検出することにより得られた事象データから生成されたリストデータを順投影演算および逆投影演算する核医学用データ処理方法および核医学診断装置に関する。
　上述した核医学診断装置、すなわちＥＣＴ(Emission Computed Tomography)装置として、ＰＥＴ(Positron Emission Tomography)装置を例に採って説明する。ＰＥＴ装置は、陽電子（Positron）、すなわちポジトロンの消滅によって発生する複数本の光子を検出して複数個の検出器でγ線を同時に検出したときのみ（つまり同時計数したときのみ）被検体の断層画像を再構成するように構成されている。
　このＰＥＴ装置では、放射性薬剤を被検体に投与した後、対象組織における薬剤蓄積の過程を経時的に測定することで、様々な生体機能の定量測定が可能である。したがって、ＰＥＴ装置によって得られる画像は機能情報を有する。
　具体的には、被検体として人体を例に採って説明すると、被検体の体内にポジトロン（陽電子）放射性の同位元素（例えば１５Ｏ、１８Ｆ、１１Ｃなど）を注入し、これらから放出されるポジトロンが電子と結合する際に発生するγ線を検出する。このγ線の検出を、被検体の長手方向の軸である体軸周りを取り囲むようにしてリング状に配置された多数のγ線検出器からなる検出器列により行う。そして、コンピュータにより通常のＸ線ＣＴ(Computed Tomography)と同様の手法で計算を行って面内で特定し、被検体のイメージを作成する。
　ＭＬ－ＥＭ法(Maximum Likelihood Expectation Maximization)，ＯＳＥＭ法(Ordered Subset Expectation Maximization)，ＲＡＭＬＡ法(Row-Action Maximum Likelihood Algorithm)，ＤＲＡＭＡ法(Dynamic Row-Action Maximum Likelihood Algorithm)等の逐次近似法（例えば、非特許文献１参照）は、ＰＥＴ装置のような核医学診断装置において画像を再構成する際には不可欠な技術となっている。
　なお、γ線を検出する事象を「イベント」と言い、同時計数されたデータを「コインシデンスデータ」と言い、同時計数されなかったデータを「シングルイベントデータ」と言う。ＰＥＴ装置のガントリから転送されてくるイベント毎のデータを「リストデータ」と呼ぶが、このリストデータは、イベント毎のＸ座標，Ｙ座標，Ｚ座標などの位置情報や時間情報などがあって、これらの情報が時系列に並んでいる。
　これに対して、このリストデータに対して時間で積分するヒストグラム処理を行うことでそのリストデータをヒストグラム化して、ヒストグラムデータを取得する。このヒストグラムデータは、所定の時間で積分されたデータである。このヒストグラムデータとしては、縦軸を投影方向，横軸を画素としたサイノグラムがある。
　近年では、このリストデータを用いた逐次近似法も行われている（例えば、非特許文献２参照）。しかし、リストデータを用いたモード（「リストモード」と呼ぶ）の場合には、事象（イベント）の数だけデータが増大し、その演算量の多さが問題となる。このリストデータを用いた逐次近似法では演算量がさらに膨大となる。
　そこで、リストモードでは、並列演算の必要性が生じる。これまで多く行われていた並列演算はＭＰＩ(Message Passing Interface)などのＭＩＭＤ(Multiple Instruction Multiple Data)型演算である。ＭＩＭＤ型演算では、複数のプロセッサが複数の異なるデータを並行処理することから、メモリ領域が多く必要となる上、多数のプロセッサが必要となる。
　それに対して、近年では、ＧＰＵ(Graphics Processing Unit)を用いた逐次近似法の並列演算として、画像再構成の核とも言える順投影処理や逆投影処理の並列化の研究も行われている（例えば、非特許文献３参照）。ＧＰＵを用いた並列演算はＳＩＭＤ(Single Instruction Multiple Data)型演算に分類される。
　ＳＩＭＤ型演算では、１回の命令で複数データに対する処理を同時に行うことができるので、大量の演算に対して同一の演算処理を施す処理に向いている。しかし、ＧＰＵを用いて順投影処理や逆投影処理の並列演算を行う場合、同一メモリへの書き込みが発生（いわゆる「メモリの競合」である）する可能性があるので、並列化するためには工夫が必要である。特に、非特許文献３ではサイノグラムを用いており、サイノグラムでは、同一画素のデータが重複するので、メモリの競合が生じる。そこで、同じＧＰＵを用いた並列演算として、リストデータを用いた逐次近似法の並列の研究も行われている（例えば、非特許文献４参照）。
田中栄一，　「ＰＥＴ画像の再構成法の現状と展望」，日本放射線技術学会雑誌，　浜松ホトニクス株式会社，　Vol.62, No.6, 771-777,(2006). A. Fukano, F. Nakayama, H. Kubo, "Performance evaluation of relaxed block-iterative algorithms for 3-D PET reconstruction," IEEE Trans. Nucl. Sci., Vol.5, pp.2830-2834, 2004. H. Yang, M. Li, K. Koizumi, and H. Kubo, "Accelerating backprojections via CUDA architecture," in Proceedings of the 9th International Meeting on Fully Three-Dimensional Image Reconstruction in Radiology and Nuclear Medicine, pp.52-55, Lindau, Germany, July 2007. Guillem Pratx, Craig S. Levin et al, "Fast, Accurate and Shift-Varying Line Projections for Iterative Reconstruction Using the GPU", IEEE Trans. Med. Imaging., Vol.28, No 3, pp.435-445, 2009.
　しかしながら、上述の非特許文献４に報告されている手法は、ＧＰＵ独特の機構を利用するものであり、ＧＰＵ以外の他のＳＩＭＤ型機構への応用は困難であるという問題点がある。ＧＰＵはそもそも画像処理のためにあるものであったので、ＧＰＵ独特の機構の利用とは、データをテクスチャ画像として扱った上で、ピクセル（画素）シェーダやバーテックスシェーダを用いることを指す。
　この発明は、このような事情に鑑みてなされたものであって、汎用性が高く、演算の高速化を図ることができる核医学用データ処理方法および核医学診断装置を提供することを目的とする。
　この発明は、このような目的を達成するために、次のような構成をとる。
　すなわち、この発明の核医学用データ処理方法は、放射線を検出することにより得られた事象データから生成されたリストデータをＬＯＲ毎に並列に順投影演算して順投影処理後データを得る順投影処理工程と、前記順投影処理後データあるいはリストデータを画像空間の断面を分割した領域毎に並列に逆投影演算する逆投影処理工程とを備えることを特徴とするものである。
　この発明の核医学用データ処理方法によれば、順投影処理工程は、放射線を検出することにより得られた事象データから生成されたリストデータをＬＯＲ(Line Of Response)毎に並列に順投影演算して順投影処理後データを得る。逆投影処理工程は、順投影処理後データあるいはリストデータを画像空間の断面を分割した領域毎に並列に逆投影演算する。ＬＯＲ毎にリストデータについて並列に演算処理を行っているので、例えばＧＰＵ独特の機構を利用せずとも、順投影演算よび逆投影演算が実現可能で、他の演算機構に適用することができ汎用的である。また、リストデータについて並列処理を行い、その並列処理を順投影演算および逆投影演算でそれぞれ行っているので、メモリの競合を防止することができ、高速化が実現できる。また、逆投影演算では、順投影演算と相違して、メモリの競合が発生する可能性があるので、それを防止するために、逆投影処理工程は、画像空間の断面を分割した領域毎に並列に逆投影演算している。その結果、汎用性が高く、演算の高速化を図ることができる。
　各ＬＯＲの主方向はそれぞれランダムに並んでいるので、最初に主方向の分類を行うのが好ましい。すなわち、この発明の核医学用データ処理方法において、ＬＯＲのベクトルと放射線を検出する検出器の配置方向とに基づいてＬＯＲの主方向を主方向決定工程は決定するのが好ましい。そして、その主方向決定工程で決定された主方向に基づいてＬＯＲが通過する画像空間の断面を断面決定工程は決定する。その断面決定工程で決定された断面に基づいてその断面におけるＬＯＲが横断する演算領域を演算領域決定工程は演算領域を決定する。さらに具体的に順投影処理工程および逆投影処理工程は以下のような演算処理を行う。すなわち、順投影処理工程は、その演算領域毎に並列に演算し、逆投影処理工程は、演算領域内を複数の領域に分割した領域毎に順投影処理後データあるいはリストデータを並列に演算する。このように、主方向を決定することで主方向を分類して、それぞれの演算処理が行いやすくなる。
　また、この発明の核医学用データ処理方法では、上述の逆投影処理工程は、上述の断面を複数の領域に分割して、分割されて分散された複数の領域を１つの領域に所属させて、所属されたそれぞれの領域毎に並列に演算を行うのが好ましい。このように領域についても分類することで、特定のスレッドに局在することを防ぎ、演算効率を向上させることができる。ここで、「スレッド(Thread)」とは並列データの処理単位を示す。
　より好ましい一例は、所定間隔離れた複数の行領域を上述の１つの領域として所属させることである。このように領域をまんべんなく分散させて分類することで、特定のスレッドに局在することをより一層防止することができる。
　上述したこれらの発明の核医学用データ処理方法の一例は、順投影処理工程および逆投影処理工程を複数回実行して、画像を逐次に近似して更新する逐次近似法を用いることである。核医学診断の分野では、上述のように逐次近似法を用いて画像再構成が行われる。したがって、この発明においても逐次近似法に適用してもよい。もちろん、逆投影処理工程において、フィルタード・バックプロジェクション（FBP: Filtered Back Projection）（「フィルタ補正逆投影法」とも呼ばれる）を適用してもよい。
　また、この発明の核医学診断装置は、放射線を検出することにより得られた事象データから生成されたリストデータをＬＯＲ毎に並列に順投影演算して順投影処理後データを得る順投影処理手段と、前記順投影処理後データを画像空間の断面を分割した領域毎に並列に逆投影演算する逆投影処理手段とを備えることを特徴とするものである。
　この発明の核医学診断装置によれば、順投影処理手段は、放射線を検出することにより得られた事象データから生成されたリストデータをＬＯＲ毎に並列に順投影演算して順投影処理後データを得て、逆投影処理手段は、順投影処理後データあるいはリストデータを画像空間の断面を分割した領域毎に並列に逆投影演算するので、核医学用データ処理方法と同様に、汎用性が高く、演算の高速化を図ることができる。
　この発明に係る核医学用データ処理方法および核医学診断装置によれば、放射線を検出することにより得られた事象データから生成されたリストデータをＬＯＲ毎に並列に順投影演算して順投影処理後データを得て、順投影処理後データあるいはリストデータを画像空間の断面を分割した領域毎に並列に逆投影演算するので、汎用性が高く、演算の高速化を図ることができる。
実施例に係るＰＥＴ(Positron Emission Tomography)装置の側面図およびブロック図 γ線検出器の概略斜視図である。 一連のリストデータの処理方法の流れを示すフローチャートである。 検出確率の説明に供するγ線検出器での同時計数を示した模式図である。 リング状に配置されたγ線検出器の配置と３次元座標との関係を定義した模式図である。 ＬＯＲの方向ベクトルを投影の概念図に併記した模式図である。 ＬＯＲが通過する画像空間の断面における演算領域の模式図である。 ＬＯＲおよび１サブセット分の投影値の並列演算の説明に供する模式図である。 （ａ）はＬＯＲが通過する画像空間の断面における演算領域の分割を示した模式図、（ｂ）は（ａ）の拡大図である。 一連のステップＳ１～Ｓ５でのリストデータの流れを示す模式図である。
　３　…　γ線検出器
　１０　…　演算処理部
　ｄ　…　（方向）ベクトル
　Ｍ　…　被検体
　以下、図面を参照してこの発明の実施例を説明する。図１は、実施例に係るＰＥＴ(Positron Emission Tomography)装置の側面図およびブロック図であり、図２は、γ線検出器の概略斜視図である。
　本実施例に係るＰＥＴ装置は、図１に示すように、被検体Ｍを載置する天板１を備えている。この天板１は、上下に昇降移動、被検体Ｍの体軸Ｚに沿って平行移動するように構成されている。このように構成することで、天板１に載置された被検体Ｍは、後述するガントリ２の開口部２ａを通って、頭部から順に腹部、足部へと走査されて、被検体Ｍの画像を得る。なお、走査される部位や各部位の走査順序については特に限定されない。
　天板１の他に、本実施例に係るＰＥＴ装置は、開口部２ａを有したガントリ２と、γ線検出器３とを備えている。γ線検出器３は、被検体Ｍの体軸Ｚ周りを取り囲むようにしてリング状に配置されており、ガントリ２内に埋設されている。γ線検出器３は、この発明における検出器に相当する。
　その他にも、本実施例に係るＰＥＴ装置は、天板駆動部４とコントローラ５と入力部６と出力部７とメモリ部８と同時計数回路９と演算処理部１０とを備えている。天板駆動部６は、天板１の上述した移動を行うように駆動する機構であって、図示を省略するモータなどで構成されている。演算処理部１０は、この発明における順投影処理手段および逆投影処理手段に相当する。
　コントローラ５は、本実施例に係るＰＥＴ装置を構成する各部分を統括制御する。コントローラ５および演算処理部１０は、中央演算処理装置（ＣＰＵ）などで構成されている。特に、本実施例では、演算処理部１０は、ＧＰＵなどに代表されるＳＩＭＤ型機構で構成されている。もちろん、演算処理部１０は、ＧＰＵ以外の他のＳＩＭＤ型機構であってもよく、並列演算が可能であれば、特に限定されない。
　入力部６は、オペレータが入力したデータや命令をコントローラ５に送り込む。入力部６は、マウスやキーボードやジョイスティックやトラックボールやタッチパネルなどに代表されるポインティングデバイスで構成されている。出力部７はモニタなどに代表される表示部やプリンタなどで構成されている。
　メモリ部８は、ＲＯＭ（Read-only Memory）やＲＡＭ（Random-Access Memory）などに代表される記憶媒体で構成されている。本実施例では、同時計数回路９で同時計数された計数値（カウント）や同時計数した２つのγ線検出器３からなる検出器対やＬＯＲといった同時計数に関するデータや、演算処理部１０で演算処理された各種のデータなどについてはＲＡＭに書き込んで記憶し、必要に応じてＲＡＭから読み出す。ＲＯＭには、各種の核医学診断を含めて撮像を行うためのプログラム等を予め記憶しており、そのプログラムをコントローラ５および演算処理部１０が実行することでそのプログラムに応じた核医学診断をそれぞれ行う。特に、本実施例では、ＧＰＵを用いた並列演算を「ＣＵＤＡ（ＮＶｉｄｉａ社提供）」と呼ばれる並列コンピューティングアーキテクチャに実行させるべく、ＣＵＤＡに関するプログラムをＲＯＭに予め記憶しており、そのＣＵＤＡに関するプログラムを演算処理部１０が実行することで後述するステップＳ１～Ｓ５の処理を実行する。
　放射性薬剤が投与された被検体Ｍから発生したγ線をγ線検出器３のシンチレータブロック３１（図２を参照）が光に変換して、変換されたその光をγ線検出器３の光電子増倍管(PMT: Photo Multiplier Tube)３３（図２を参照）は増倍させて電気信号に変換する。その電気信号をイベントとして同時計数回路９に送り込む。
　具体的には、被検体Ｍに放射性薬剤を投与すると、ポジトロン放出型のＲＩのポジトロンが消滅することにより、２本のγ線が発生する。同時計数回路９は、シンチレータブロック３１（図２を参照）の位置とγ線の入射タイミングとをチェックし、被検体Ｍの両側にある２つのシンチレータブロック３１でγ線が同時に入射したときのみ、送り込まれたイベントを適正なデータと判定する。一方のシンチレータブロック３１のみにγ線が入射したときには、同時計数回路９は棄却する。つまり、同時計数回路９は、上述した電気信号に基づいて、２つのγ線検出器３においてγ線が同時観測されたことを検出する。
　同時計数回路９に送り込まれたイベントを、演算処理部１０に送り込む。演算処理部１０は順投影処理や逆投影処理による画像再構成を行って、被検体Ｍの画像を求める。画像を、コントローラ５を介して出力部７に送り込む。このようにして、演算処理部１０で得られた画像に基づいて核医学診断を行う。演算処理部１０の具体的な機能については後述する。
　γ線検出器３は、図２に示すようにシンチレータブロック３１と、そのシンチレータブロック３１に対して光学的に結合されたライトガイド３２と、そのライトガイド３２に対して光学的に結合された光電子増倍管（以下、単に「ＰＭＴ」と略記する）３３とを備えている。シンチレータブロック３１を構成する各シンチレータ素子は、γ線の入射に伴って発光することでγ線から光に変換する。この変換によってシンチレータ素子はγ線を検出する。シンチレータ素子において発光した光がシンチレータブロック３１で十分に拡散されて、ライトガイド３２を介してＰＭＴ３３に入力される。ＰＭＴ３３は、シンチレータブロック３１で変換された光を増倍させて電気信号に変換する。その電気信号は、上述したようにイベントとして同時計数回路９（図１を参照）に送り込まれる。
　次に、演算処理部１０の具体的な機能について、図３～図１０を参照して説明する。図３は、一連のリストデータの処理方法の流れを示すフローチャートであり、図４は、検出確率の説明に供するγ線検出器での同時計数を示した模式図であり、図５は、リング状に配置されたγ線検出器の配置と３次元座標との関係を定義した模式図であり、図６は、ＬＯＲの方向ベクトルを投影の概念図に併記した模式図であり、図７は、ＬＯＲが通過する画像空間の断面における演算領域の模式図であり、図８は、ＬＯＲおよび１サブセット分の投影値の並列演算の説明に供する模式図であり、図９は、ＬＯＲが通過する画像空間の断面における演算領域の分割を示した模式図であり、図１０は、一連のステップＳ１～Ｓ５でのリストデータの流れを示す模式図である。図４、図５ではγ線検出器３として、シンチレータブロック３１のみを図示して、ライトガイド３２やＰＭＴ３３については図示を省略する。また、本実施例では、３次元のボクセル(voxel)で構成される画素を例に採って説明する。
　本実施例では、逐次近似法アルゴリズムとしてリストデータを用いたlist-mode　３Ｄ－ＤＲＡＭＡ法(Dynamic Row-Action Maximum Likelihood Algorithm)を採用し、吸収補正、感度補正、ランダム補正等の補正を取り入れた上で、ＧＰＵを用いた並列演算を上述のＣＵＤＡで実装している。list-mode　３Ｄ－ＤＲＡＭＡ法の更新式は、下記（１）式で表される。
Figure JPOXMLDOC01-appb-M000001
　なお、上記（１）式中のｉ（ｔ）はｔ番目の同時計数イベントを示す。総イベント数をサブセット数Ｌで除算し、それぞれのサブセットをＳｌ（ｌ＝０，１，…，Ｌ－１）（ｌは小文字のＬ）と表す。ｋ（ｋ＝０，１，…）は逐次近似回数、ｌ（小文字のＬ）はサブセット数であり、ｘｊ （ｋ，ｌ）は更新後のボクセルｊ（ｊ＝０，１，…，Ｊ－１）の画素値である。図４に示すように、ａｉ（ｔ）ｊはボクセルｊから放出されたγ線フォトンがＬＯＲｉ（ｔ）に検出される確率であり、「システムマトリックス」とも呼ばれている。システムマトリックスの演算では検出器応答関数も考慮している。Ｃｊは正規化マトリックス、ｐｌｊはブロックングファクタを示す。
　ブロッキングファクタについては、上述の非特許文献２で２種紹介されているが、本検討では収束性が高く、すべてのリストデータを用いて演算する手法を採用している。λ（ｋ，ｌ）は緩和係数であり、収束速度に関与する。γとβ０はλ（ｋ，ｌ）を制御するパラメータである。また、ａ´ｉ（ｔ）ｊは吸収補正後のａｉ（ｔ）ｊであり、下記（２）式に補正関連の定義を示す。
Figure JPOXMLDOC01-appb-M000002
　なお、上記（２）式中のＡｉ（ｔ）ｊは吸収補正係数を表し、μｊは吸収係数の分布を表す。また、ｒｉはシングル計数の計数率ｓｉ０，ｓｉ１と同時計数のタイムウィンドウ２τよりランダム計数率を予測した値である。
　また、上記（１）式中のＦＰはForward Projection、すなわち順投影を表し、上記（１）式中の部分で順投影処理での演算が行われる。一方、上記（１）式中のＢＰはBack Projection、すなわち逆投影を表し、上記（１）式中の部分で逆投影処理での演算が行われる。
　先ず、初期画像であるｘｊ （０，０）を適宜に設定する。初期画像ｘｊ （０，０）については、例えば一様な画素値を有する画像であればよく、ｘｊ （０，０）＞０とする。設定された初期画像ｘｊ （０，０）と、上記（２）式で補正されたａ´ｉ（ｔ）ｊとを用いて、上記（１）式に繰り返し代入することで、ｘｊ （０，０），…，ｘｊ （０，Ｌ－１）が逐次に求められ、最終的に求められたｘｊ （０，Ｌ－１）をｘｊ （１，０）とすることでｘｊ （１，０）に繰り上げる。以下、同様に、ｘｊを順に繰り上げる（ｘｊ （０，０），ｘｊ （１，０）…，ｘｊ （ｋ，０））。反復を表すｋの回数については特に限定されず、適宜に設定すればよい。このように最終的に求められたｘｊをそれに対応するボクセルνｊごとに並べることで、順投影処理や逆投影処理による画像再構成を行い、被検体Ｍの画像を求める。
　続いて、各々の順投影処理および逆投影処理での具体的な並行演算について詳しく述べる。
　（ステップＳ１）主方向決定
　図５に示すように、水平軸をｘ、垂直軸をｙとし、被検体Ｍの体軸をｚとする。リング状に配置されたγ線検出器３をｘｙ平面上に並べたときには、体軸ｚは奥行き方向の軸になる。図６に示すように、ＬＯＲに平行な方向ベクトルをｄ（ｄｘ，ｄｙ，ｄｚ）と定義し（ただし||ｄ||＝１、単位ベクトル）、ＬＯＲのベクトルｄとγ線検出器３の配置方向とに基づいて、ＬＯＲの主方向を決定する。γ線検出器３はｘｙ平面上で配置されているので、ベクトルｄ（ｄｘ，ｄｙ，ｄｚ）のうちｄｘとｄｙとの大きさを比較して大きい方を主方向とする。図６では、ｘ方向を主方向としているときの模式図である。このステップＳ１は、この発明における主方向決定工程に相当する。
　（ステップＳ２）断面決定
　ステップＳ１で決定された主方向に基づいてＬＯＲが通過する画像空間の断面を決定する。本実施例では、主方向が直交する面を、ＬＯＲが通過する画像空間の断面として決定する。図６では、ｘ方向を主方向としている図であるので、図６、図７に示すように断面はｙｚ平面となる。このステップＳ２は、この発明における断面決定工程に相当する。
　（ステップＳ３）演算領域決定
　ステップＳ２で決定された断面（図６ではｙｚ平面）に基づいてその断面におけるＬＯＲが横断する演算領域を決定する。まず、ＬＯＲの始点となる断面（図６では「slice_start」で表記）およびＬＯＲの終点となる断面（図６では「slice_end」で表記）を決定して、始点となる断面上の座標と終点となる断面上の座標とを決定する。図６、図７では、断面をｙｚ平面としているので、始点となる断面上の座標はｙ，ｚ座標となり、同じく終点となる断面上の座標をｙ，ｚ座標となる。ベクトルは一方向に延びているので、始点・終点のｙ座標のいずれか一方が最小値の座標（図７では「Min_0」で表記）となり、他方が最大値の座標（図７では「Max_0」で表記）となる。同じく、始点・終点のｚ座標のいずれか一方が最小値の座標（図７では「Min_1」で表記）となり、他方が最大値の座標（図７では「Max_1」で表記）となる。つまり、断面上では、ＬＯＲは各断面上の最小値の座標位置から最大値の座標位置を横切ることになる。したがって、ＬＯＲが横断する演算領域は、図７に示すようにMin_0，Max_0，Min_1およびMax_1で区切られた領域となる。このステップＳ３は、この発明における演算領域決定工程に相当する。
　（ステップＳ４）順投影処理
　ステップＳ３で決定された演算領域でリストデータについて並列に演算処理を行う。具体的には、それぞれのボクセルｊについてシステムマトリックスを求めて、決定された演算領域の範囲で上記（１）式中のＦＰで表記された部分において画素値と乗算して加算することにより投影値の更新を行う。上記（１）式での並列演算をＬＯＲ毎に行う。投影値を代入するメモリについては、１サブセット当たり、１イベント当たり１要素を書き込む構成となり、例えば図８のように書き込まれる。図８中の「Thread_0」、「Thread_1」、「Thread_2」、「Thread_3」、…は、書き込まれる並列データの処理単位を示す。各ＬＯＲは、１サブセット当たり各スレッドに一対一に対応している。
　例えば、総イベント数が3276800イベントの処理を行う場合に、サブセット数を128とすると、１サブセット当たりイベント数が25600（＝3276800/128）となる。したがって、順投影処理での並列演算を行う場合には、25600イベントを25600スレッドで演算することから、投影値の書き込み先はスレッド毎に異なり、メモリの競合は発生しない。このステップＳ４は、この発明における順投影処理工程に相当する。
　（ステップＳ５）逆投影処理（領域分割）
　ステップＳ４で順投影処理が行われた順投影処理後データあるいはリストデータについて逆投影処理を行う。この場合には、ＬＯＲが通過する画像空間の各断面を複数の領域に分割して、分割されて分散された複数の領域を１つの領域に所属させて、所属されたそれぞれの領域毎に並列に演算を行う。具体的には、所属されたそれぞれの領域で上記（１）式中のＢＰで表記された部分において加算することにより逆投影値の更新を行う。
　図９を参照してより詳しく説明すると、所定間隔離れた行領域を１つの領域として所属させる。例えば、演算領域が256×256ピクセルの断面（スライス）の場合に、３１行おきの行を１つの領域と定義すると、１スライス当たり３２の領域に分割したことになり、８（＝256/32）行で１つの領域となる。したがって、この場合には、３１行間隔離れた行領域毎に１つの領域として所属させることになる。図９の場合には、図示の便宜上、７行おきの行の１つの領域として定義している。図９（ａ）を拡大したのが図９（ｂ）である。図９に示すように、７行おきの行の１つの領域として定義すると、１スライス当たり８の領域に分割することになる。すなわち、図９（ｂ）で表記された「Thread_0」、「Thread_1」、…「Thread_7」でそれぞれ分割することになる。このステップＳ５は、この発明における逆投影処理工程に相当する。
　また、一連のステップＳ１～Ｓ５でのリストデータの流れは、図１０に示す通りである。すなわち、各イベント毎（すなわち各ＬＯＲ毎）のリストデータ（図１０では「Event_0」、「Event_1」、「Event_2」、…で表記）を、ステップＳ１で主方向を決定することで主方向を分類する。本実施例では、ｄｘとｄｙとの大きさを比較して大きい方を主方向としているので、主方向がｘ方向のリストデータ（図１０では「主方向ｘのリストデータ」で表記）および主方向がｙ方向のリストデータ（図１０では「主方向ｙのリストデータ」で表記）に分類して、それぞれステップＳ５で逆投影処理について並列に演算処理を行う。
　なお、上述した逐次近似法としては、上述したＤＲＡＭＡ法に限定されず、スタティックな（つまり静的な）ＲＡＭＬＡ法(Row-Action Maximum Likelihood Algorithm)でもよいし、ＭＬ－ＥＭ法(Maximum Likelihood Expectation Maximization)でもよいし、ＯＳＥＭ法(Ordered Subset ML-EM)でもよい。
　上述の構成を備えた本実施例に係るＰＥＴ装置および処理方法によれば、ステップＳ４（順投影処理）は、γ線を検出することにより得られた事象データから生成されたリストデータをＬＯＲ毎に並列に順投影演算して順投影処理後データを得る。ステップＳ５（逆投影処理）は、順投影処理後データあるいはリストデータを画像空間の断面を分割した領域毎に並列に逆投影演算する。ＬＯＲ毎にリストデータについて並列に演算処理を行っているので、例えばＧＰＵ独特の機構を利用せずとも、順投影演算よび逆投影演算が実現可能で、他の演算機構に適用することができ汎用的である。また、リストデータについて並列処理を行い、その並列処理を順投影演算および逆投影演算でそれぞれ行っているので、メモリの競合を防止することができ、高速化が実現できる。また、逆投影演算では、順投影演算と相違して、メモリの競合が発生する可能性があるので、それを防止するために、ステップＳ５（逆投影処理）は、画像空間の断面を分割した領域毎に並列に逆投影演算している。その結果、汎用性が高く、演算の高速化を図ることができる。
　各ＬＯＲの主方向はそれぞれランダムに並んでいるので、本実施例では、好ましくは最初に主方向の分類を行っている。すなわち、ＬＯＲのベクトルｄとγ線を検出するγ線検出器３の配置方向とに基づいてＬＯＲの主方向をステップＳ１（主方向決定）は好ましくは決定している。そして、そのステップＳ１（主方向決定）で決定された主方向に基づいてＬＯＲが通過する画像空間の断面をステップＳ２（断面決定）は決定する。そのステップＳ２（断面決定）で決定された断面に基づいてその断面におけるＬＯＲが横断する演算領域をステップＳ３（演算領域決定）は演算領域を決定する。さらに具体的にステップＳ４（順投影処理）およびステップＳ５（逆投影処理）は上述したような演算処理を行っている。すなわち、ステップＳ４（順投影処理）は、その演算領域毎に並列に演算し、ステップＳ５（逆投影処理）は、演算領域内を複数の領域に分割した領域毎に順投影処理後データあるいはリストデータを並列に演算する。このように、主方向を決定することで主方向を分類して、それぞれの演算処理が行いやすくなる。
　また、本実施例では、上述のステップＳ５（逆投影処理）は、上述の断面を複数の領域に分割して、分割されて分散された複数の領域を１つの領域に所属させて、所属されたそれぞれの領域毎に並列に演算を好ましくは行っている。このように領域についても分類することで、特定のスレッドに局在することを防ぎ、演算効率を向上させることができる。
　より好ましくは、所定間隔離れた行領域毎に１つの領域として所属させている（図９を参照）。このように領域をまんべんなく分散させて分類することで、特定のスレッドに局在することをより一層防止することができる。
　本実施例では、上述のステップＳ４（順投影処理）およびステップＳ５（逆投影処理）を複数回実行して、画像を逐次に近似して更新する逐次近似法に適用している。すなわち、逐次近似法を用いて並列に上述の演算処理を行っている。核医学診断の分野では、上述のように逐次近似法を用いて画像再構成が行われる。
　また、本実施例に係るＰＥＴ装置は、γを検出することにより得られた事象データから生成されたリストデータをＬＯＲ毎に並列に順投影演算して順投影処理後データを得て、順投影処理後データを画像空間の断面を分割した領域毎に並列に逆投影演算する上述のステップＳ１～Ｓ５の演算処理を行う演算処理部１０とを備えている。本実施例のステップＳ１～Ｓ５の方法をＰＥＴ装置に適用することで、本実施例のステップＳ１～Ｓ５の方法と同様に、汎用性が高く、演算の高速化を図ることができる。
　本実施例では、上述のlist-mode　３Ｄ－ＤＲＡＭＡ法のＣＵＤＡ実装を行い、画像再構成の高速化の検討を行った結果、３．２Ｍイベントのリストデータについて約１８倍の高速化が実現したのを確認した。また、ＧＰＵを用いて演算して得られた再構成画像は一般のＣＰＵを用いた場合とほぼ同一であることも確認した。
　この発明は、上記実施形態に限られることはなく、下記のように変形実施することができる。
　（１）上述した実施例では、逐次近似法を用いて上述のステップＳ４（順投影処理）およびステップＳ５（逆投影処理）で並列に上述の演算処理を行ったが、逆投影処理工程において、フィルタード・バックプロジェクションを用いたフェルドカンプ(Feldkamp)法を適用してもよい。
　（２）上述した実施例では、主方向を決定することで主方向を分類したが、必ずしも主方向を決定する必要はない。ただし、それぞれの演算処理が行いやすくなるためには、実施例のように、主方向を決定することで主方向を分類するのがより好ましい。
　（３）上述した実施例では、水平軸をｘ、垂直軸をｙとし、被検体Ｍの体軸をｚとした上で主方向を決定したが、必ずしもこれらの軸の方向に限定されない。
　（４）上述した実施例では、リング状に配置されたγ線検出器３を図５に示すようにｘｙ平面に並べた上で主方向を決定したが、γ線検出器の配置状態に応じて主方向を適宜決定すればよい。非リング状の検出器の場合には、同時計数する２つの検出器上の平面に基づいて主方向を決定すればよい。
　（５）上述した実施例では、主方向を決定する際に大きい方を主方向としたが、必ずしもこれに限定されない。一定の規則に基づいて主方向を決定するのであれば、例えば小さい方を主方向としてもよい。
　（６）上述した実施例では、主方向が直交する面を、ＬＯＲが通過する画像空間の断面として決定したが、必ずしも主方向が直交する面を当該断面として決定する必要はない。９０°以外の角度で横切る面を当該断面として決定してもよい。
　（７）上述した実施例では、所定間隔離れた行領域毎に１つの領域として所属させたが、行領域以外のブロックの領域で区分させて、それらのいくつかをまとめて１つの領域として所属させてもよい。
　（９）上述した実施例では、核医学診断装置について述べたが、コリメータ機構のあるガンマカメラ、ＳＰＥＣＴにも応用が可能である。

Claims (6)

  1. 　放射線を検出することにより得られた事象データから生成されたリストデータをＬＯＲ毎に並列に順投影演算して順投影処理後データを得る順投影処理工程と、
    　前記順投影処理後データあるいはリストデータを画像空間の断面を分割した領域毎に並列に逆投影演算する逆投影処理工程と
    　を備えることを特徴とする核医学用データ処理方法。
  2. 　請求項１に記載の核医学用データ処理方法において、
    　前記ＬＯＲのベクトルと前記放射線を検出する検出器の配置方向とに基づいて前記ＬＯＲの主方向を決定する主方向決定工程と、
    　その主方向決定工程で決定された前記主方向に基づいて前記ＬＯＲが通過する画像空間の断面を決定する断面決定工程と、その断面決定工程で決定された前記断面に基づいてその断面における前記ＬＯＲが横断する演算領域を決定する演算領域決定工程とを備え、
    　前記順投影処理工程は、前記演算領域毎に並列に演算し、
    　前記逆投影処理工程は、前記演算領域内を複数の領域に分割した領域毎に並列に演算することを特徴とする核医学用データ処理方法。
  3. 　請求項１または請求項２に記載の核医学用データ処理方法において、
    　前記逆投影処理工程は、前記断面を複数の領域に分割して、分割されて分散された複数の領域を１つの領域に所属させて、所属されたそれぞれの領域毎に並列に演算を行うことを特徴とする核医学用データ処理方法。
  4. 　請求項３に記載の核医学用データ処理方法において、
    　所定間隔離れた複数の行領域を前記１つの領域として所属させることを特徴とする核医学用データ処理方法。
  5. 　請求項１から請求項４のいずれかに記載の核医学用データ処理方法において、前記順投影処理工程および前記逆投影処理工程を複数回実行して、画像を逐次に近似して更新する逐次近似法を用いることを特徴とする核医学用データ処理方法。
  6. 　放射線を検出することにより得られた事象データから生成されたリストデータをＬＯＲ毎に並列に順投影演算して順投影処理後データを得る順投影処理手段と、
    　前記順投影処理後データを画像空間の断面を分割した領域毎に並列に逆投影演算する逆投影処理手段と
    　を備えることを特徴とする核医学診断装置。
PCT/JP2009/004381 2009-09-04 2009-09-04 核医学用データ処理方法および核医学診断装置 WO2011027402A1 (ja)

Priority Applications (4)

Application Number Priority Date Filing Date Title
JP2011529700A JP5263402B2 (ja) 2009-09-04 2009-09-04 核医学用データ処理方法および核医学診断装置
CN200980161257.2A CN102483459B (zh) 2009-09-04 2009-09-04 核医学用数据处理方法以及核医学诊断装置
US13/388,792 US8987674B2 (en) 2009-09-04 2009-09-04 Data processing method for nuclear medicine, and a nuclear medicine diagnostic apparatus
PCT/JP2009/004381 WO2011027402A1 (ja) 2009-09-04 2009-09-04 核医学用データ処理方法および核医学診断装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/JP2009/004381 WO2011027402A1 (ja) 2009-09-04 2009-09-04 核医学用データ処理方法および核医学診断装置

Publications (1)

Publication Number Publication Date
WO2011027402A1 true WO2011027402A1 (ja) 2011-03-10

Family

ID=43648964

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2009/004381 WO2011027402A1 (ja) 2009-09-04 2009-09-04 核医学用データ処理方法および核医学診断装置

Country Status (4)

Country Link
US (1) US8987674B2 (ja)
JP (1) JP5263402B2 (ja)
CN (1) CN102483459B (ja)
WO (1) WO2011027402A1 (ja)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR101356881B1 (ko) 2012-04-02 2014-02-12 한국과학기술원 고해상도 양전자 방출 단층 촬영에서 병렬 처리를 위해 영상을 재구성하는 방법 및 장치
WO2015056299A1 (ja) * 2013-10-15 2015-04-23 株式会社島津製作所 断層画像処理方法およびそれを用いた放射型断層撮影装置
JP2019120530A (ja) * 2017-12-28 2019-07-22 浜松ホトニクス株式会社 画像処理装置および画像処理方法
WO2021186504A1 (ja) * 2020-03-16 2021-09-23 株式会社島津製作所 リストモード画像再構成方法および核医学診断装置

Families Citing this family (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9111381B2 (en) * 2010-01-27 2015-08-18 Koninklijke Philips N.V. Shift-varying line projection using graphics hardware
NL2010267C2 (en) * 2013-02-07 2014-08-11 Milabs B V High energy radiation detecting apparatus and method.
US9839401B2 (en) * 2016-01-19 2017-12-12 Shimadzu Corporation Radiation tomography apparatus
EP3350776B1 (en) 2016-04-20 2021-09-22 Shanghai United Imaging Healthcare Co., Ltd. System and method for image reconstruction
WO2017214766A1 (zh) * 2016-06-12 2017-12-21 上海联影医疗科技有限公司 正电子发射断层成像系统及其图像重建方法
JP6932253B2 (ja) * 2017-09-21 2021-09-08 コーニンクレッカ フィリップス エヌ ヴェＫｏｎｉｎｋｌｉｊｋｅ Ｐｈｉｌｉｐｓ Ｎ．Ｖ． 陽電子放出断層撮影の偶発同時推定のための緩和された反復的最大尤度期待値最大化
CN107908363B (zh) * 2017-11-10 2021-12-07 湖北锐世数字医学影像科技有限公司 基于cuda的pet符合事件筛选的方法、系统及装置
CN107945203A (zh) * 2017-11-24 2018-04-20 中国科学院高能物理研究所 Pet图像处理方法及装置、电子设备、存储介质
US11232612B2 (en) * 2019-03-15 2022-01-25 University Of Florida Research Foundation, Incorporated Highly accurate and efficient forward and back projection methods for computed tomography
CN111881412A (zh) * 2020-07-28 2020-11-03 南京航空航天大学 一种基于cuda的pet系统矩阵计算方法

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2007232548A (ja) * 2006-02-28 2007-09-13 Hitachi Ltd 核医学診断装置

Family Cites Families (15)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5224037A (en) * 1991-03-15 1993-06-29 Cti, Inc. Design of super-fast three-dimensional projection system for Positron Emission Tomography
JP4764050B2 (ja) * 2005-03-31 2011-08-31 株式会社日立製作所 核医学診断装置および核医学診断装置の冷却方法
US7498581B2 (en) * 2005-10-05 2009-03-03 Koninklijke Philips Electronics N.V. Distributed iterative image reconstruction
JP4679348B2 (ja) * 2005-11-22 2011-04-27 ジーイー・メディカル・システムズ・グローバル・テクノロジー・カンパニー・エルエルシー Ｘ線ｃｔ装置
WO2007100954A2 (en) * 2006-02-22 2007-09-07 Koninklijke Philips Electronics, N.V. Image reconstruction using data ordering
US8314796B2 (en) * 2006-02-24 2012-11-20 The Board Of Trustees Of The Leland Stanford Junior University Method of reconstructing a tomographic image using a graphics processing unit
JP4656233B2 (ja) * 2006-03-10 2011-03-23 株式会社島津製作所 核医学診断装置およびそれに用いられる診断システム
JP2007271509A (ja) * 2006-03-31 2007-10-18 Shimadzu Corp 核医学診断装置およびそれに用いられる診断システム
DE102007020879A1 (de) * 2006-05-10 2009-04-02 Gachon University Of Medicine & Science Industry-Academic Cooperation Foundation Verfahren und Vorrichtung für die äußerst schnelle Symmetrie- und SIMD- gestützte Projektion/Rückprojektion für die 3D-PET-Bildrekonstruktion
US7876941B2 (en) * 2007-03-09 2011-01-25 Siemens Medical Solutions Usa, Inc. Incorporation of axial system response in iterative reconstruction from axially compressed data of cylindrical scanner using on-the-fly computing
US7983465B2 (en) * 2007-05-09 2011-07-19 Société De Commercialisation Des Produits De La Recherche Appliquée - Socpra Sciences Santé Et Humaines, S.E.C. Image reconstruction methods based on block circulant system matrices
EP2156408B1 (en) * 2007-05-30 2021-03-17 Koninklijke Philips N.V. Pet local tomography
US8359345B2 (en) * 2008-05-09 2013-01-22 Siemens Medical Solutions Usa, Inc. Iterative algorithms for variance reduction on compressed sinogram random coincidences in PET
US8000513B2 (en) * 2008-09-22 2011-08-16 Siemens Medical Solutions Usa, Inc. System and method for 3D time of flight PET forward projection based on an exact axial inverse rebinning relation in fourier space
US8265365B2 (en) * 2010-09-20 2012-09-11 Siemens Medical Solutions Usa, Inc. Time of flight scatter distribution estimation in positron emission tomography

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2007232548A (ja) * 2006-02-28 2007-09-13 Hitachi Ltd 核医学診断装置

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR101356881B1 (ko) 2012-04-02 2014-02-12 한국과학기술원 고해상도 양전자 방출 단층 촬영에서 병렬 처리를 위해 영상을 재구성하는 방법 및 장치
WO2015056299A1 (ja) * 2013-10-15 2015-04-23 株式会社島津製作所 断層画像処理方法およびそれを用いた放射型断層撮影装置
JPWO2015056299A1 (ja) * 2013-10-15 2017-03-09 株式会社島津製作所 断層画像処理方法およびそれを用いた放射型断層撮影装置
JP2019120530A (ja) * 2017-12-28 2019-07-22 浜松ホトニクス株式会社 画像処理装置および画像処理方法
WO2021186504A1 (ja) * 2020-03-16 2021-09-23 株式会社島津製作所 リストモード画像再構成方法および核医学診断装置
JPWO2021186504A1 (ja) * 2020-03-16 2021-09-23
JP7302737B2 (ja) 2020-03-16 2023-07-04 株式会社島津製作所 リストモード画像再構成方法および核医学診断装置

Also Published As

Publication number Publication date
US20120126125A1 (en) 2012-05-24
US8987674B2 (en) 2015-03-24
CN102483459B (zh) 2016-08-24
JPWO2011027402A1 (ja) 2013-01-31
CN102483459A (zh) 2012-05-30
JP5263402B2 (ja) 2013-08-14

Similar Documents

Publication Publication Date Title
JP5263402B2 (ja) 核医学用データ処理方法および核医学診断装置
JP5152202B2 (ja) ポジトロンｃｔ装置
US8208599B2 (en) Iterative reconstruction with enhanced noise control filtering
US9155514B2 (en) Reconstruction with partially known attenuation information in time of flight positron emission tomography
US20100074500A1 (en) System and method for 3d time of flight pet forward projection based on an exact axial inverse rebinning relation in fourier space
Yamaya et al. First human brain imaging by the jPET-D4 prototype with a pre-computed system matrix
US8975587B2 (en) Positron CT apparatus and a reconstruction method
US9480450B2 (en) Scatter component estimating method
WO1998047103A1 (en) Direct tomographic reconstruction
US9916670B1 (en) Fast, efficient, and list-mode compatible tomographic image reconstruction using a novel quadratic surrogate
JP6014738B2 (ja) 三次元画像の投影方法
Szlávecz et al. GPU-based acceleration of the MLEM algorithm for SPECT parallel imaging with attenuation correction and compensation for detector response
KR101493683B1 (ko) 콘-빔 기반 반응선 재구성을 이용한 초고해상도 ｐｅｔ 영상 재구성 장치 및 그 방법
JP2010266235A (ja) ポジトロンｃｔ装置および再構成方法
US10354417B2 (en) Medical image processing apparatus and medical image diagnosis apparatus and medical image processing method
JP2010203806A (ja) 画像処理装置、画像再構成システム、画像処理方法およびプログラム
US20230206516A1 (en) Scatter estimation for pet from image-based convolutional neural network
WO2023228910A1 (ja) 画像処理装置および画像処理方法
JP6206501B2 (ja) 断層画像処理方法およびそれにおける各工程を実行するための手段を備えた放射線放出型断層撮影装置
Al-Saffar Decomposing of projected Footprints of Voxels in Separable Footprints 3D CT Reconstruction Method
Zou et al. Image Reconstruction for Source Trajectories with Kinks

Legal Events

Date Code Title Description
WWE Wipo information: entry into national phase

Ref document number: 200980161257.2

Country of ref document: CN

121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 09848933

Country of ref document: EP

Kind code of ref document: A1

ENP Entry into the national phase

Ref document number: 2011529700

Country of ref document: JP

Kind code of ref document: A

WWE Wipo information: entry into national phase

Ref document number: 13388792

Country of ref document: US

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 09848933

Country of ref document: EP

Kind code of ref document: A1