JP5342228B2

JP5342228B2 - Image processing apparatus and three-dimensional PET apparatus

Info

Publication number: JP5342228B2
Application number: JP2008330727A
Authority: JP
Inventors: 栄一田中; 卓志磯部; 希望大手
Original assignee: Hamamatsu Photonics KK
Current assignee: Hamamatsu Photonics KK
Priority date: 2008-12-25
Filing date: 2008-12-25
Publication date: 2013-11-13
Anticipated expiration: 2028-12-25
Also published as: JP2010151653A

Abstract

<P>PROBLEM TO BE SOLVED: To provide an image processing unit and a three-dimensional PET device for reducing required time for reconfiguring images by a block repetition method. <P>SOLUTION: The image processing unit uses projection data obtained by the three-dimensional PET device, where a plurality of detector rings including each of a plurality of radiation detectors are laminated axially, to reconfigure images by the block repetition method. The image processing unit creates a three-dimensional initial image by reconfiguring images using partial projection data including ring differences 0 and 1 in projection data (S10), reconfigures three-dimensional images mutually in parallel by the block repetition method using one portion of the projection data (S20<SB>0</SB>-S20<SB>N-1</SB>) starting from the initial images, and integrates the obtained three-dimensional images (S30). <P>COPYRIGHT: (C)2010,JPO&INPIT

Description

本発明は、３次元ＰＥＴ装置において得られた投影データを用いて画像再構成を行う画像処理装置に関するものである。 The present invention relates to an image processing apparatus that performs image reconstruction using projection data obtained in a three-dimensional PET apparatus.

ＰＥＴ装置は、陽電子放出アイソトープ（ＲＩ線源）が投入された生体（被検体）内における電子・陽電子の対消滅に伴って発生し互いに逆方向に飛行する光子対を同時計数法により検出することにより、その被検体内の極微量物質の挙動を画像化することができる装置である。ＰＥＴ装置は、被検体が置かれる測定空間の周囲に配列された多数の小型の放射線検出器を有する検出部を備えており、電子・陽電子の対消滅に伴って発生するエネルギ５１１ｋｅＶの光子対を同時計数法で検出し、この同時計数情報を蓄積してヒストグラムを作成する。そして、この作成されたヒストグラムに基づいて、測定空間における光子対の発生頻度の空間分布を表す画像を再構成する。このＰＥＴ装置は核医学分野等で重要な役割を果たしており、これを用いて例えば生体機能や脳の高次機能の研究を行うことができる。このようなＰＥＴ装置は、２次元ＰＥＴ装置と３次元ＰＥＴ装置とに大別される。 The PET device detects coincident photon pairs generated in the living body (subject) into which a positron emitting isotope (RI radiation source) has been introduced, which are generated when electrons and positrons are annihilated in opposite directions. Thus, the device can image the behavior of a trace amount substance in the subject. The PET apparatus includes a detection unit having a large number of small radiation detectors arranged around a measurement space in which a subject is placed, and generates a photon pair of energy 511 keV generated as a result of electron-positron pair annihilation. Detection is performed by the coincidence method, and the coincidence count information is accumulated to create a histogram. Based on the created histogram, an image representing the spatial distribution of the occurrence frequency of photon pairs in the measurement space is reconstructed. This PET apparatus plays an important role in the field of nuclear medicine and the like, and can be used to study, for example, biological functions and higher-order brain functions. Such PET apparatuses are roughly classified into two-dimensional PET apparatuses and three-dimensional PET apparatuses.

そのうち、２次元ＰＥＴ装置は、検出部が軸方向に積層された複数の検出器リングを含み、各検出器リングが複数の放射線検出器を含むものであって、検出器リング間にシールド板を有している。２次元ＰＥＴ装置の検出部は、中心軸との角度が約９０度の方向から飛来した光子対のみを同時計数することができる。すなわち、２次元ＰＥＴ装置の検出部により得られ蓄積された同時計数情報すなわち２次元投影データは、同一の検出器リングまたは隣接する（若しくは極めて近い）検出器リングに含まれる１対の放射線検出器によるものに限られる。したがって、２次元ＰＥＴ装置は、測定空間外の位置で発生した光子対が散乱された散乱線を効率よく除外することができ、また、２次元投影データに対して吸収補正や感度補正を容易に行うことができる。 Among them, the two-dimensional PET apparatus includes a plurality of detector rings in which the detectors are stacked in the axial direction, each detector ring includes a plurality of radiation detectors, and a shield plate is provided between the detector rings. Have. The detection unit of the two-dimensional PET apparatus can simultaneously count only photon pairs flying from a direction whose angle with the central axis is about 90 degrees. That is, the coincidence counting information obtained by the detection unit of the two-dimensional PET apparatus, that is, the two-dimensional projection data, is a pair of radiation detectors included in the same detector ring or adjacent (or very close) detector rings. Limited to Therefore, the two-dimensional PET apparatus can efficiently exclude scattered rays in which photon pairs generated at positions outside the measurement space are scattered, and absorption correction and sensitivity correction can be easily performed on two-dimensional projection data. It can be carried out.

一方、３次元ＰＥＴ装置は、検出部が軸方向に積層された複数の検出器リングを含み、各検出器リングが複数の放射線検出器を含むものであって、検出器リング間にシールド板を有していない。３次元ＰＥＴ装置の検出部は、あらゆる方向から飛来した光子対を同時計数することができる。すなわち、３次元ＰＥＴ装置の検出部により得られ蓄積される同時計数情報すなわち３次元投影データは、任意の検出器リングに含まれる１対の放射線検出器によるものが可能である。したがって、３次元ＰＥＴ装置は、２次元ＰＥＴ装置と比較して５倍〜１０倍程度の高い感度で、光子対を同時計数することができる。 On the other hand, the three-dimensional PET apparatus includes a plurality of detector rings in which the detectors are stacked in the axial direction, and each detector ring includes a plurality of radiation detectors, and a shield plate is provided between the detector rings. I don't have it. The detection unit of the three-dimensional PET apparatus can simultaneously count photon pairs flying from all directions. That is, the coincidence counting information obtained by the detection unit of the three-dimensional PET apparatus, that is, the three-dimensional projection data, can be obtained by a pair of radiation detectors included in an arbitrary detector ring. Therefore, the three-dimensional PET apparatus can simultaneously count photon pairs with a sensitivity about 5 to 10 times higher than that of the two-dimensional PET apparatus.

それ故、近年では、３次元ＰＥＴ装置の普及が急速に進みつつあり、また、３次元ＰＥＴ装置で得られた投影データに基づく画像再構成技術についての研究も為されている。３次元ＰＥＴ装置における画像再構成技術としては、ＭＬ-ＥＭ（maximum likelihood expectation maximization）法や、これを改良したブロック反復法による逐次近似的画像再構成技術が知られている。特に、ブロック反復法による逐次近似的画像再構成技術であるＯＳＥＭ（ordered subset ＭＬ-ＥＭ）法、ＲＡＭＬＡ（row-action maximum likelihood algorithm）法およびＤＲＡＭＡ（dynamic ＲＡＭＬＡ）法などが注目されている（非特許文献１を参照）。 Therefore, in recent years, the spread of 3D PET apparatuses is rapidly progressing, and research on image reconstruction technology based on projection data obtained by 3D PET apparatuses is also being conducted. As an image reconstruction technique in a three-dimensional PET apparatus, an ML-EM (maximum likelihood expectation maximization) method and a successive approximation image reconstruction technique based on a block iteration method improved therefrom are known. In particular, the OSEM (ordered subset ML-EM) method, the RAMLA (row-action maximum likelihood algorithm) method, the DRAMA (dynamic RAMLA) method, etc., which are successive approximate image reconstruction techniques based on the block iteration method, are attracting attention (non- (See Patent Document 1).

ＯＳＥＭ法を改良したものがＲＡＭＬＡ法であり、このＲＡＭＬＡ法を更に改良したものがＤＲＡＭＡ法である。これらＯＳＥＭ法，ＲＡＭＬＡ法およびＤＲＡＭＡ法は、何れも、投影データを複数のサブセットに分割してサブセット毎に画像を修正し、全てのサブセットについて修正を行って、これを１回の近似とし、この近似を反復することで再構成画像を得る。その反復式は下記（１）式で表される。 The RAMLA method is an improvement of the OSEM method, and the DRAMA method is a further improvement of the RAMLA method. In any of these OSEM method, RAMLA method, and DRAMA method, the projection data is divided into a plurality of subsets, the image is corrected for each subset, all the subsets are corrected, and this is regarded as one approximation. A reconstructed image is obtained by repeating the approximation. The iterative formula is represented by the following formula (1).

ここで、Ａは、所定の係数であり、ＯＳＥＭ法，ＲＡＭＬＡ法およびＤＲＡＭＡ法の何れであるかによって相違する。ｎ（＝０，１，…）は、逐次近似の反復回数である。Ｍは、サブセットの個数である。ｍ（＝０，１，２，…，Ｍ−１）は、サブセットの番号である。Ｓ_ｍは、第ｍのサブセットに含まれる投影データの集合である。ｘ_ｊ（ｊ＝１，２，…，Ｊ）は、３次元画像における第ｊ画素から放出された光子の数の平均である。ｘ_ｊ ^(n,m) は、第ｎ回近似のうちの第ｍサブセットについての修正の際に用いる第ｊ画素の値である。ｙ_ｉ（ｉ＝１，２，…，Ｉ）は、検出部における第ｉ放射線検出器対により検出された光子の数である。ａ_ｉ,ｊは, 第ｊ画素から放出された光子が第ｉ放射線検出器対により検出される確率である。 Here, A is a predetermined coefficient and differs depending on whether it is the OSEM method, the RAMLA method, or the DRAMA method. n (= 0, 1,...) is the number of iterations of successive approximation. M is the number of subsets. m (= 0, 1, 2,..., M−1) is a subset number. S _m is a set of projection data included in the m-th subset. x _j (j = 1, 2,..., J) is an average of the number of photons emitted from the j-th pixel in the three-dimensional image. x _j ^{(n, m)} is the value of the j-th pixel used when correcting the m-th subset of the n-th approximation. y _i (i = 1, 2,..., I) is the number of photons detected by the i-th radiation detector pair in the detector. a _{i, j} is the probability that the photon emitted from the j th pixel will be detected by the i th radiation detector pair.

上記（１）式において、右辺第２項の括弧内にある分数式の分母（Σａ_ｉ,ｋｘ_ｋ ^(n,m)）は、各時点での画素値ｘ_ｊ ^(n,m) に基づいて投影データを求める計算式（順投影計算）である。右辺第２項の括弧内は、順投影計算で得られた投影データと実測の投影データｙ_ｉとの差を求める計算式である。右辺第２項は、上記の差に基づいて画素値ｘ_ｊ ^(n,m) の修正量を求める計算式（逆投影計算）である。そして、右辺における第１項と第２項との加算は、上記修正量に基づいて画素値ｘ_ｊ ^(n,m) が修正された後の画素値ｘ_ｊ ^(n,m+1) を求める計算式（画像更新計算）である。つまり、上記（１）式は、順投影計算で得られる投影データと実測の投影データｙ_ｉとの比が値１に近づくように、順投影計算，逆投影計算および画像更新計算を繰り返し行って、画素値ｘ_ｊ ^(n,m) を逐次に修正していく。 In the above equation (1), the denominator (Σa _{i, k} x _k ^{(n, m)} ) in the parenthesis of the second term on the right side is based on the pixel value x _j ^{(n, m)} at each time point. This is a calculation formula for obtaining projection data (forward projection calculation). In the parentheses in the second term on the right side is a calculation formula for obtaining the difference between the projection data obtained by forward projection calculation and the actually measured projection data y _i . The second term on the right side is a calculation formula (back projection calculation) for obtaining the correction amount of the pixel value x _j ^{(n, m)} based on the above difference. The addition of the first term and the second term in the right side, obtains the pixel values x _{j ^(n,} ^m) pixels after is corrected values _{^{x j (n, m + 1}} ) based on the correction amount This is a calculation formula (image update calculation). That is, the above formula (1) is obtained by repeatedly performing forward projection calculation, back projection calculation, and image update calculation so that the ratio of the projection data obtained by forward projection calculation and the actually measured projection data y _i approaches value 1. The pixel value x _j ^{(n, m)} is corrected sequentially.

ブロック反復法による逐次近似的画像再構成では、先ず、第０回近似において、画素値ｘ_ｊ ^(0,0) および第０サブセットＳ_０の投影データを用いた計算により画素値ｘ_ｊ ^(0,1) が求められ、この画素値ｘ_ｊ ^(0,1) および第１サブセットＳ_１の投影データを用いた計算により画素値ｘ_ｊ ^(0,2) が求められ、・・・、という計算が各サブセットについて行われて、ｘ_ｊ ^(0,1)，ｘ_ｊ ^(0,2)，・・・，ｘ_ｊ ^(0,M)（＝ｘ_ｊ ^(1,0)）が逐次に求められる。次に、第１回近似において、同様の計算が各サブセットについて行われて、ｘ_ｊ ^(1,1)，ｘ_ｊ ^(1,2)，・・・，ｘ_ｊ ^(1,M)（＝ｘ_ｊ ^(2,0)）が逐次に求められる。反復回数ｎを適当な数とすることで、再構成画像を表す画素値ｘ_ｊ ^(n,0) が得られる。 In the successive approximation image reconstruction by block iteration method, first, in the zero approximation, the pixel values x _{j ^(0,0)} and the 0 subset S pixel values x _{j ⁽⁰} by calculation using the projection data of ^{_0, 1)} is obtained, and the pixel value x _j ^(0,2) is obtained by calculation using the pixel value x _j ^(0,1) and the projection data of the first subset S _1. For each subset, x _j ^(0,1) , x _j ^(0,2) ,..., X _j ^{(0, M)} (= x _j ^(1,0) ) are sequentially obtained. Next, in the first approximation, a similar calculation is performed for each subset, and x _j ^(1,1) , x _j ^(1,2) ,..., X _j ^{(1, M)} (= x _j ^(2,0) ) is obtained sequentially. By setting the number of iterations n to an appropriate number, a pixel value x _j ^{(n, 0)} representing a reconstructed image is obtained.

このように、ブロック反復法による逐次近似的画像再構成では、順投影計算，逆投影計算および画像更新計算を含む副反復（sub-iteration）をＭ回行い、このＭ回の副反復を１回の主反復（main-iteration）として、この主反復を複数回行う。したがって、ブロック反復法による逐次近似的画像再構成では、膨大な計算量および長大な計算時間を必要とする。そこで、計算時間の短縮を意図した画像再構成技術の研究が行われている。 As described above, in the successive approximation image reconstruction by the block iteration method, the sub-iteration including the forward projection calculation, the back projection calculation, and the image update calculation is performed M times, and the M sub iterations are performed once. This main iteration is performed multiple times as a main-iteration. Therefore, the successive approximation image reconstruction by the block iteration method requires a huge amount of calculation and a long calculation time. Therefore, research on image reconstruction technology intended to shorten the calculation time has been conducted.

非特許文献２に記載された画像再構成技術では、各サブセットが更に複数のグループに分割され、その分割された後の各グループについて順投影計算および逆投影計算が並列的に行われ、各グループについての逆投影計算の結果が集約された後に画像更新計算が行われる。このように、各グループについての順投影計算および逆投影計算が並列的に行われることで各々の副反復の計算時間の短縮化が図られるとする。
田中栄一、「ＰＥＴ画像の再構成法の現状と展望」、日本放射線技術学会雑誌、第６２巻、第６号、第７７１頁〜第７７７頁、２００６年６月 M. D. Jones, R. Yao, C. P. Bhole, “Hybrid MPI-OpenMP Programming forParallel OSEM PET Reconstruction”, IEEE Trans. Nucl. Sci., Vol.53, No.5,pp.2752-2758, 2006. In the image reconstruction technique described in Non-Patent Document 2, each subset is further divided into a plurality of groups, and the forward projection calculation and the back projection calculation are performed in parallel for each group after the division. The image update calculation is performed after the results of the backprojection calculation for are collected. As described above, it is assumed that the forward projection calculation and the back projection calculation for each group are performed in parallel to reduce the calculation time of each sub-repetition.
Eiichi Tanaka, “Current Status and Prospects of PET Image Reconstruction Methods”, Journal of Japanese Society of Radiological Technology, Vol. 62, No. 6, pp. 771-777, June 2006 MD Jones, R. Yao, CP Bhole, “Hybrid MPI-OpenMP Programming for Parallel OSEM PET Reconstruction”, IEEE Trans. Nucl. Sci., Vol. 53, No. 5, pp. 2752-2758, 2006.

しかしながら、非特許文献２に記載された画像再構成技術では、順投影計算および逆投影計算が並列的に行われることで計算時間の短縮化が図られるとしても、順投影計算および逆投影計算を並列処理するための複数の処理部に対するデータ転送の回数が増加してデータ転送時間が長くなる。並列処理する為の処理部の個数の増加に応じて、計算時間が短縮されるものの、データ転送時間は長くなる。したがって、並列処理する為の処理部の個数を増やしたとしても、計算時間およびデータ転送時間を含む画像再構成処理時間の短縮化には限界がある。 However, in the image reconstruction technique described in Non-Patent Document 2, even if the forward projection calculation and the back projection calculation are performed in parallel to reduce the calculation time, the forward projection calculation and the back projection calculation are performed. The number of data transfers for a plurality of processing units for parallel processing increases, and the data transfer time becomes longer. Although the calculation time is shortened as the number of processing units for parallel processing increases, the data transfer time becomes longer. Therefore, even if the number of processing units for parallel processing is increased, there is a limit to shortening the image reconstruction processing time including the calculation time and the data transfer time.

本発明は、上記問題点を解消する為になされたものであり、ブロック反復法による画像再構成の所要時間を短縮することができる画像処理装置、および、このような画像処理装置を含む３次元ＰＥＴ装置を提供することを目的とする。 The present invention has been made to solve the above problems, and an image processing apparatus capable of reducing the time required for image reconstruction by the block iteration method, and a three-dimensional image including such an image processing apparatus. An object is to provide a PET apparatus.

本発明に係る画像処理装置は、複数の放射線検出器を各々含む複数の検出器リングが軸方向に積層された３次元ＰＥＴ装置において得られた投影データを用いてブロック反復法により画像再構成を行う画像処理装置であって、(1) 投影データのうちリング差０及び１を含む部分投影データを用いて画像再構成を行って３次元の初期画像を作成する初期画像作成部と、(2) 各々、初期画像作成部により作成された初期画像から出発して、投影データの一部を用いてブロック反復法により３次元画像の再構成を行う複数の処理部と、(3) 複数の処理部それぞれによる画像再構成により得られた３次元画像を統合する画像統合部と、(4) 初期画像作成部，複数の処理部および画像統合部それぞれの動作を制御する制御部と、を備えることを特徴とする。さらに、本発明に係る画像処理装置では、制御部は、投影データがリング差および方位角によって分割された複数のサブセットのうち共通のリング差に対応するサブセット群の投影データを複数の処理部のうちの何れかの処理部に割り当て、複数の処理部は、各々に割り当てられたサブセット群の投影データを用いて順投影計算，逆投影計算および画像更新計算を繰り返すブロック反復法による画像再構成を互いに並列的に行うことを特徴とする。 An image processing apparatus according to the present invention performs image reconstruction by a block repetition method using projection data obtained in a three-dimensional PET apparatus in which a plurality of detector rings each including a plurality of radiation detectors are stacked in the axial direction. (1) an initial image creation unit that creates a three-dimensional initial image by performing image reconstruction using partial projection data including ring differences 0 and 1 among projection data; and (2 A plurality of processing units, each starting from an initial image created by the initial image creation unit, and reconstructing a three-dimensional image by a block iteration method using a part of the projection data; and (3) a plurality of processings An image integration unit that integrates three-dimensional images obtained by image reconstruction by each unit, and (4) a control unit that controls operations of the initial image creation unit, the plurality of processing units, and the image integration unit. It is characterized by. Furthermore, in the image processing apparatus according to the present invention, the control unit is configured to output projection data of a subset group corresponding to a common ring difference among a plurality of subsets obtained by dividing the projection data by the ring difference and the azimuth angle, to the plurality of processing units. Assigned to any one of the processing units, the plurality of processing units perform image reconstruction by a block iteration method that repeats forward projection calculation, back projection calculation, and image update calculation using projection data of a subset group assigned to each of the processing units. It is characterized by being performed in parallel with each other.

本発明に係る画像処理装置では、３次元ＰＥＴ装置において得られた投影データが用いられて、ブロック反復法により以下のようにして画像再構成が行われる。先ず、初期画像作成部により、投影データのうちリング差０及び１を含む部分投影データが用いられて、画像再構成が行われ３次元の初期画像が作成される。次に、制御部により、投影データがリング差および方位角によって分割された複数のサブセットのうち共通のリング差に対応するサブセット群の投影データが、複数の処理部のうちの何れかの処理部に割り当てられる。複数の処理部それぞれにより、初期画像作成部により作成された初期画像から出発して、各々に割り当てられたサブセット群の投影データが用いられて、ブロック反復法により３次元画像の再構成が互いに並列的に行われる。そして、画像統合部により、複数の処理部それぞれによる画像再構成により得られた３次元画像が統合される。これにより、並列的に処理を行う各処理部における処理時間が短縮され、かつ、各処理部に関わるデータ転送時間も短縮されるので、３次元投影データから３次元再構成画像を得るための全体の所要時間が短縮される。 In the image processing apparatus according to the present invention, the projection data obtained by the three-dimensional PET apparatus is used, and image reconstruction is performed as follows by the block iteration method. First, the initial image creation unit uses partial projection data including ring differences 0 and 1 in the projection data to perform image reconstruction and create a three-dimensional initial image. Next, the projection data of the subset group corresponding to the common ring difference among the plurality of subsets obtained by dividing the projection data by the ring difference and the azimuth by the control unit is processed by any one of the plurality of processing units. Assigned to. Starting from the initial image created by the initial image creation unit by each of the plurality of processing units, the projection data of the subset group assigned to each is used, and the reconstruction of the three-dimensional image is performed in parallel by the block iteration method. Done. Then, the three-dimensional image obtained by the image reconstruction by each of the plurality of processing units is integrated by the image integration unit. As a result, the processing time in each processing unit that performs processing in parallel is shortened, and the data transfer time related to each processing unit is also shortened. Therefore, the whole for obtaining a three-dimensional reconstructed image from three-dimensional projection data The required time is shortened.

本発明に係る画像処理装置では、初期画像作成部は、部分投影データが軸方向の位置によって分割された各データを用いてスライス画像の再構成を互いに並列的に行い、これにより得られた軸方向の各位置の画像を合成して初期画像を作成することが好ましい。或いは、初期画像作成部は、部分投影データを用いてＦＢＰ法により軸方向の各位置の２次元の画像を作成し、これにより得られた軸方向の各位置の画像を合成して初期画像を作成することも好ましい。 In the image processing apparatus according to the present invention, the initial image creation unit reconstructs slice images in parallel using each data obtained by dividing the partial projection data according to the position in the axial direction, and obtains the axis obtained thereby. It is preferable to create an initial image by combining images at respective positions in the direction. Alternatively, the initial image creation unit creates a two-dimensional image at each position in the axial direction by the FBP method using the partial projection data, and synthesizes the images at each position in the axial direction thus obtained to generate the initial image. It is also preferable to create it.

本発明に係る３次元ＰＥＴ装置は、(1) 複数の放射線検出器を各々含む複数の検出器リングが軸方向に積層された検出部と、(2) この検出部による放射線検出に基づいて得られた投影データを用いてブロック反復法により画像再構成を行う上記の本発明に係る画像処理装置と、を備えることを特徴とする。 The three-dimensional PET apparatus according to the present invention is obtained based on (1) a detection unit in which a plurality of detector rings each including a plurality of radiation detectors are stacked in the axial direction, and (2) radiation detection by the detection unit. And an image processing apparatus according to the present invention that performs image reconstruction by a block iteration method using the projection data obtained.

本発明によれば、ブロック反復法による画像再構成の所要時間を短縮することができる。 According to the present invention, the time required for image reconstruction by the block iteration method can be shortened.

以下、添付図面を参照して、本発明を実施するための最良の形態を詳細に説明する。なお、図面の説明において同一の要素には同一の符号を付し、重複する説明を省略する。 The best mode for carrying out the present invention will be described below in detail with reference to the accompanying drawings. In the description of the drawings, the same elements are denoted by the same reference numerals, and redundant description is omitted.

図１は、本実施形態に係る３次元ＰＥＴ装置１および画像処理装置２の構成を示す図である。３次元ＰＥＴ装置１は、検出部１０，信号処理部２０，収集部３０，スイッチングハブ４０，Ｎ個の処理部５０_０〜５０_Ｎ−１，制御部６０および記憶部７０を備える。ただし、Ｎは２以上の整数である。これらのうちＮ個の処理部５０_０〜５０_Ｎ−１，制御部６０および記憶部７０は、スイッチングハブ４０を介して相互に接続されていて、画像処理装置２を構成している。また、スイッチングハブ４０には収集部３０も接続されている。 FIG. 1 is a diagram illustrating a configuration of a three-dimensional PET apparatus 1 and an image processing apparatus 2 according to the present embodiment. The three-dimensional PET apparatus 1 includes a detection unit 10, a signal processing unit 20, a collection unit 30, a switching hub 40, N processing units 50 ₀ to 50 _N−1 , a control unit 60, and a storage unit 70. However, N is an integer of 2 or more. Among these, the N processing units 50 ₀ to 50 _N−1 , the control unit 60, and the storage unit 70 are connected to each other via the switching hub 40 and constitute the image processing apparatus 2. The collection unit 30 is also connected to the switching hub 40.

検出部１０は、軸方向に積層された複数の検出器リングを含む。各検出器リングは、複数の放射線検出器を含んでおり、検出器リング間にシールド板が存在しない。各放射線検出器は、検出部１０の内側の測定空間から飛来した光子を検出して、その光子エネルギに応じた値の光子検出データを出力する。検出部１０の詳細については後に図２を用いて説明する。 The detection unit 10 includes a plurality of detector rings stacked in the axial direction. Each detector ring includes a plurality of radiation detectors, and there is no shield plate between the detector rings. Each radiation detector detects a photon flying from the measurement space inside the detection unit 10 and outputs photon detection data having a value corresponding to the photon energy. Details of the detection unit 10 will be described later with reference to FIG.

信号処理部２０は、検出器１０に含まれる各放射線検出器から出力される光子検出データを入力して、この光子検出データに基づいて、電子・陽電子の対消滅に伴って発生し互いに逆方向に飛行する光子対を１対の放射線検出器が検出したか否かを判定する。そして、信号処理部２０は、１対の放射線検出器が光子対を検出したと判定したとき（すなわち、同時計数イベントが発生したと判定したとき）に、その１対の放射線検出器を識別するデータ（同時計数情報）を出力する。この１対の放射線検出器を結ぶ直線は同時計数ラインと呼ばれる。 The signal processing unit 20 receives photon detection data output from each radiation detector included in the detector 10, and is generated along with the pair annihilation of electrons and positrons based on the photon detection data and is in the opposite direction to each other. It is determined whether a pair of radiation detectors has detected a pair of photons flying in Then, the signal processing unit 20 identifies the pair of radiation detectors when it is determined that the pair of radiation detectors has detected the photon pair (that is, when it is determined that a coincidence event has occurred). Output data (simultaneous counting information). A straight line connecting the pair of radiation detectors is called a coincidence line.

収集部３０は、信号処理部２０から同時計数情報を受け取り、その同時計数情報を記憶部７０に保存させる。Ｎ個の処理部５０_０〜５０_Ｎ−１は、記憶部７０により保存された同時計数情報の集合である投影データを用いてブロック反復法により画像再構成を行う。制御部６０は、３次元ＰＥＴ装置１の全体の動作を制御する。また、制御部６０は、画像再構成処理の一部の計算を担ってもよい。 The collection unit 30 receives the coincidence counting information from the signal processing unit 20 and causes the storage unit 70 to store the coincidence counting information. The N processing units 50 ₀ to 50 _N−1 perform image reconstruction by the block iteration method using projection data that is a set of coincidence counting information stored in the storage unit 70. The control unit 60 controls the overall operation of the three-dimensional PET apparatus 1. In addition, the control unit 60 may perform a part of the image reconstruction process.

図２は、本実施形態に係る３次元ＰＥＴ装置１の検出部１０の断面図である。この図は、中心軸を含む面で検出部１０を切断したときの断面を示している。３次元ＰＥＴ装置１の検出部１０の構成は、シールド１１とシールド１２との間に一定ピッチで順に積層された検出器リングＲ_１〜Ｒ_７を有している。検出器リングＲ_１〜Ｒ_７それぞれは、中心軸に垂直な面上にリング状に配された複数の放射線検出器を有している。各放射線検出器は、例えばＢＧＯ（Ｂｉ_４Ｇｅ_３Ｏ_１２）等のシンチレータと光電子増倍管とを組み合わせたシンチレーション検出器であり、中心軸を含む測定空間から飛来して到達した光子を検出する。２次元ＰＥＴ装置の場合と異なり、３次元ＰＥＴ装置１ではスライスセプタが備えられていない。このように構成される３次元ＰＥＴ装置１の検出部１０は、あらゆる方向から飛来した光子対を同時計数することができる。すなわち、３次元ＰＥＴ装置１の検出部１０により得られ蓄積される同時計数情報は、任意の検出器リングに含まれる１対の放射線検出器によるものが可能である。 FIG. 2 is a cross-sectional view of the detection unit 10 of the three-dimensional PET apparatus 1 according to the present embodiment. This figure has shown the cross section when the detection part 10 is cut | disconnected by the surface containing a central axis. The configuration of the detection unit 10 of the three-dimensional PET apparatus 1 includes detector rings R _{1 to} R ₇ that are sequentially stacked at a constant pitch between the shield 11 and the shield 12. Each of the detector rings R _{1 to} R ₇ has a plurality of radiation detectors arranged in a ring shape on a plane perpendicular to the central axis. Each radiation detector is a scintillation detector that combines a scintillator such as BGO (Bi ₄ Ge ₃ O ₁₂ ) and a photomultiplier tube, and detects photons that have arrived from a measurement space including the central axis. . Unlike the case of the two-dimensional PET apparatus, the three-dimensional PET apparatus 1 is not provided with a slice ceptor. The detection unit 10 of the three-dimensional PET apparatus 1 configured as described above can simultaneously count photon pairs flying from all directions. That is, the coincidence information obtained and accumulated by the detection unit 10 of the three-dimensional PET apparatus 1 can be obtained by a pair of radiation detectors included in an arbitrary detector ring.

３次元ＰＥＴ装置１における同時計数ラインすなわち同時計数情報は、４つの変数ｒ，θ，ｚ，δにより表される。変数ｒは、中心軸からの同時計数ラインの距離を表す。変数θは、同時計数ラインの方位角を表す。変数ｚは、同時計数ラインの中点の中心軸方向位置を表す。また、変数δは、光子対を検出した各放射線検出器が含まれる検出器リングの間の中心軸方向距離（リング差）を表す。３次元の投影データは、４つの変数ｒ，θ，ｚ，δに対して、同時計数情報を取得した頻度（同時計数イベントの発生頻度）を表すヒストグラムとして表現される。記憶部７０は、このような投影データを保存する。 A coincidence line in the three-dimensional PET apparatus 1, that is, coincidence information is represented by four variables r, θ, z, and δ. The variable r represents the distance of the coincidence line from the central axis. The variable θ represents the azimuth angle of the coincidence counting line. The variable z represents the central axis direction position of the midpoint of the coincidence counting line. Further, the variable δ represents a central axis direction distance (ring difference) between detector rings including each radiation detector that detects a photon pair. The three-dimensional projection data is expressed as a histogram representing the frequency of acquiring coincidence count information (occurrence frequency of coincidence count events) for the four variables r, θ, z, and δ. The storage unit 70 stores such projection data.

以下では、検出部１０に含まれる検出器リングの個数は、処理部５０_０〜５０_Ｎ−１の個数と等しいＮであるとする。このとき、リング差δは、０からＮ−１までのＮ個の値をとり得る。また、検出部１０の中心軸方向に垂直なスライスの個数は２Ｎ−１個であり、そのうち、リング差０のデータから得られるダイレクトスライスの個数はＮ個であり、リング差１のデータから得られるクロススライスの個数は(Ｎ−１)個である。 Hereinafter, it is assumed that the number of detector rings included in the detection unit 10 is _N equal to the number of processing units 50 ₀ to 50 _N−1 . At this time, the ring difference δ can take N values from 0 to N−1. In addition, the number of slices perpendicular to the central axis direction of the detection unit 10 is 2N−1, of which the number of direct slices obtained from data with a ring difference of 0 is N and obtained from data with a ring difference of 1. The number of cross slices to be obtained is (N-1).

本実施形態においては、３次元投影データは、リング差δおよび方位角θによって複数のサブセットＳ_δ,θに分割される。リング差δがＮ_δとおりの値をとり、方位角θがＮ_θとおりの値をとるとすれば、３次元投影データはＮ_δＮ_θ個のサブセットに分割される。サブセットＳ_δ,θには、３次元投影データのうち特定のリング差δおよび特定の方位角θの投影データが含まれる。 In the present embodiment, the three-dimensional projection data is divided into a plurality of subsets S _{δ, θ} by the ring difference δ and the azimuth angle θ. If the ring difference δ takes N _δ values and the azimuth angle θ takes N _θ values, the three-dimensional projection data is divided into N _δ N _θ subsets. The subset S _{δ, θ} includes projection data of a specific ring difference δ and a specific azimuth angle θ among the three-dimensional projection data.

また、本実施形態においては、共通のリング差δに対応するサブセットＳ_δ,θの集合は、サブセット群ＳＧ_δと表される。３次元投影データはＮ_δ個のサブセット群に分割される。サブセット群ＳＧ_δには、特定のリング差δであって全ての方位角のＮ_θ個のサブセットＳ_δ,θ、すなわち、３次元投影データのうち特定のリング差δの投影データが含まれる。 In the present embodiment _, a set of subsets S _{δ, θ} corresponding to a common ring difference _δ is represented as a subset group SG _δ . The three-dimensional projection data is divided into N _δ subset groups. The subset group SG _δ includes a specific ring difference δ and N _θ subsets S _{δ, θ of} all azimuth angles, that is, projection data of a specific ring difference δ among the three-dimensional projection data.

そして、画像処理装置２は、Ｎ_θＮ_δ個のサブセットまたはＮ_δ個のサブセット群に分割される３次元投影データを用いて、ブロック反復法により画像再構成を行う。図３〜５は、本実施形態に係る画像処理装置２における画像再構成処理のフローチャートである。図３に示されるように、画像処理装置２は、３次元の初期画像を作成し（ステップＳ１０）、この作成した３次元の初期画像から出発して各サブセット群ＳＧ_δの投影データを用いて画像生構成処理して（ステップＳ２０_０〜Ｓ２０_Ｎ−１）、これにより得られた各々の３次元画像を統合する（ステップＳ３０）。 Then, the image processing apparatus 2 performs image reconstruction by the block iteration method using the three-dimensional projection data divided into N _θ N _δ subsets or N _δ subset groups. 3 to 5 are flowcharts of image reconstruction processing in the image processing apparatus 2 according to the present embodiment. As shown in FIG. 3, the image processing apparatus 2 creates a three-dimensional initial image (step S10), and uses the projection data of each subset group SG _δ starting from the created three-dimensional initial image. Image raw composition processing is performed (steps S20 _{0 to} S20 _N-1 ), and the respective three-dimensional images obtained thereby are integrated (step S30).

初期画像作成処理（ステップＳ１０）では、３次元投影データのうちリング差０及び１を含む部分投影データを用いて画像再構成を行って３次元の初期画像を作成する。このとき用いられる部分投影データは、リング差０及び１だけでなくリング差２以上を含んでいてもよい。初期画像作成処理（ステップＳ１０）を行う初期画像作成部として、処理部５０_０〜５０_Ｎ−１または制御部６０が用いられる。初期画像作成処理（ステップＳ１０）の詳細については後に図４を用いて説明する。 In the initial image creation process (step S10), image reconstruction is performed using partial projection data including ring differences 0 and 1 among the three-dimensional projection data to create a three-dimensional initial image. The partial projection data used at this time may include not only the ring differences 0 and 1, but also the ring difference 2 or more. The processing units 50 ₀ to 50 _N-1 or the control unit 60 is used as an initial image creating unit that performs the initial image creating process (step S10). Details of the initial image creation processing (step S10) will be described later with reference to FIG.

各サブセット群ＳＧ_δの投影データを用いた画像再構成処理（ステップＳ２０_δ）では、ステップＳ１０で作成された３次元の初期画像から出発して、サブセット群ＳＧ_δの投影データを用いてブロック反復法により３次元画像の再構成を行う。これに際して、制御部６０は、記憶部７０に保存されている３次元投影データのうちサブセット群ＳＧ_δを処理部５０_δに割り当てる。そして、各処理部５０_δは、割り当てられたサブセット群ＳＧ_δの投影データを用いて画像再構成を互いに並列的に行う。すなわち、Ｎ個の処理部５０_０〜５０_Ｎ−１は互いに並列的に動作して、ステップＳ２０_０〜Ｓ２０_Ｎ−１は互いに並列的に行われる。この画像再構成処理（ステップＳ２０_δ）の詳細については後に図５を用いて説明する。 In the image reconstruction process using the projection data of each subset group SG _δ (step S20 _δ ), starting from the three-dimensional initial image created in step S10, block iteration is performed using the projection data of the subset group SG _δ. The 3D image is reconstructed by the method. At this time, the control unit 60 assigns the subset group SG _δ among the three-dimensional projection data stored in the storage unit 70 to the processing unit 50 _δ . Then, each processing unit 50 _δ performs image reconstruction in parallel with each other using the projection data of the assigned subset group SG _δ . That is, the N processing units 50 ₀ to 50 _N-1 operate in parallel with each other, and steps S20 _{0 to} S20 _N-1 are performed in parallel with each other. Details of this image reconstruction process (step S20 _δ ) will be described later with reference to FIG.

なお、リング差０，１のサブセット群の投影データについては、初期画像作成処理（ステップＳ１０）において既に処理しているので、画像再構成処理（ステップＳ２０_δ）において再度の処理を行わないこととしてもよい。 Note that the projection data subset group of rings difference 0,1, since in the initial image generation processing (step S10) has already been treated, as is not performed again processed by the image reconstruction processing (step S20 _[delta]) Also good.

統合処理（ステップＳ３０）では、Ｎ個の処理部５０_０〜５０_Ｎ−１それぞれによる画像再構成で得られた３次元画像を統合する。ここで、統合とは、Ｎ個の処理部５０_０〜５０_Ｎ−１により得られたＮ個の３次元画像を単純に加算すること、Ｎ個の３次元画像を加算平均すること、Ｎ個の３次元画像を加重平均すること、等の何れであってもよい。統合処理（ステップＳ３０）を行う統合部として処理部５０_０〜５０_Ｎ−１または制御部６０が用いられる。 In the integration process (step S30), three-dimensional images obtained by image reconstruction by each of the _N processing units 50 ₀ to 50 _N−1 are integrated. Here, integration refers to simply adding N three-dimensional images obtained by _N processing units 50 ₀ to 50 _N−1 , adding and averaging N three-dimensional images, N The three-dimensional image may be weighted averaged. The processing units 50 ₀ to 50 _N-1 or the control unit 60 is used as an integration unit that performs the integration process (step S30).

なお、本実施形態では各処理部５０_δが１つのサブセット群ＳＧ_δの投影データを処理することとしたが、各処理部５０_δが処理するサブセット群の個数が一定でない場合には、このことを考慮して加重平均をするのが好適である。また、各処理部５０_δにより得られた再構成画像の統計精度が大きく異なっている場合にも、このことを考慮して加重平均をするのが好適である。また、必要に応じて、統合処理（ステップＳ３０）で得られた３次元画像をスライス面および中心軸方向に平滑化して雑音を低減することが好ましい。 In the present embodiment, each processing unit 50 _δ processes projection data of one subset group SG _δ . However, if the number of subset groups processed by each processing unit 50 _δ is not constant, this is the case. It is preferable to take a weighted average in consideration of the above. Further, even when the statistical accuracy of the reconstructed images obtained by the respective processing units 50 _δ is greatly different, it is preferable to perform the weighted average in consideration of this fact. Moreover, it is preferable to reduce noise by smoothing the three-dimensional image obtained by the integration process (step S30) in the slice plane and the central axis direction as necessary.

このようにして統合処理（ステップＳ３０）により得られた３次元画像は、検出部１０の内部の測定空間における光子対の発生頻度の空間分布を表す画像を再構成したものとなる。 Thus, the three-dimensional image obtained by the integration process (step S30) is a reconstructed image representing the spatial distribution of the occurrence frequency of photon pairs in the measurement space inside the detection unit 10.

図４は、本実施形態に係る画像処理装置２における画像再構成処理のうち初期画像作成処理（ステップＳ１０）のフローチャートである。３次元の初期画像を作成する初期画像作成処理（ステップＳ１０）は、３次元投影データのうちリング差０及び１を含む部分投影データを用いて各スライスの画像を再構成する処理（ステップＳ１１_０〜Ｓ１１_Ｎ−１）、各々のスライス画像を合成する処理（ステップＳ１２）、および、合成により得られた画像を平滑化する処理（ステップＳ１３）を含む。 FIG. 4 is a flowchart of the initial image creation process (step S10) in the image reconstruction process in the image processing apparatus 2 according to the present embodiment. The initial image creation process (step S10) for creating a three-dimensional initial image reconstructs an image of each slice using partial projection data including ring differences 0 and 1 in the three-dimensional projection data (step S11 _0). -S11N _-1 ), the process (step S12) which synthesize | combines each slice image, and the process (step S13) which smoothes the image obtained by the synthesis | combination.

スライス画像再構成処理は、互いに並列的に行われるＮ個の処理に分けられる（ステップＳ１１_０〜Ｓ１１_Ｎ−１）。ステップＳ１１_０では、処理部５０_０により、スライス０（ダイレクトスライス）のスライス画像、および、スライス１（クロススライス）のスライス画像が再構成される。ステップＳ１１_１では、処理部５０_１により、スライス２（ダイレクトスライス）のスライス画像、および、スライス３（クロススライス）のスライス画像が再構成される。ステップＳ１１_２〜Ｓ１１_Ｎ−２も同様である。ステップＳ１１_Ｎ−１では、処理部５０_Ｎ−１により、スライス２Ｎ−１（ダイレクトスライス）のスライス画像が再構成される。 The slice image reconstruction process is divided into N processes performed in parallel with each other (steps S11 _{0 to} S11 _N-1 ). In step S11 _0, the processing unit 50 _0, the slice image of the slice 0 (direct slices), and is reconstructed slice image of the slice 1 (cross slice). In step S11 _1, the processing unit 50 _1, the slice image of the slice 2 (Direct slices), and reconstructed slice images of slice 3 (cross-slice). The same applies to steps S11 _{2 to} S11 _N-2 . In step S11 _N-1 , the processing unit 50 _N-1 reconstructs a slice image of the slice 2N-1 (direct slice).

なお、これらステップＳ１１_０〜Ｓ１１_Ｎ−１では、３次元投影データのうちリング差０及び１を含む部分投影データが軸方向の位置によって分割された各データを用いてスライス画像の再構成を互いに並列的に行うものであるが、これに替えて、部分投影データを用いて２次元ＦＢＰ法により軸方向の各位置の２次元の画像を作成してもよい。 Note that in these steps S11 _{0 to} S11 _N−1 , reconstruction of slice images is performed using each data obtained by dividing partial projection data including ring differences 0 and 1 among the three-dimensional projection data according to the positions in the axial direction. Although it is performed in parallel, instead of this, a two-dimensional image of each position in the axial direction may be created by the two-dimensional FBP method using partial projection data.

２次元ＦＢＰ法により軸方向の各位置の２次元の画像を作成する場合には、解析的に画像を再構成することができ、また、高速に処理をすることができるが、逐次近似法と比較して統計雑音が大きく、また、再構成された画像に負の値を含むことが多い。これを逐次近似法の初期画像とするためには、画像を充分に平滑して統計雑音を低減した後、ゼロおよび負の画素値を小さな正の値に修正しておく必要がある。 When creating a two-dimensional image of each position in the axial direction by the two-dimensional FBP method, the image can be reconstructed analytically and can be processed at high speed. In comparison, the statistical noise is large, and the reconstructed image often includes a negative value. In order to make this an initial image of the successive approximation method, it is necessary to correct the zero and negative pixel values to small positive values after sufficiently smoothing the image to reduce statistical noise.

合成処理（ステップＳ１２）では、ステップＳ１１_０〜Ｓ１１_Ｎ−１で得られたスライス０〜２Ｎ−１の各スライス画像を合成して、３次元画像を作成する。続く平滑化処理（ステップＳ１３）では、合成により得られた３次元画像をスライス面および中心軸方向に平滑化して、以降の処理に用いられる３次元の初期画像を作成する。これら合成処理（ステップＳ１２）および平滑化処理（ステップＳ１３）は、処理部５０_０〜５０_Ｎ−１または制御部６０により行われる。 In the synthesis process (step S12), the slice images of slices _{0 to} 2N _-1 obtained in steps S11 _{0 to} S11 _N-1 are synthesized to create a three-dimensional image. In the subsequent smoothing process (step S13), the three-dimensional image obtained by the synthesis is smoothed in the slice plane and the central axis direction, and a three-dimensional initial image used for the subsequent processes is created. The synthesis process (step S12) and the smoothing process (step S13) are performed by the processing units 50 ₀ to 50 _N-1 or the control unit 60.

スライス画像再構成処理（ステップＳ１１_０〜Ｓ１１_Ｎ−１）は２次元ＰＥＴ装置における画像再構成処理と同様の処理であるが、合成処理（ステップＳ１２）により全体として３次元の再構成画像を得ることができる。ただし、全体の３次元投影データのうちの一部の投影データのみが計算に用いられるので、合成処理（ステップＳ１２）により得られる３次元画像には統計雑音が多いという欠点がある。そこで、この統計雑音が後の画像再構成処理（ステップＳ２０_δ）に与える影響を低減するために、平滑化処理（ステップＳ１３）において３次元画像を平滑化しておくことが好ましい。 The slice image reconstruction process (steps S11 _{0 to} S11 _N-1 ) is the same process as the image reconstruction process in the two-dimensional PET apparatus, but a three-dimensional reconstruction image is obtained as a whole by the synthesis process (step S12). be able to. However, since only a part of the entire three-dimensional projection data is used for calculation, the three-dimensional image obtained by the synthesizing process (step S12) has a drawback that there are many statistical noises. Therefore, in order to reduce the influence of this statistical noise on the subsequent image reconstruction process (step S20 _δ ), it is preferable to smooth the three-dimensional image in the smoothing process (step S13).

図５は、本実施形態に係る画像処理装置２における画像再構成処理のうちサブセット群ＳＧ_δの投影データを用いた画像再構成処理（ステップＳ２０_δ）のフローチャートである。この図は、図３中のステップＳ２０_０〜Ｓ２０_Ｎ−１それぞれの処理のフローチャートを示すものである。画像再構成処理（ステップＳ２０_δ）は、上記のようにして作成された３次元の初期画像から出発して、特定のリング差δであって全ての方位角のＮ_θ個のサブセットＳ_δ,θ（すなわち、３次元投影データのうち特定のリング差δの投影データ）を含むサブセット群ＳＧ_δの投影データを用いて、ブロック反復法により３次元画像の再構成を行う処理である。ブロック反復法として、ＯＳＥＭ法，ＲＡＭＬＡ法およびＤＲＡＭＡ法の何れが採用されてもよい。 FIG. 5 is a flowchart of the image reconstruction process (step S20 _δ ) using the projection data of the subset group SG _δ among the image reconstruction processes in the image processing apparatus 2 according to the present embodiment. This figure shows the flowchart of each process of step S20 ₀ -S20 _N-1 in FIG. The image reconstruction process (step S20 _δ ) starts from the three-dimensional initial image created as described above, and has a specific ring difference δ and N _θ subsets S _{δ, of} all azimuth angles _{. This} is a process for reconstructing a three-dimensional image by the block iteration method using projection data of a subset group SG _δ including _θ (that is, projection data of a specific ring difference δ among the three-dimensional projection data). Any of the OSEM method, the RAMLA method, and the DRAMA method may be adopted as the block iteration method.

画像再構成処理（ステップＳ２０_δ）は、副反復回数ｍを初期化する処理（ステップＳ２１）、順投影計算を行う処理（ステップＳ２２）、逆投影計算を行う処理（ステップＳ２３）、画像更新計算を行う処理（ステップＳ２４）、副反復回数ｍを判断する処理（ステップＳ２５）および副反復回数ｍに値１を加算する処理（ステップＳ２６）を含む。 The image reconstruction process (step S20 _δ ) includes a process for initializing the sub-repetition count m (step S21), a process for performing forward projection calculation (step S22), a process for performing backprojection calculation (step S23), and an image update calculation. (Step S24), a process for determining the number of sub-repetitions m (step S25), and a process for adding the value 1 to the number of sub-repetitions m (step S26).

副反復回数初期化処理（ステップＳ２１）では、副反復回数ｍを値０に初期化する。順投影計算処理（ステップＳ２２）では前述の（１）式のうちの順投影計算を行い、逆投影計算処理（ステップＳ２３）では前述の（１）式のうちの逆投影計算を行い、また、画像更新計算処理（ステップＳ２４）では前述の（１）式のうちの画像更新計算を行う。これにより、順投影計算，逆投影計算および画像更新計算を含む副反復が１回行われて、各画素値が修正される。 In the sub iteration number initialization process (step S21), the sub iteration number m is initialized to 0. In the forward projection calculation process (step S22), the forward projection calculation of the above-described expression (1) is performed, and in the back-projection calculation process (step S23), the back-projection calculation of the above-described expression (1) is performed. In the image update calculation process (step S24), the image update calculation in the above-described equation (1) is performed. Thereby, the sub-iteration including forward projection calculation, back projection calculation, and image update calculation is performed once, and each pixel value is corrected.

副反復回数判断処理（ステップＳ２５）では副反復回数ｍが (Ｎ_θ−１)より小さいか否かを判断し、副反復回数ｍが (Ｎ_θ−１)より小さいときには副反復回数加算処理（ステップＳ２６）で副反復回数ｍに値１を加算してステップＳ２２に進み、副反復回数ｍが (Ｎ_θ−１)に等しいときには終了する。 In the sub iteration number determination process (step S25), it is determined whether or not the sub iteration number m is smaller than (N _θ −1). When the sub iteration number m is smaller than (N _θ −1), the sub iteration number addition process ( In step S26), a value of 1 is added to the number of sub-repetitions m, and the process proceeds to step S22. When the number of sub-repetitions m is equal to (N _θ −1), the process ends.

このようにして、サブセット群ＳＧ_δに含まれるサブセットＳ_δ,θの個数Ｎ_θだけ副反復が繰り返される。Ｎ_θ回の副反復のうちの第ｍ回の副反復における順投影計算処理（ステップＳ２２）および逆投影計算処理（ステップＳ２３）では、サブセット群ＳＧ_δに含まれる特定のサブセットＳ_δ,ｍの投影データが用いられる。なお、このＮ_θ回の副反復を１回の主反復として、この主反復を複数回行ってもよい。 In this way, the sub-repetition is repeated by the number N _θ of the subsets S _{δ, θ} included in the subset group SG _δ . In the forward projection calculation process (step S22) and backprojection calculation process (step S23) in the mth sub-repetition among N _θ sub-repetitions, a specific subset S _{δ, m} included in the subset group SG _δ Projection data is used. The N _θ sub-repetitions may be regarded as one main repetition, and the main repetition may be performed a plurality of times.

そして、統合処理（ステップＳ３０）では、各ステップＳ２０_δで処理部５０_δによる画像再構成で得られたＮ個の３次元画像を統合する。このようにして統合処理（ステップＳ３０）により得られた３次元画像は、検出部１０の内部の測定空間における光子対の発生頻度の空間分布を表す画像を再構成したものとなる。また、この３次元再構成画像は統計雑音が少ないものとなる。 In the integration process (step S30), N three-dimensional images obtained by image reconstruction by the processing unit 50 _δ are integrated in each step S20 _δ . Thus, the three-dimensional image obtained by the integration process (step S30) is a reconstructed image representing the spatial distribution of the occurrence frequency of photon pairs in the measurement space inside the detection unit 10. Further, this three-dimensional reconstructed image has little statistical noise.

なお、図３および図５に示されたフローチャートでは、各処理部が画像再構成処理（ステップＳ２０_δ）において１つのサブセット群を処理することとしたが、何れかの処理部または全ての処理部が２以上のサブセット群の投影データを処理することとしてもよい。例えば、リング差δが１より大きいサブセット群の投影データを用いた画像再構成処理の所要時間はリング差δに略比例するので、各処理部の計算時間が均等になるようにサブセット群の投影データが各処理部に割り当てられることが好ましい。 In the flowcharts shown in FIG. 3 and FIG. 5, each processing unit processes one subset group in the image reconstruction processing (step S20 _δ ). However, projection data of two or more subset groups may be processed. For example, the time required for the image reconstruction processing using the projection data of the subset group in which the ring difference δ is greater than 1 is substantially proportional to the ring difference δ, so that the projection of the subset group is performed so that the calculation time of each processing unit is equal. Preferably, data is assigned to each processing unit.

或る処理部が２以上のＮ_SG個のサブセット群の投影データを処理する場合には、その処理部は、１主反復のうちにＮ_SGＮ_θ個のサブセットの処理を直列的に行うことになる。このようにすることで画像更新の回数が増えるので、画像再構成の収束性が向上する。また、主反復を複数回行うことによっても画像再構成の収束性が向上する。 When a certain processing unit processes projection data of two or more N _SG subset groups, the processing unit performs processing of N _SG N _θ subsets in series in one main iteration. become. By doing so, the number of times of image update increases, so that the convergence of image reconstruction is improved. Further, the convergence of the image reconstruction is improved by performing the main iteration a plurality of times.

また、３次元投影データをリング差δまたは方位角θに関してスパニング処理を行って、その処理後の３次元投影データをサブセットまたはサブセット群に分割してもよい。３次元投影データをリング差δに関してスパニング処理を行う場合の一例として、リング差０〜３の投影データをリング差０のセグメントの投影データとして纏め、リング差４〜１０の投影データをリング差７のセグメントの投影データとして纏め、リング差８〜１７の投影データをリング差１４のセグメントの投影データとして纏める、という処理を行う。このようなスパニング処理を行った後の３次元投影データを用いて、初期画像作成処理（ステップＳ１０），画像再構成処理（ステップＳ２０_δ）および統合処理（ステップＳ３０）を行ってもよい。この場合、初期画像作成処理（ステップＳ１０）ではリング差０のセグメントの投影データを用いて画像再構成を行って３次元の初期画像を作成し、画像再構成処理（ステップＳ２０_δ）では各リング差セグメントの投影データを用いてブロック反復法により３次元画像の再構成を行う。 Further, spanning processing may be performed on the three-dimensional projection data with respect to the ring difference δ or the azimuth angle θ, and the three-dimensional projection data after the processing may be divided into subsets or subset groups. As an example of performing the spanning process on the three-dimensional projection data with respect to the ring difference δ, the projection data of the ring difference 0 to 3 are collected as the projection data of the segment of the ring difference 0, and the projection data of the ring difference 4 to 10 is combined with the ring difference 7 And the projection data of the ring difference 8 to 17 are collected as the projection data of the segment of the ring difference 14. The initial image creation process (step S10), the image reconstruction process (step S20 _δ ), and the integration process (step S30) may be performed using the three-dimensional projection data after performing such a spanning process. In this case, in the initial image creation process (step S10), image reconstruction is performed using the projection data of the segment with the ring difference 0 to create a three-dimensional initial image, and in the image reconstruction process (step S20 _δ ), each ring A three-dimensional image is reconstructed by the block iteration method using the projection data of the difference segment.

以上のように構成されて画像再構成処理を行う本実施形態に係る３次元ＰＥＴ装置１または画像処理装置２は、初期画像作成処理（ステップＳ１０），サブセット群ＳＧ_δの投影データを用いた画像再構成処理（ステップＳ２０_δ）および統合処理（ステップＳ３０）のうち、最も長い処理時間を要する画像再構成処理（ステップＳ２０_δ）を処理部５０_δにより並列的に処理することができる。 The three-dimensional PET apparatus 1 or the image processing apparatus 2 according to the present embodiment configured as described above and performing image reconstruction processing is an initial image creation process (step S10), an image using projection data of the subset group SG _δ. Of the reconstruction process (step S20 _δ ) and the integration process (step S30), the image reconstruction process (step S20 _δ ) that requires the longest processing time can be processed in parallel by the processing unit 50 _δ .

また、本実施形態に係る３次元ＰＥＴ装置１または画像処理装置２では、並列処理する各処理部５０_δに関わるデータ転送は、初期画像作成処理（ステップＳ１０）で得られた３次元の初期画像のデータの各処理部５０_δへの転送と、サブセット群ＳＧ_δの投影データの各処理部５０_δへの転送と、画像再構成処理（ステップＳ２０_δ）で得られた再構成画像のデータの各処理部５０_δからの転送と、を含む。しかし、画像再構成処理（ステップＳ２０_δ）における各処理部５０_δによる処理の際には、副反復回数および主反復回数が幾ら多くとも、他の処理部および制御部６０との間でデータ転送の必要がない。 In the three-dimensional PET apparatus 1 or the image processing apparatus 2 according to the present embodiment, the data transfer related to each processing unit 50 _δ to be processed in parallel is the three-dimensional initial image obtained in the initial image creation process (step S10). of the transfer to the processing unit 50 _[delta] of the data, and transfer to the processing unit 50 _[delta] of projection data subsets group SG _[delta], the image reconstruction processing of the data of the obtained reconstructed image (step S20 _[delta]) Transfer from each processing unit 50 _δ . However, during the processing by each processing unit 50 _δ in the image reconstruction processing (step S20 _δ ), data transfer is performed between the other processing units and the control unit 60 regardless of the number of sub-repetitions and main iterations. There is no need for.

したがって、本実施形態に係る３次元ＰＥＴ装置１または画像処理装置２は、各処理部５０_δおよび制御部６０における処理時間が短縮され、かつ、各処理部５０_δおよび制御部６０の間におけるデータ転送時間も短縮されるので、３次元投影データから３次元再構成画像を得るための全体の所要時間が短縮される。 Therefore, in the three-dimensional PET apparatus 1 or the image processing apparatus 2 according to the present embodiment, the processing time in each processing unit 50 _δ and the control unit 60 is reduced, and data between each processing unit 50 _δ and the control unit 60 is reduced. Since the transfer time is also shortened, the total time required to obtain a three-dimensional reconstructed image from the three-dimensional projection data is shortened.

次に、本実施形態に係る３次元ＰＥＴ装置１または画像処理装置２における画像再構成処理の計算機シミュレーション結果について、図６〜図９を用いて説明する。 Next, computer simulation results of image reconstruction processing in the three-dimensional PET apparatus 1 or the image processing apparatus 2 according to the present embodiment will be described with reference to FIGS.

図６〜図８では、複雑なＲＩ線源分布を有するファントムを仮定して、統計雑音がない条件で再構成した画像の誤差（ファントムにおけるＲＩ線源分布からの二乗平均誤差）を求めて、これを再構成誤差とした。この再構成誤差は収束の度合いを示すものである。検出器リングの直径を８０ｃｍとした。検出器リングピッチを８ｍｍとした。画像マトリックスのサイズを１２８×１２８とした。画素サイズを４ｍｍとした。方位角数Ｎ_θを１２８とした。また、最終画像の解像度を安定させるための平滑フィルタ（ポストスムージング）の半値幅を２画素とした。 6 to 8, assuming a phantom having a complicated RI source distribution, an error of a reconstructed image under a condition without statistical noise (root mean square error from the RI source distribution in the phantom) is obtained. This was regarded as a reconstruction error. This reconstruction error indicates the degree of convergence. The diameter of the detector ring was 80 cm. The detector ring pitch was 8 mm. The size of the image matrix was 128 × 128. The pixel size was 4 mm. The azimuth number _Nθ was 128. Further, the half-value width of the smoothing filter (post smoothing) for stabilizing the resolution of the final image is set to 2 pixels.

また、図６〜図８それぞれで、para-DRAMA(s) は、各処理部が１つのサブセット群の投影データを処理する場合の結果を示す。para-DRAMA(t) は、各処理部が２つのサブセット群の投影データを処理する場合の結果を示す。これら para-DRAMA(s) および para-DRAMA(t) は本実施形態に係るものである。3D-DRAMA は、並列処理を行わない従来のＤＲＡＭＡ法により処理する場合の結果を示す。 Further, in each of FIGS. 6 to 8, para-DRAMA (s) indicates a result when each processing unit processes projection data of one subset group. para-DRAMA (t) indicates a result when each processing unit processes projection data of two subset groups. These para-DRAMA (s) and para-DRAMA (t) relate to this embodiment. 3D-DRAMA shows the result of processing by the conventional DRAMA method without parallel processing.

図６は、リング差δの最大値と再構成画像の解像度（半値全幅）との関係を示すグラフである。この図に示されるように、para-DRAMA(s) および para-DRAMA(t) それぞれの場合の再構成画像の解像度は、3D-DRAMA の場合の再構成画像の解像度と比べて幾分低いが、その差は僅かである。したがって、統合処理（ステップＳ３０）の後の平滑化処理の際の半値幅により解像度を調整することで、para-DRAMA(s) および para-DRAMA(t) それぞれの場合においても充分な解像度を有する再構成画像を得ることができる。 FIG. 6 is a graph showing the relationship between the maximum value of the ring difference δ and the resolution (full width at half maximum) of the reconstructed image. As shown in this figure, the resolution of the reconstructed image in each of para-DRAMA (s) and para-DRAMA (t) is somewhat lower than that of 3D-DRAMA. The difference is slight. Therefore, by adjusting the resolution according to the half width at the time of the smoothing process after the integration process (step S30), sufficient resolution is obtained in each case of para-DRAMA (s) and para-DRAMA (t). A reconstructed image can be obtained.

図７は、リング差δの最大値と再構成画像の再構成誤差との関係を示すグラフである。この図に示されるように、para-DRAMA(s) および para-DRAMA(t) それぞれの場合の再構成画像の再構成誤差は、3D-DRAMA の場合の再構成画像の再構成誤差と比べて大きいが、実用上では何ら問題がない。para-DRAMA(t) の場合の再構成誤差は、para-DRAMA(s) の場合の再構成誤差と比べて小さい。これは、para-DRAMA(t) の場合には、各処理部が２つのサブセット群の投影データを処理するので、画像更新回数が２倍になり、収束が進んだことを意味している。収束についてはpara-DRAMA(t) で実用上十分であるが、もし必要があれば、各処理部が処理するサブセット群の個数または主反復回数を増大することによって、更に充分な収束を得ることも可能である。3D-DRAMA の場合には、最大リング差の増加とともに再構成誤差は小さくなっている。これは、3D-DRAMA の場合には、最大リング差に比例して画像更新回数が増大することによるが、その代わり並列演算ができない欠点がある。 FIG. 7 is a graph showing the relationship between the maximum value of the ring difference δ and the reconstruction error of the reconstructed image. As shown in this figure, the reconstruction error of the reconstructed image in each case of para-DRAMA (s) and para-DRAMA (t) is compared with the reconstruction error of the reconstructed image in 3D-DRAMA. Large, but no problem in practical use. The reconstruction error in the case of para-DRAMA (t) is small compared to the reconstruction error in the case of para-DRAMA (s). This means that in the case of para-DRAMA (t), since each processing unit processes projection data of two subset groups, the number of image updates has doubled and convergence has advanced. For convergence, para-DRAMA (t) is practically sufficient, but if necessary, more convergence can be obtained by increasing the number of subsets or the number of main iterations processed by each processing unit. Is also possible. In the case of 3D-DRAMA, the reconstruction error decreases as the maximum ring difference increases. This is because in 3D-DRAMA, the number of image updates increases in proportion to the maximum ring difference.

また、図８は、リング差δの最大値と再構成画像の相対統計雑音（ＲＭＳ値）との関係を示すグラフである。この図に示されるように、para-DRAMA(s)，para-DRAMA(t) および 3D-DRAMA の何れの場合にも、最大リング差の増大とともに、全計数値が増大する結果、統計雑音が減少する。para-DRAMA(t) の場合の統計雑音がpara-DRAMA(s) の場合の統計雑音に比して大きいのは、前者の再構成誤差が後者より小さいことに起因する。 FIG. 8 is a graph showing the relationship between the maximum value of the ring difference δ and the relative statistical noise (RMS value) of the reconstructed image. As shown in this figure, in all cases of para-DRAMA (s), para-DRAMA (t) and 3D-DRAMA, as the maximum ring difference increases, the total count value increases, resulting in statistical noise. Decrease. The statistical noise in the case of para-DRAMA (t) is larger than that in the case of para-DRAMA (s) because the former reconstruction error is smaller than the latter.

図９は、比較例における処理部の台数（ノード数）Ｎと処理時間との関係を示すグラフである。ノード数Ｎを１〜１６の各値とした。比較例では、ノード数Ｎを増加させると、処理部における計算処理に要する時間Ｔ_２は減少する傾向にあるが、その一方で、処理部の間でデータ転送に要する時間Ｔ_１は増加する傾向にあり、ノード数Ｎに依らない処理時間Ｔ_３を含め全体の処理時間（Ｔ_１＋Ｔ_２＋Ｔ_３）は減少する傾向にある。しかし、比較例では、ノード数Ｎを更に増加させても、３次元投影データから３次元再構成画像を得るための全体の所要時間の短縮化の効果は小さい。これに対して、本実施形態では、ノード数Ｎを増加させると処理部における計算処理に要する時間Ｔ_２が減少する傾向にあるだけでなく、比較例に比べて、処理部の間でデータ転送に要する時間Ｔ_１も短縮されるので、３次元投影データから３次元再構成画像を得るための全体の所要時間が短縮される。 FIG. 9 is a graph showing the relationship between the number of processing units (number of nodes) N and processing time in the comparative example. The number of nodes N was set to 1 to 16 values. In the comparative example, when the number N of nodes is increased, the time T ₂ required for the calculation process in the processing unit tends to decrease, while the time T ₁ required for data transfer between the processing units tends to increase. The total processing time (T ₁ + T ₂ + T ₃ ) including the processing time T ₃ that does not depend on the number of nodes N tends to decrease. However, in the comparative example, even if the number of nodes N is further increased, the effect of shortening the overall required time for obtaining a 3D reconstructed image from 3D projection data is small. In contrast, in the present embodiment, not only the time T ₂ required for the calculation processing tends to decrease in a processing unit increases the number of nodes N, in comparison with the comparative example, the data transfer between the processing unit Since the time T ₁ required for the process is also shortened, the overall time required for obtaining a three-dimensional reconstructed image from the three-dimensional projection data is shortened.

本実施形態に係る３次元ＰＥＴ装置１および画像処理装置２の構成を示す図である。It is a figure which shows the structure of the three-dimensional PET apparatus 1 and the image processing apparatus 2 which concern on this embodiment. 本実施形態に係る３次元ＰＥＴ装置１の検出部１０の断面図である。It is sectional drawing of the detection part 10 of the three-dimensional PET apparatus 1 which concerns on this embodiment. 本実施形態に係る画像処理装置２における画像再構成処理のフローチャートである。It is a flowchart of the image reconstruction process in the image processing apparatus 2 according to the present embodiment. 本実施形態に係る画像処理装置２における画像再構成処理のうち初期画像作成処理のフローチャートである。It is a flowchart of an initial image creation process among the image reconstruction processes in the image processing apparatus 2 according to the present embodiment. 本実施形態に係る画像処理装置２における画像再構成処理のうちサブセット群ＳＧ_δの投影データを用いた画像再構成処理のフローチャートである。It is a flowchart of an image reconstruction processing using projection data subsets group SG _[delta] of the image reconstruction processing in the image processing apparatus 2 according to the present embodiment. リング差δの最大値と再構成画像の解像度（半値全幅）との関係を示すグラフである。It is a graph which shows the relationship between the maximum value of ring difference (delta), and the resolution (full width at half maximum) of a reconstruction image. リング差δの最大値と再構成画像の再構成誤差との関係を示すグラフである。It is a graph which shows the relationship between the maximum value of ring difference (delta), and the reconstruction error of a reconstruction image. リング差δの最大値と再構成画像の相対統計雑音（ＲＭＳ値）との関係を示すグラフである。It is a graph which shows the relationship between the maximum value of ring difference (delta), and the relative statistical noise (RMS value) of a reconstruction image. 比較例における処理部の台数（ノード数）Ｎと処理時間との関係を示すグラフである。It is a graph which shows the relationship between the number of processing parts (number of nodes) N and processing time in a comparative example.

Explanation of symbols

１…３次元ＰＥＴ装置、２…画像処理装置、３…被検体、１０…検出部、２０…信号処理部、３０…収集部、４０…スイッチングハブ、５０_０〜５０_Ｎ−１…処理部、６０…制御部、７０…記憶部。
1 ... 3 dimensional PET apparatus, 2 ... image processing apparatus, 3 ... subject, 10 ... detection unit, 20 ... signal processing unit, 30 ... collecting portion, 40 ... switching _hub, 50 0 to 50 _{N-1 ...} processing unit, 60 ... control unit, 70 ... storage unit.

Claims

An image processing apparatus that performs image reconstruction by a block iteration method using projection data obtained in a three-dimensional PET apparatus in which a plurality of detector rings each including a plurality of radiation detectors are stacked in the axial direction,
An initial image creation unit for creating a three-dimensional initial image by performing image reconstruction using partial projection data including ring differences 0 and 1 in the projection data;
A plurality of processing units each configured to reconstruct a three-dimensional image by a block iteration method using a part of the projection data, starting from the initial image created by the initial image creation unit;
An image integration unit that integrates three-dimensional images obtained by image reconstruction by each of the plurality of processing units;
A control unit that controls operations of the initial image creation unit, the plurality of processing units, and the image integration unit;
With
The control unit outputs projection data of a subset group corresponding to a common ring difference among a plurality of subsets obtained by dividing the projection data by a ring difference and an azimuth to any one of the plurality of processing units. allocation,
The plurality of processing units perform image reconstruction by a block iterative method that repeats forward projection calculation, backprojection calculation, and image update calculation in parallel with each other using projection data of a subset group allocated to each other in parallel.
An image processing apparatus.

The initial image creation unit reconstructs slice images in parallel using the data obtained by dividing the partial projection data according to the position in the axial direction, and obtains an image at each position in the axial direction obtained thereby. The image processing apparatus according to claim 1, wherein the initial image is generated by synthesis.

The initial image creation unit creates a two-dimensional image at each position in the axial direction by the FBP method using the partial projection data, and synthesizes the images at the respective positions in the axial direction obtained thereby, to generate the initial image. The image processing apparatus according to claim 1, wherein:

A detector in which a plurality of detector rings each including a plurality of radiation detectors are stacked in the axial direction;
The image processing apparatus according to any one of claims 1 to 3, wherein image reconstruction is performed by a block iteration method using projection data obtained based on radiation detection by the detection unit,
A three-dimensional PET apparatus comprising: