JP2010204755A - Image processor, image reconstruction system, image processing method, and program - Google Patents

Abstract

<P>PROBLEM TO BE SOLVED: To provide an image processor for efficiently reconstructing an image and for significantly reducing calculation speed. <P>SOLUTION: In the image processor 3, an iterative operation part 32 performs arithmetic processing including first MAC processing of detection probability C<SB>ij</SB>read from a parameter storage part 39 and a voxel value x<SB>j</SB>of the j-th voxel Bx<SB>j</SB>, and a second MAC processing of the detection probability C<SB>ij</SB>and projection data on the i-th detector element. Then, the iterative operation part 32 skips execution of either the first or second MAC processing for a voxel Bx corresponding to a scintillation detector 11 whose projection data are equal to or less than a predetermined value. <P>COPYRIGHT: (C)2010,JPO&INPIT

Description

  The present invention relates to an image processing technique for reconstructing a three-dimensional image of an internal structure of a measurement object based on an amount of light emitted from the measurement object or transmitted through the measurement object.

  Radiation such as X-rays transmitted through the measurement object such as the human body or γ-rays emitted from the measurement object is detected, and the two-dimensional or three-dimensional image of the internal structure of the measurement object is reproduced based on the detection result. Research and development of image processing technology for configuration is underway. Such an image processing technique is used in an X-ray CT apparatus (X-ray computed tomography apparatus) and a positron emission tomography apparatus (PET).

  As image reconstruction algorithms, statistical image reconstruction methods such as ML-EM (Maximum Likelihood-Expectation Maximization) method and OS-EM (Ordered Subsets-Expectation Maximization) method are widely used. The OS-EM method is a kind of ML-EM method. The ML-EM method is disclosed in, for example, Non-Patent Document 1 and Non-Patent Document 2, and the OS-EM method is disclosed in, for example, Non-Patent Document 3.

  The statistical image reconstruction method is a successive approximation algorithm and has a problem that it requires a lot of calculation time. Various algorithms have been proposed as high-speed arithmetic algorithms for shortening calculation time and improving throughput. For example, Non-Patent Document 4 proposes a high-speed calculation algorithm using the geometric symmetry of the detector system.

L. A. Shepp and Y. Vardi, "Maximum likelihood reconstruction for emission tomography," IEEE Trans. Med. Imaging, 1982, Vol. MI-1, No. 2, pp. 113-122. K. Lange and R. Carson, "EM reconstruction algorithms for emission and transmission tomography," J. Comput. Assist. Tomogra., 1984, Vol. 8, No. 2, pp. 306-316. H. M. Hudson and R. S. Larkin, "Accelerated image reconstruction using ordered subsets of projection data," IEEE Trans. Med. Imaging, 1994, Vol. 13, pp. 601-609. A. Motta et al., "Fast 3D-EM reconstruction using Planograms for stationary planar positron emission mammography camera," Computerized Medical Imaging and Graphics, 2005, Vol. 29, pp. 587-596.

  However, even if a conventional high-speed calculation algorithm is used, since the calculation amount is enormous, there is a problem that it is difficult to say that the calculation speed is clinically practical. For example, if the calculation time until the image is reconstructed is long, the time required for diagnosis becomes long, so it is difficult to say that it is clinically practical.

  The present invention has been made in view of the above problems, and provides an image processing apparatus, an image reconstruction system, an image processing method, and a program that can efficiently execute image reconstruction processing and can greatly improve the calculation speed. To do.

According to the present invention, projection data indicating the amount of photons flying from the object to be measured to each of the plurality of detector elements, and a plurality of voxels arranged in a three-dimensional manner within the detection region where the object to be measured is arranged And a third order of the internal structure of the object to be measured by repeatedly executing a calculation process for estimating a voxel value for each voxel from the projection data using a statistical image reconstruction method based on the analyzed model. An iterative operation unit for reconstructing the original image;
a parameter storage unit storing data indicating a detection probability C _{ij at} which the photon that has passed through the j th voxel is detected by the i th detector element;
An image processing apparatus is provided.
In this image processing apparatus, the iterative operation unit is
(A) a first product-sum operation process of the detection probability C _ij read from the parameter storage unit and the voxel value of the j-th voxel;
A second product-sum operation process of the detection probability C _ij and the projection data applied to the i-th detector element;
And performing the arithmetic processing including
(B) The execution of at least one of the first or second product-sum operation processing is skipped for the voxel corresponding to the detector element whose projection data is equal to or less than a predetermined value.

  In the image processing apparatus of the present invention, as a more specific aspect, the iterative operation unit sets a voxel value of the voxel corresponding to a detector element whose projection data is equal to or less than the predetermined value to zero. May be.

In the image processing apparatus of the present invention, as a more specific aspect, a zero region that is a set of voxels corresponding to detector elements whose projection data is equal to or less than the predetermined value, and the projection data is the predetermined value. A two-dimensional image generation unit for generating a two-dimensional image of voxel values, including a non-zero region that is a set of voxels corresponding to the detector elements larger than
Generating a three-dimensional initial image including the two-dimensional image, and determining a voxel corresponding to the zero region in the three-dimensional initial image;
Further comprising
The iterative operation unit may skip execution of at least one of the first or second product-sum operation processing for voxels corresponding to the zero region in the three-dimensional initial image.

Further, in the image processing apparatus of the present invention, as a more specific aspect, the arithmetic processing is repeatedly executed by reducing the size of the voxel in a plurality of stages.
For the voxels after miniaturization included in the voxels corresponding to the zero region in the three-dimensional initial image before miniaturization, the iterative operation unit includes at least one of the first or second product-sum operation processing Execution may be skipped.

  According to the present invention, there is provided the image processing apparatus as described above including a pair of scintillation detectors arranged opposite to each other as the detector element, and the scintillation detector is arranged in a predetermined direction from the object to be measured. An image reconstruction system is provided that detects photon pairs emitted on both sides.

Further, according to the present invention, based on the projection data indicating the amount of photons flying from the object to be measured to each of the plurality of detector elements, the internal structure of the object to be measured using a statistical image reconstruction method. An image processing method for reconstructing a three-dimensional image is provided.
The image processing method includes (a) a model generation step of generating an analysis model by three-dimensionally arranging a plurality of voxels within a detection region where the object to be measured is arranged;
(B) A forward projection step of performing a first product-sum operation process of a detection probability C _{ij in} which the photon that has passed through the j-th voxel is detected by the i-th detector element and a voxel value of the j-th voxel When,
(C) a back projection step of performing a second product-sum operation process of the detection probability C _ij and the projection data applied to the i-th detector element;
(D) an iterative calculation step of repeatedly executing the forward projection step and the backward projection step to estimate a voxel value for each voxel;
Including
For the voxel corresponding to the detector element whose projection data is equal to or less than a predetermined value, the execution of at least one of the first or second product-sum operation processing is skipped.

Further, according to the present invention, based on the projection data indicating the amount of photons flying from the object to be measured to each of the plurality of detector elements, the internal structure of the object to be measured using a statistical image reconstruction method. A program for causing an information processing apparatus to execute an image reconstruction process for reconstructing a three-dimensional image is provided.
In this program, the image reconstruction process is
(A) a model generation process for generating an analysis model by three-dimensionally arranging a plurality of voxels within a detection region in which the object to be measured is arranged;
(B) A forward projection process in which a first product-sum operation process is performed between the detection probability C _{ij at} which the photon that has passed through the j-th voxel is detected by the i-th detector element and the voxel value of the j-th voxel. When,
(C) a back projection process for performing a second product-sum operation process on the detection probability C _ij and the projection data applied to the i-th detector element;
(D) an iterative calculation process of repeatedly executing the forward projection process and the backward projection process to estimate a voxel value for each voxel;
Including
For the voxel corresponding to the detector element whose projection data is equal to or less than a predetermined value, the execution of at least one of the first or second product-sum operation processing is skipped.

  In the above invention, the fact that photons fly from the object to be measured to the detector element means that the photon is detected by the interaction between the photon and the detector element. Also, when a photon comes from the object to be measured, a photon emitted from the back of the object to be measured passes through the object to be measured and comes to the detector element, and a photon is emitted using the object to be measured as a light source. Including cases. Hereinafter, when expressing that a photon is emitted from the object to be measured, it means to include both of the above aspects.

The various components of the present invention need only be formed so as to realize their functions. For example, dedicated hardware that exhibits a predetermined function, image processing in which a predetermined function is provided by a computer program It can be realized as an apparatus, a predetermined function realized in the image processing apparatus by a computer program, an arbitrary combination thereof, or the like.
Further, the various components of the present invention do not have to be individually independent, that a plurality of components are formed as one member, and one component is formed of a plurality of members. That a certain component is a part of another component, a part of a certain component overlaps a part of another component, and the like.

In the image processing method of the present invention, a plurality of steps are described in order, but the description order does not limit the order in which the plurality of steps are executed. For this reason, when carrying out the image processing method of the present invention, the order of the plurality of steps can be changed within a range that does not hinder the contents.
Furthermore, the image processing method of the present invention is not limited to being executed at a timing at which a plurality of steps are individually different. For this reason, another process may occur during execution of a certain process, or a part or all of the execution timing of a certain process and the execution timing of another process may overlap.

Further, the repetitive arithmetic processing referred to in the present invention is a CPU (Central Processing Unit), a ROM (Read Only Memory), a RAM (Random Access Memory), an I, so that a computer program can be read and the corresponding data processing can be executed. It can be implemented using hardware constructed by general-purpose devices such as / F (Interface) units, dedicated logic circuits constructed to execute predetermined data processing, combinations thereof, and the like.
Note that causing the information processing apparatus to execute various operations corresponding to the program in the present invention includes causing the information processing apparatus to control operations of various devices.

Furthermore, in the present invention, storing or storing data means that the apparatus of the present invention has a function of holding the data at least when executing arithmetic processing. Yes. Therefore, for example, projection data indicating the light amount and data indicating the detection probability C _ij are already registered when the apparatus of the present invention is shipped, and are not registered when shipped. However, it is allowed to be registered at the time of the arithmetic processing or during the arithmetic processing.

  According to the image processing technique of the present invention, the execution of at least one of the first and second sum-of-products operations is skipped for voxels corresponding to detector elements whose projection data indicating the amount of light of incoming photons is equal to or less than a predetermined value. To do. This is because the voxel value can be regarded as a predetermined value (for example, zero) uniformly for the voxels existing in the region where the photons do not fly. The first and second product-sum operation processes are processes with a large amount of calculation in an iterative operation process using a statistical image reconstruction method. Therefore, with respect to the above voxels, the total time of the arithmetic processing can be greatly shortened by skipping the product-sum arithmetic processing.

It is a functional block diagram which shows schematic structure of the image reconstruction system of one Embodiment which concerns on this invention. It is a schematic diagram for demonstrating measured projection data. It is a flowchart which shows roughly the process sequence by ML-EM method. It is a flowchart which shows the specific procedure of the production process of a two-dimensional image. It is a flowchart which shows the specific procedure of the production output of a three-dimensional initial image. It is a flowchart which shows the specific procedure of an image reconstruction process. (A) is an analysis model and a detection area, (b) is an analysis model and an object to be measured, and (c) is a perspective view showing an example of a two-dimensional image. It is a perspective view which shows a three-dimensional initial image. (A) is a schematic model used for a rough reconstruction process, (b) is a perspective view of the fine model used for a fine reconstruction process.

Embodiments according to the present invention will be described below with reference to the drawings. In all the drawings, the same components are denoted by the same reference numerals, and detailed description thereof is appropriately omitted so as not to overlap.
FIG. 1 is a functional block diagram showing a schematic configuration of an image reconstruction system 1 and an image processing apparatus 3 according to an embodiment of the present invention.

The image reconstruction system 1 of this embodiment includes an image processing device 3 and a detection unit 2.
The image processing device 3 includes an iterative calculation unit 32 and a parameter storage unit 39.
The iterative calculation unit 32 includes projection data indicating the amount of photons flying from the DUT 20 to each of a plurality of detector elements (crystal pair 14: see FIG. 2), and a detection region in which the DUT 20 is arranged. Based on an analysis model in which a plurality of voxels Bx (see FIG. 7) are arranged three-dimensionally in 22, a voxel value x for each voxel Bx is estimated from projection data using a statistical image reconstruction method. The three-dimensional image of the internal structure of the DUT 20 is reconstructed by repeatedly executing the arithmetic processing.
The parameter storage unit 39 stores data indicating a detection probability C _{ij in} which a photon that has passed through the j th voxel Bx _j is detected by the i th scintillation detector 11 i.

Then, the iterative operation unit 32 of the present embodiment includes a first product-sum operation process of the detection probability C _ij read from the parameter storage unit 39 and the voxel value x _j of the j-th voxel Bx _j , the detection probability A calculation process including a second product-sum calculation process of C _ij and the projection data applied to the i-th detector element is performed.
In addition, the iterative operation unit 32 of the present embodiment skips at least one of the first and second product-sum operation processes for the voxel Bx corresponding to the scintillation detector 11 whose projection data is equal to or less than a predetermined value.

  The image reconstruction system 1 constitutes an opposed PET apparatus having a pair of detection plates 10 (10A, 10B) facing each other, measures the dose of radiation emitted from the measured object 20, and measures the measured object 20 A three-dimensional image of the internal structure is estimated (reconstructed) together with the outer shape of the.

  The detection unit 2 includes a pair of detection plates 10A and 10B having detection surfaces 12 (12A and 12B) facing each other, and a drive mechanism 41 that translates or rotates the detection plates 10A and 10B. . A measurement object 20 such as a human body is sandwiched between the detection plates 10A and 10B.

On the detection plate 10, a large number of scintillation detectors 11 (11A, 11B) are arranged in a two-dimensional plane.
The arrangement of the scintillation detectors 11 is not particularly limited, and may be a grid, a staggered pattern, or any other pattern. In the present embodiment, the detection plate 10 in which the scintillation detectors 11 are arranged in a grid (matrix) along the XY plane is illustrated.
The number of scintillation detectors 11 arranged in each direction is arbitrary, but in this embodiment, m scintillation detectors 11 are arranged along the X and Y axis directions to form the detection plate 10. It shall be.

The scintillation detector 11 includes a scintillator crystal (not shown) and a light receiving element (not shown) that receives fluorescence generated by the disappearance of radiation inside the scintillator crystal.
The opposing detection surfaces 12 of the detection plates 10A and 10B are constituted by a set of end faces of scintillator crystals. Each scintillator crystal is arranged flush with each other, and the detection surface 12 is formed flat.
The scintillator crystal of the present embodiment is a quadrangular prism whose end surface constituting the detection surface 12 forms a square. However, the shape and dimensions of the end face of the scintillator crystal are not limited to this.
And this embodiment detects the radiation of radiation from the to-be-measured object 20 for every combination of the scintillation detectors 11A and 11B which oppose.

  That is, the image reconstruction system 1 according to the present embodiment includes a pair of scintillation detectors 11 (crystal pair 14) arranged to face each other as a detector element. The scintillation detector 11 detects a pair of photons emitted from the device under test 20 to both sides in a predetermined direction.

The drive mechanism 41 has a function of driving the pair of detection plates 10 </ b> A and 10 </ b> B in accordance with a control signal from the control unit 40 incorporated in the image processing apparatus 3. In the present embodiment, the drive mechanism 41 narrows the interval between the detection plates 10A and 10B, rotates the pair of detection plates 10A and 10B around the X axis, or moves the pair of detection plates 10A and 10B to the X axis. Or move in the direction.
Therefore, in the image reconstruction system 1 of the present embodiment, the distance between the pair of detection surfaces 12A and 12B is variable by driving the drive mechanism 41. Accordingly, the object to be measured 20 sandwiched between the detection plates 10A and 10B is pressed by the detection surfaces 12A and 12B and widened. In the image reconstruction system 1, the radiation flying from the pressed measurement object 20 is received.
The user can cause the drive mechanism 41 to drive the detection plates 10A and 10B by giving an instruction to the control unit 40 through the input operation unit 51 such as a keyboard or a pointing device.

  Radioactive isotopes (RI) that generate positrons are distributed inside the object to be measured 20. For example, if the object to be measured 20 is an object such as a human body or an animal, the radioisotope can be distributed in the object to be measured 20 by administering a radioisotope to the object. Positrons are generated in the decay process of radioisotopes. When a positron interacts with an electron in the object to be measured 20 and disappears, two gamma rays (photon pairs) having an energy of about 511 keV are generated, and these gamma rays are emitted in opposite directions. Therefore, the pair of scintillation detectors 11A and 11B having point symmetry with respect to the point where the positron has disappeared respectively detect the two gamma rays almost simultaneously. Therefore, the point where the positron has disappeared exists on the line connecting the pair of scintillation detectors 11A and 11B that respectively detected the photon pair.

The scintillator crystal used in the scintillation detectors 11A and 11B is a crystal body that converts photons of high energy radiation such as X-rays and gamma rays into a large number of photons (fluorescence) having lower energy. As this type of crystal, for example, a known crystal such as BGO or LSO may be used.
The light receiving elements included in the detection plates 10A and 10B amplify the fluorescence output generated by the scintillation detectors 11A and 11B and convert the fluorescence output into an electrical signal. The electric signal is transferred to the projection data generation unit 31 of the image processing apparatus 3 together with the address information of the scintillation detectors 11A and 11B that have detected the photons via a transfer means (not shown) such as a cable.

  The projection data generation unit 31 counts the photons detected for each crystal pair 14 based on the electrical signals transferred from the detection plates 10A and 10B, and for each crystal pair 14 (detector element) based on the counting result. Projection data (actual measurement projection data) is generated. The actually measured projection data is stored in the data storage unit 38 in the memory 37.

FIG. 2 is a schematic diagram for explaining the actually measured projection data y. Photons emitted from the positron annihilation point (ray source) 21 on the projection line in the object to be measured 20 to both sides in a predetermined direction are detected almost simultaneously by the scintillation detectors 11A and 11B of the detection plates 10A and 10B.
The radiation source 21 is distributed over almost the entire object to be measured 20, and photons are emitted from almost the entire object to be measured 20 in various directions. Since the radiation source 21 is concentrated at high density in a highly active site such as a cancer cell, many photons are detected in the crystal pair 14 sandwiching the cancer cell.
As shown in FIG. 2, when focusing on a specific scintillation detector 11B, among the scintillation detectors 11A constituting the detection plate 10A, the scintillation detector 11A1 that sandwiches the radiation source 21 and the scintillation that does not sandwich the radiation source 21 There may be a detector 11A2.
The same applies to the other scintillation detectors 11B.

When the projection data generation unit 31 (see FIG. 1) detects that the photons are detected within a predetermined time difference by the opposing scintillation detectors 11A and 11B, the crystal pair corresponding to the combination of the scintillation detectors 11A and 11B. For 14 _i (i is the number of the crystal pair), the light quantity of the photon is integrated and counted. By performing such measurement over a predetermined time, the actual projection data y (actual projection data y _i ) for each crystal pair number can be obtained.
The distribution density of the radiation source 21 is reconstructed by using a statistical image reconstruction method based on the accumulated distribution of the actually measured projection data y _i , and the presence and position of the cancer cell 23 and its position are estimated. .

The iterative operation unit 32 shown in FIG. 1 includes a plurality of voxels Bx (a small three-dimensional element such as a cubic shape or a rectangular shape, see FIG. 7) in a three-dimensional manner within the detection region 22 where the device under test 20 is arranged. Image reconstruction processing is executed using the arranged analysis model 25.
The iterative calculation unit 32 reads the measured projection data y _i from the data storage unit 38, and uses the measured projection data y _i to calculate the voxel value x corresponding to the radioactivity concentration for each voxel Bx according to the statistical image reconstruction method. The estimation processing is repeatedly executed.
This has a function of reconstructing a tomographic image or a three-dimensional image of the internal structure of the DUT 20. In the present embodiment, the three-dimensional image reconstruction includes outputting a tomographic image that is a two-dimensional image based on the three-dimensional image of the voxel value.

  In the present embodiment, a typical ML-EM method can be used as the statistical image reconstruction method. However, the statistical image reconstruction method that can be used is not limited to this.

  The iterative calculation unit 32 includes functional blocks of a front projection calculation unit 33, a rear projection calculation unit 34, and a density calculation unit 35. These functional blocks repeatedly execute processing according to the ML-EM method. In the ML-EM method, the following estimation formula (1) is used.

However, C _j in the estimation formula (1) is given by the following formula (1A).

In the estimation formula (1), k is the number of iterations, and x _j ^(k) is the radioactivity concentration (voxel value x) of the j-th voxel Bx _j in the k-th process. y _i is a value (projection value) of the above-mentioned measured projection data obtained based on the detection result by the i-th crystal pair (crystal pair 14 _i ).
In the case of the image reconstruction system 1 of the present embodiment, the detection plate 10 is configured by arranging m scintillation detectors 11 along the X and Y axis directions. Therefore, each detection plate 10 has m ² scintillation detectors 11. Then, there are only m ⁴ crystal pairs 14 that are combinations of the scintillation detectors 11, and the maximum value of i is m ⁴ .

X _j ^{(k + 1)} is the radioactivity concentration in the j-th voxel Bx _j in the (k + 1) -th processing. Therefore, the k-th voxel value _x ^{j (k),} is (k + 1) -th voxel values _x ^{j (k + 1)} is calculated.
Here, the initial value x _j ⁽⁰⁾ of x _j ^(k) is set to any value as long as it is a positive real number. By doing so, a three-dimensional image with a predetermined accuracy can be reconstructed.

In the estimation formula (1), C _ij is a probability (detection probability) that a photon pair flying from the j-th voxel Bx _j is detected by the i-th crystal pair 14 _i . It is known that the number of photon pairs flying from the voxel Bx _j and incident on the crystal pair 14 _i depends on the detection probability C _ij . The values of the detection probabilities C _ij and C _j are calculated in advance and stored in the parameter storage unit 39.

The detection probability C _ij can be regarded as an element in the i-th row and j-th column of the matrix (referred to as “system matrix”). The value of the detection probability C _ij is susceptible to incorporate physical factors such as scattering and attenuation, in the present embodiment, the value of the detection probability C _ij is, j th for the i-th crystal pairs 14 _i It can be determined only by the geometric position of the voxel Bx _j .
As the value of the detection probability C _ij, a value determined by a ray-tracing method, a Monte Carlo simulation method, or the like can be set.

  Next, the image processing method of the present embodiment (hereinafter referred to as the present method) realized as the operation of the image processing apparatus 3 will be described.

First, an outline of this method will be described.
This method uses the statistical image reconstruction method based on projection data indicating the amount of photons flying from the device under test 20 to each of a plurality of detector elements (crystal pairs 14). A method for reconstructing a three-dimensional image of a structure.
This method is
(A) a model generation step for generating an analysis model by three-dimensionally arranging a plurality of voxels Bx in the detection region 22 where the device under test 20 is arranged;
(B) First product-sum operation processing of a detection probability C _{ij in} which a photon that has passed through the j-th voxel Bx _j is detected by the i-th crystal pair 14 _i and the voxel value x _j of the j-th voxel Bx _j Performing a forward projection step,
(C) a back projection step of performing a second product-sum operation process of the detection probability C _ij and the projection data _concerning the i-th crystal pair 14 _i ;
(D) an iterative operation step of repeatedly executing the forward projection step and the backward projection step to estimate a voxel value x for each voxel Bx;
including.
Then, this method skips execution of at least one of the first and second product-sum operation processes for the voxel Bx corresponding to the crystal pair 14 whose projection data is equal to or less than a predetermined value.

  The method further includes (e) an actual measurement step of measuring the amount of photons in a state where the device under test 20 is pressed and widened.

Next, the method will be described in more detail with reference to the drawings.
FIG. 3 is a flowchart schematically showing a processing procedure by the ML-EM method.
FIG. 4 is a flowchart showing a specific procedure of the two-dimensional image creation process (step S12).
FIG. 5 is a flowchart showing a specific procedure for generating and outputting a three-dimensional initial image (step S13).
FIG. 6 is a flowchart showing a specific procedure of the image reconstruction process (steps S14 and S15).

In step S _{< b} > 10, the iterative operation unit 32 reads the detection probability {C _ij } and the coefficient {C _j } from the parameter storage unit 39.
Subsequently, the iterative operation unit 32 reads the actual projection data {y _i } from the data storage unit 38 (step S11). In this embodiment, the detection probability {C _ij }, the coefficient {C _j }, and the actually measured projection data {y _i } are first read out at a time. However, the present invention is not limited to this, and is read out when necessary. May be.

Subsequently, the iterative operation unit 32 initializes the value of the voxel value {x _j } corresponding to the radioactivity concentration, based on the two-dimensional image created based on the actually measured projection data y _i and the two-dimensional image. A three-dimensional initial image created in this way is created (steps S13 and S14).
Thereafter, the iterative operation unit 32 executes an image reconstruction process (step S15) using a rough model with a large voxel size and an image reconstruction process (step S16) using a fine model with a finer voxel size. .
In this method, an example in which the image reconstruction processing is performed by reducing the voxel size in two stages will be described as an example. However, the stage of miniaturization is not limited to this and may be three or more stages.
Further, miniaturization of the voxel size is optional in the present invention, and image reconstruction processing may be performed using only a single analysis model.

  First, specific means for the image reconstruction process (steps S15 and S16) will be described with reference to FIG. The three-dimensional initial image created in this embodiment as an initial image that is preferably used for such image reconstruction processing will be described later.

Referring to FIG. 6, first, the iterative operation unit 32 sets the upper limit of the number of iterations to an appropriate value (step S140).
The number of iterations k1 set in the image reconstruction process based on the schematic model (hereinafter sometimes referred to as “schematic reconstruction process”) and the image reconstruction process based on the fine model (hereinafter referred to as “fine reconstruction process”). The number of iterations k2 set in step 2) may be the same as or different from each other.
In addition, the values of the number of iterations k1 and k2, and the magnitude relationship between them, are arbitrary, but from the viewpoint of reducing the total time required for the arithmetic processing and performing image reconstruction with the same level of accuracy, k1 is larger than k2. It is good to set.

The forward projection calculation unit 33 in FIG. 1 performs a product-sum operation (first product-sum operation) between the detection probability C _ij and the voxel value x _j ^(k) corresponding to the estimated density according to the following equation (2). Then, the front projection value p _i is calculated (step S141).

A large part of the calculation amount in the estimation formula (1) is the calculation of the forward projection data value (forward projection value) p _i given by the formula (2).
The forward projection calculation unit 33 may decompose the first product-sum operation of the above equation (2) into a plurality of threads (groups) and execute it by parallel operation.

After the process of step S141 is completed, the back projection calculation unit 34 in FIG. 1 calculates a product q _i of the actually measured projection data y _i and the reciprocal of the front projection value p _i according to the following equation (3): A second product-sum operation (second product-sum operation) between the detection probability C _ij and the product q _i is executed to calculate a rear projection value D _j (step S142).

The calculation of the rear projection value D _j is also a part with a large amount of calculation, like the calculation of the front projection value.

Then, the density calculation unit 35 calculates the voxel value x _j ^{(k + 1)} (= D _j × x _j ^(k) ) corresponding to the (k + 1) -th estimated density by using the rear projection value D _j (step) S143). Thereafter, it is determined whether or not the number of iterations k (k1, k2) has reached the upper limit (step S144), and the procedures of steps S141 to S144 are repeatedly executed until the number of iterations k reaches the upper limit. When the number of iterations k finally reaches the upper limit (YES in step S144), the process returns to the flowchart of FIG. 3 to end the rough reconstruction process (step S14) or the fine reconstruction process (step S15).

And the density | concentration calculation part 35 outputs voxel value { _xj } (step S16).

The display control unit 36 generates a three-dimensional image from the voxel values {x _j } output from the density calculation unit 35. The generated image can be displayed on the display 50 as a two-dimensional tomographic image or a three-dimensional image. Alternatively, the voxel value {x _j } output from the density calculation unit 35 may be stored in the memory 37.

<Three-dimensional initial image creation processing>
Here, in the conventional statistical image reconstruction method, a positive value is set for all voxels as the initial value of the voxel value, and the iterative operation of the voxel value x _j ^(k) is performed for all the voxels Bx _j . I was doing it.
However, as described above, the voxel value x _j ^{(k + 1)} corresponding to the ^{(k + 1)} -th estimated density is obtained by multiplying the voxel value x _j ^(k) corresponding to the k-th estimated density by the rear projection value D _j. Calculated. Therefore, once a voxel in which x _j ^(k) has become zero, even if the subsequent iterative operation is performed, the voxel value is zero.

Therefore, in this method, iterative computation unit 32 sets the voxel value x _j of the voxels Bx which measured projection data y _i corresponds to the crystal pair 14 is less than the predetermined value to zero. For the voxel Bx _j in which the voxel value x _j is set to zero, the first product-sum operation process (the above equation (2)) or the rear projection step by the forward projection calculation unit 33 in the forward projection step (step S141). The execution of at least one of the second product-sum operation processing (the above expression (3)) by the back projection calculation unit 34 in (Step S142) is skipped.
As a result, the total time required for the calculation processing by the iterative calculation unit 32 is reduced.

As shown in FIG. 1, the image processing apparatus 3 according to the present embodiment has an initial value x _j ⁽ of a voxel value x _j for a voxel Bx corresponding to a crystal pair 14 whose measured projection data y _i is equal to or less than a predetermined value. ⁰⁾ has an initial value setting unit 60 that gives zero.

The initial value setting unit 60 includes a two-dimensional image generation unit 61 and a three-dimensional initial image generation unit 62.
The function of the initial value setting unit 60 will be described with reference to FIG.

As the model generation step (a), the initial value setting unit 60 creates an analysis model 25 in which a plurality of voxels Bx are three-dimensionally arranged inside the detection region 22 sandwiched between the detection plates 10A and 10B. In creating the three-dimensional initial image generation unit 62, in this method, the cell size of the voxel Bx and the size of the scintillation detector 11 are matched.
In this method, the entire area sandwiched between the detection plates 10A and 10B is set as the detection area 22, and the entire detection area 22 is divided by the voxel Bx to constitute the analysis model 25. It is not limited to this. A part of the region sandwiched between the detection plates 10 </ b> A and 10 </ b> B may be the detection region 22, and the analysis model 25 may be generated corresponding to the partial region of the detection region 22.

FIG. 7A is a perspective view of the analysis model 25 and the detection region 22. As shown in the figure, in this embodiment, the shape of the voxel Bx is a cube. The detection plates 10A and 10B are installed in the XY plane, and the arrangement directions of the scintillation detectors 11A and 11B are the X-axis and Y-axis directions. That is, the detection plates 10A and 10B face each other in the Z-axis direction.
In this method, it is assumed that m voxels Bx are arranged in the X, Y, and Z axis directions with respect to the detection region 22.

FIG. 7B is a perspective view of the analysis model 25 and the DUT 20. As shown in the figure, the DUT 20 is set inside the detection region 22. In the figure, the detection plates 10A and 10B are not shown.
The measurement object 20 is exemplified by a subject's breast. The breast is a site characterized by being substantially rotationally symmetric. Specifically, it forms a hemisphere or a spheroid (half body).
However, in the image reconstruction system 1 of the present embodiment, a human body such as a region including a subject's lymph node, a region including a visceral organ such as a stomach or a liver, or another object can be used as the DUT 20.
The DUT 20 is sandwiched between the detection plates 10A and 10B with the rotation axis AX aligned with the X-axis or Y-axis direction. In the figure, the rotation axis AX is made to coincide with the X-axis direction.

Here, as described above with reference to FIG. 2, the projection data generation unit 31 measures the actual projection data y _i by integrating the quantity of photons emitted from the measurement object 20 for each crystal pair 14 _i (Measurement step (e)). Such measurement is performed in a state where the drive mechanism 41 is driven to adjust the facing distance between the detection plates 10A and 10B and the object to be measured 20 is pressed and widened by the detection plates 10A and 10B.
As a result, photons are detected in an area wider than the Z-axis direction projection of the DUT 20 in the natural state. On the other hand, the cancer cell 23 exists inside the object to be measured 20. For this reason, the measured projection data y _i is measured in a state where the area around the cancer cell 23 in the measured object 20 is expanded by pressing the measured object 20. Therefore, the non-zero region 26NZ of the three-dimensional initial image 26 to be described later sufficiently includes the cancer cell 23, and the region to be measured 20 in which the radiation source 21 (see FIG. 2) may be included. It never becomes set erroneously zero initial value of the voxel values x _j for.

The two-dimensional image generation unit 61 includes a zero region 24ZA that is a set of voxels Bx corresponding to the crystal pair 14 _{i in} which the actual projection data y _i is equal to or less than a predetermined value, and a crystal pair in which the actual projection data y _i is greater than the predetermined value. The two-dimensional image 24 of the voxel value x _j is generated including the non-zero region 24NZ that is a set of voxels Bx corresponding to 14 _i .

  With reference to FIG. 4, the two-dimensional image creation processing (step S12 in FIG. 3) by the two-dimensional image generation unit 61 will be specifically described.

The two-dimensional image generation unit 61 designates the crystal pairs 14 in order (step S120). The detection plate 10 of the present embodiment includes m ² scintillation detectors 11 each, and there are m ⁴ crystal pairs 14.
Next, the two-dimensional image generation unit 61 reads out the actually measured projection data y _i related to the designated crystal pair 14 _i from the data storage unit 38 (step S121).

The two-dimensional image generation unit 61 performs zero determination on the actually measured projection data y _i (step S122). If zero (step S122: No), the two-dimensional image generation unit 61 performs processing on the next crystal pair 14 _i .
In the zero determination of the measured projection data y _i , a predetermined minute threshold is set, and when the measured projection data y _i is equal to or less than the threshold, the measured projection data y _i is close to zero so that the influence on the reconstructed image can be ignored. The measured projection data y _i may be determined to be zero as being (substantially zero). Hereinafter, the value of the measured projection data y _i or the voxel value x _j being zero includes the case where the value is zero and the case where the value is substantially zero.

Then, when the actually measured projection data y _i is non-zero (step S122: Yes), it is understood that the radiation source 21 exists between the scintillation detectors 11A and 11B constituting the crystal pair 14 _i . The two-dimensional image generation unit 61 adds the value of the actually measured projection data y _i to the voxel Bx _j sandwiched between the crystal pairs 14 _i that is present at the intermediate position in the Z-axis direction (step S123). .
In the detection plate 10 of this embodiment, since the scintillation detector 11 is arranged flat in a plane, the intermediate position of the crystal pair 14 _i is the Z height of the detection plate 10A and the Z height of the detection plate 10B. A two-dimensional plane 27 existing in the middle is drawn.
Thereby, a non-zero voxel value x _j is set for the voxel Bx _j constituting the two-dimensional plane 27.

The two-dimensional image generation unit 61 repeats the above processing until i reaches m ⁴ (step S124).
As a result, a large voxel value x _j is set for a voxel Bx _j arranged at an intermediate height between the scintillation detectors 11A and 11B and having a high detection probability in the crystal pair 14.

Then, the two-dimensional image generation unit 61 performs an appropriate operation such as dividing the accumulated voxel value x _j by the coefficient {C _j } to calculate the density value of the two-dimensional image 24 (step S125).

FIG. 7C is a perspective view showing an example of the two-dimensional image 24 generated by the two-dimensional creation process.
The two-dimensional image 24 includes a zero area 24ZA where the density value is zero and a non-zero area 24NZ where the density value is non-zero.
The plurality of voxel values x _j in the non-zero region 24NZ of the two-dimensional image 24 are different from each other according to the value of the measured projection data y _i . That is, in this method, the value of the measured projection data y _i is added to the voxel Bx on the two-dimensional plane 27 in step S123, so that the magnitude of the measured projection data y _i is the pixel density of the two-dimensional image 24. It is reflected as a value.
In the non-zero region 24NZ of the two-dimensional image 24, a region with a high pixel density value (high concentration region 28) projects cancer cells 23 (see FIG. 2) in which the radiation sources 21 are densely projected onto the two-dimensional plane 27. It corresponds to the shape.

  Next, the three-dimensional initial image generation unit 62 creates a three-dimensional initial image 26 that includes the two-dimensional image 24, and determines a voxel corresponding to the zero region 26ZA in the three-dimensional initial image 26.

  The shape of the three-dimensional initial image 26 is not particularly limited. In the present method, as a first example, a case will be described in which the three-dimensional initial image generation unit 62 creates a three-dimensional initial image 26 obtained by rotating the two-dimensional image 24 about its axis. FIG. 8A is a perspective view showing a three-dimensional initial image 26 created by this method.

The three-dimensional initial image 26 is obtained by rotating the voxel Bx of the two-dimensional image 24 formed by the non-zero region 24NZ and the zero region 24ZA on the two-dimensional plane 27 by 180 degrees or 360 degrees around the rotation axis AX.
Of sterically present voxel Bx in the detection area 22, the one corresponding to the rotational position of the voxels constituting the zero region 24ZA in the two-dimensional plane 27 are both set to voxel values x _j to zero.
Then, for those corresponding to the rotational position of the voxels constituting the non-zero regions 24NZ in the two-dimensional plane 27 are both set to voxel values x _j to a non-zero value. The non-zero set value is not limited, but by reflecting the density of the pixel density value in the two-dimensional image 24 also in the three-dimensional initial image 26, image reconstruction calculation based on the statistical image reconstruction method is performed. Convergence can be accelerated.

  With reference to FIG. 5, the three-dimensional initial image creation processing (step S13 in FIG. 3) by the three-dimensional initial image generation unit 62 will be specifically described.

The three-dimensional initial image generation unit 62 calculates, for each pixel of the two-dimensional image 24, for each pixel in the X-axis direction and the Y-axis direction, a sum of integrated values of the pixel coordinate value and the pixel density value (step S130). .
Next, the three-dimensional initial image generation unit 62 divides the sum by the total amount of the measured projection data y _{i in} the X-axis direction or the Y-axis direction, so that the barycentric coordinates (Gx, Gy) is obtained (step S131).
Then, the three-dimensional initial image generation unit 62 sets a straight line that passes through the barycentric coordinates (Gx, Gy) and is parallel to the X axis as the rotation axis AX of the two-dimensional image 24. Specifically, a straight line with y = Gy and z = m / 2 becomes the rotation axis AX.

The three-dimensional initial image generation unit 62 sequentially specifies the voxels Bx, and calculates the distance between the voxel Bx _j and the rotation axis AX in the YZ plane (step S132).
In this method, since m voxels Bx are arranged in the X, Y, and Z axis directions with respect to the detection region 22, there are m ³ voxels Bx in total. Therefore, for example, voxels Bx are specified in order by arguments from j = 1 to m ³ .

Then, the three-dimensional initial image generation unit 62 substitutes and gives the voxel value x of the voxel Bx on the two-dimensional image 24 having the same distance from the rotation axis AX as the voxel value x _j of the voxel Bx _j ( Step S133). Thus, zero is given as an initial value x _j ⁽⁰⁾ to the voxel Bx existing at the rotation position of the zero region 24ZA in the two-dimensional image 24. A non-zero value is given as an initial value x _j ⁽⁰⁾ to the voxel Bx existing at the rotational position of the non-zero region 24NZ in the two-dimensional image 24.
With this processing, the pixel density value in the two-dimensional image 24 is given around the rotation axis AX, and the three-dimensional initial image 26 is created.
Note that there may be two voxels Bx on the two-dimensional plane 27 having the same distance from the rotation axis AX in the ± Y-axis direction with respect to an arbitrary voxel Bx _j in the analysis model 25. In such a case, with respect to the voxel Bx _j, may be given a voxel value x _j of either voxels Bx on a two dimensional plane 27.

Here, since the rotating solid formed by rotating the two-dimensional image 24 about the rotation axis AX has a cylindrical shape, the voxel Bx having a large absolute value of the Y-axis and Z-axis components in the cubic analysis model 25 is the cylinder. It will be located outside the shape. Regarding this voxel Bx, assuming that it corresponds to the external region of the DUT 20, the present method gives zero to the initial value x _j ⁽⁰⁾ .

The three-dimensional initial image generation unit 62 sets the initial value x _j ⁽⁰⁾ of the voxel value x _j for all the voxels Bx _j (step S134).

  Then, returning to step S14 in FIG. 3, the iterative operation unit 32 performs a forward projection step (see step S141: see FIG. 6) or a backward projection step (step S142) for the voxel corresponding to the zero region 26ZA in the three-dimensional initial image 26. ) At least one of the first or second product-sum operation processing is skipped.

In this method, the boundary between the zero region 26ZA and the non-zero region 26NZ of the three-dimensional initial image 26 is set to coincide with the outer surface of the object 20 to be measured, or slightly outside the outer surface of the object 20 to be measured. It is preferable. Thereby, the product-sum calculation process is repeatedly performed on the voxels including the outer surface of the object to be measured 20, and the product-sum operation process on the external area of the object to be measured 20 is skipped. You can balance and improve time.
Here, when the DUT 20 has an axisymmetric shape, the initial value x _j ⁽⁰⁾ of each voxel Bx _j of the three-dimensional initial image 26 is non-zero with respect to the minimum region including the DUT 20. Setting the value is preferable because the number of skips of the product-sum operation process is maximized.

On the other hand, when the shape of the device under test 20 is not an axisymmetric shape, or when the shape of the device under test 20 is unspecified, the 3D initial image generation unit 62 converts the 2D image 24 into the 2D image. A three-dimensional initial image 26 formed by sweeping in a direction orthogonal to the plane 27 (Z direction) may be created. As a result, the initial value x _j ⁽⁰⁾ of the voxel Bx _j including the device under test 20 is not set to zero by mistake.
FIG. 8B is a perspective view showing a second example of the three-dimensional initial image 26. In this example, as indicated by an arrow, a non-zero region 24NZ of the two-dimensional plane 27 is drawn in the Z-axis direction to create a columnar three-dimensional initial image 26. As a result, the three-dimensional shape including the DUT 20 having a general shape is specified, and the voxel Bx having the three-dimensional shape in part or in whole can be set in the non-zero region 26NZ.

Further, in this method, as the three-dimensional initial image 26 created in the three-dimensional initial image generation process (step S13), the axially symmetric shape shown in FIG. 8A and the columnar shape shown in FIG. It may be selectable.
That is, the three-dimensional initial image generation unit 62 of this embodiment includes a designation input unit (input operation unit 51) that receives input information that designates the type of the three-dimensional initial image 26 to be created, and an input that is received by the input operation unit 51. Based on the information, either a three-dimensional initial image 26 formed by axially rotating the two-dimensional image 24 or a three-dimensional initial image 26 formed by sweeping the two-dimensional image 24 in a direction orthogonal to the two-dimensional plane. You may create it.
Thus, when the DUT 20 can be approximated to an axially symmetric shape, the three-dimensional initial image 26 is set to an axially symmetric shape, and when the DUT 20 is other shapes, the three-dimensional initial image 26 is set to a columnar shape. Can do.

In the case of the two-dimensional image creation process (step S12) of the present method, the voxel value x _j constituting the two-dimensional image 24 is set in the zero region 24ZA because the crystal pair 14 facing each other with the voxel Bx interposed therebetween. This is a case where the measured value of the measured projection data y _i is zero for all of the above. That is, as shown in FIG. 2, since there is the radiation source 21 between the scintillation detector 11A1 and 11B, the measured projection data _{y i} related to the crystal pair 14 _i becomes non-zero, of the crystal pair 14 _i The voxel Bx existing at the intermediate position belongs to the non-zero region 24NZ in the two-dimensional image 24.

Meanwhile, since the radiation source 21 is not present between the scintillation detector 11A2 and 11B shown in FIG. 2, the measured projection data y _i related to the crystal pair 14 _i is zero, present in the intermediate position of the crystal pair 14 _i The voxel Bx _j to be assigned belongs to the zero region 24ZA in the two-dimensional image 24.
In the case of this method, the initial value x _j ⁽⁰⁾ of an arbitrary voxel Bx _j constituting the two-dimensional image 24 is set to zero because the measurement is performed on all crystal pairs 14 sandwiching the voxel Bx _j. This is a case where the projection data y _i is zero.
That is, when the measured projection data y _i applied to a plurality of different crystal pairs 14 _i corresponding to the voxel Bx are all equal to or less than a predetermined threshold value, at least one of the first or second sum-of-products operation processing is performed on the voxel Bx _j . Execution is skipped.

Therefore, even if no photon happens to be detected by stochastic fluctuation in one of the crystal pairs 14 sandwiching an arbitrary voxel Bx, a predetermined threshold value is set in any of the other crystal pairs 14 _i sandwiching the voxel Bx _j. When the above light quantity is detected, a non-zero value is set as the initial value x _j ⁽⁰⁾ of the voxel Bx. For this reason, the arithmetic processing regarding the voxel Bx which should perform a front projection step (step S141) and a back projection step (step S142) is not skipped accidentally.

<Image reconstruction processing>
The voxel value x _j is estimated by the statistical image reconstruction method using the detection probability C _ij , the coefficient C _j, and the actually measured projection data y _i for the voxel Bx _j for which the initial value x _j ⁽⁰⁾ is set.
In the voxel Bx _j to which zero is given as the initial value x _j ⁽⁰⁾ , the forward projection value p _i calculated by the above equation (2) is always zero.
Furthermore, since the voxel Bx _j given zero as the initial value x _j ⁽⁰⁾ is caused by the fact that the corresponding actually measured projection data y _i is zero, the back projection value calculated by the above equation (3) D _{j is} also zero.
Therefore, at the beginning of step S141 related to the first product-sum operation, the forward projection calculation unit 33 performs zero determination of the voxel Bx _j . If voxel Bx _j is zero, both the first product-sum operation (step S141) and the second product-sum operation (step S142) in FIG. 6 are skipped.

Here, in order to reconstruct a high-resolution three-dimensional image, it is required to reduce the voxel size.
In this method, in the image reconstruction process (steps S14 and S15), the iterative calculation unit 32 repetitively executes the calculation process by reducing the size of the voxel Bx in a plurality of stages.
The iterative operation unit 32 performs the first or second product-sum operation process on the voxel Bx _j after the refinement included in the voxel Bx corresponding to the zero region 26ZA in the three-dimensional initial image 26 before the refinement. Skip execution of at least one of

That is, in this method, the iterative operation of the rough reconstruction process (step S14) and the iterative operation of the fine reconstruction process (step S15) are performed using different analysis models 25.
FIG. 9A is a perspective view showing the voxel arrangement of the analysis model 25 (schematic model 251) used for the rough reconstruction process. FIG. 5B is a perspective view showing the voxel arrangement of the analysis model 25 (fine model 252) used for the fine reconstruction process.
The voxel Bx (voxel Bx2) in the fine model 252 is obtained by halving the side length of the voxel Bx (voxel Bx1) in the general model 251 in each axial direction.
Therefore, the fine model 252 is configured with eight times the number of voxels as compared with the schematic model 251, and can reconstruct a three-dimensional image of the internal structure of the DUT 20 with a resolution twice that in each axial direction. .
A total of eight voxels Bx2 constituting the fine model 252 together with other voxels adjacent in the respective axial directions correspond to the individual voxels Bx1 of the schematic model 251. That is, voxel Bx1 and voxel Bx2 correspond by 1: N (N is an integer of 2 or more).

  Further, in the fine model 252, the inside of the detection region 22 is divided by 2m voxels in each axial direction. Then, four voxels correspond to each scintillation detector 11.

In this way, a voxel value x _j of the voxels Bx1 obtained by executing the first and second product-sum operation at a high speed by schematic model 251, as the first and the initial value of the second product-sum operation by definition model 252 Give to voxel Bx2.
Therefore, regarding the voxel Bx2 included in the voxel Bx1 to which zero is given as the initial value x _j ⁽⁰⁾ of the rough reconstruction process, zero is continuously given as the initial value of the fine reconstruction process.
Thus, for the voxel Bx2 included in the zero region 26ZA of the three-dimensional initial image 26 corresponding to the crystal pair 14 _i for which the actual projection data y _i is zero, the first and second product sums in the fine reconstruction process are performed. Since the calculation is skipped, the total time required for the fine reconstruction process can be suppressed.

In the fine model 252, four voxels correspond to one scintillation detector 11. Here, regarding the actual measurement projection data y _i shown in the above equation (3), the projection value generated by the projection data generation unit 31 and subjected to the rough reconstruction process may be used as it is, or corresponds to the voxel Bx2. Interpolated data obtained by interpolating the projection values may be used.
For detection probability C _ij and coefficient C _j , values calculated in advance corresponding to voxel Bx2 of fine model 252 are used.

In this method, the case where the voxel Bx1 of the general model 251 and the voxel Bx2 of the fine model 252 correspond by 1: N is exemplified, but the present invention is not limited to this. That is, the voxel Bx1 of the general model 251 and the voxel Bx2 of the fine model 252 may be associated with each other by M: N (M is an integer of 2 or more smaller than N).
In this case, the voxel values x _j for each voxel Bx1 obtained in schematic reconstruction process, and during the interpolation appropriately based on the positional relationship of the voxel Bx1 and voxel Bx2, may determine the initial value of the voxel Bx2.

  Incidentally, all or part of the functional blocks 31 to 36 of the image processing apparatus 3 of FIG. 1 may be realized by hardware such as a semiconductor integrated circuit, or may be recorded on a recording medium such as a nonvolatile memory or an optical disk. It may be realized by a programmed program or program code. Such a program or program code causes a computer having a processor such as a CPU to execute all or part of the processing of the functional blocks 31 to 36.

  In addition, the data storage unit 38 and the parameter storage unit 39 perform recording and reading of data to and from a recording medium such as a volatile memory or a non-volatile memory (for example, a semiconductor memory or a magnetic recording medium). And a circuit or a program for this purpose. The data storage unit 38 and the parameter storage unit 39 may be configured in advance on a predetermined storage area of the recording medium, or may be configured on an appropriate storage area assigned when the image processing apparatus 3 is operated. Also good.

  As mentioned above, although embodiment of this invention was described with reference to drawings, these are the illustrations of this invention, Various structures other than the above are also employable.

For example, in the above embodiment, the three-dimensional initial image 26 is created in order to set the initial value x _j ⁽⁰⁾ of the image reconstruction process, but the present invention is not limited to this. For example, the image processing apparatus 3, the measured projection data y _i is stores all the number i of voxels Bx sandwiched crystal pairs 14 _i is equal to or less than the threshold, the first or second product-sum operation on the voxel Bx _j Processing may be skipped.

In the above embodiment has illustrated the manner of setting zero as the initial value x _{j ⁽⁰⁾} of the voxel values x _j, it is not limited thereto. A predetermined specific flag value is given as an initial value x _j ⁽⁰⁾ to a predetermined voxel Bx, and the front projection calculation unit 33 or the rear projection calculation unit 34 determines whether the voxel value x _j matches the flag value. You may do.

  Moreover, in the said embodiment, although a pair of detection plate 10A, 10B was used and the detection surface of these detection plate 10A, 10B was a plane, it is not limited to these. Instead of the detection unit 2 shown in FIG. 1, a detection device having a single detection surface 12 or three or more detection surfaces 12 distributed continuously may be used. Further, the shape of the detection surface 12 is not limited to a flat surface and may be a curved surface.

In the above embodiment, the image processing apparatus 3 is applied to the opposed PET apparatus, but is not limited to this. The image processing apparatus 3 can also be applied to other detection apparatuses such as an X-ray CT apparatus.
In the case of an X-ray CT apparatus, the detection unit 2 does not necessarily include a pair of detection plates 10 facing each other. In this case, the amount of X-rays irradiated from the light source and transmitted through the DUT 20 is measured by the detection plate 10 at a plurality of positions. Then, by detecting a detector element having a known incident intensity with respect to the measurement object 20 and a detection intensity substantially equal to the detection plate 10, the measurement object 20 is positioned between the position of the detection plate 10 and the light source. It is known that it does not exist. Then, by setting the initial value x _j ⁽⁰⁾ of the voxel Bx arranged in such a non-existing region to zero, the same effect as in the above embodiment can be obtained for the X-ray CT apparatus.

  In the above embodiment, the ML-EM method is used as the statistical image reconstruction method. However, if the statistical image reconstruction method uses the forward projection data or the backward projection data using the detection probability. Well, it is not limited to the ML-EM method. For example, an improved version of the ML-EM method such as the OS-EM method or the MAP-EM (Maximum A Posteriori-Expectation Maximization) method can be used.

Further, the measured projection data y _i used for the image reconstruction process may be measured and acquired by the detection unit 2 as in the above embodiment, or may be taken from the storage medium into the data storage unit 38 of the image processing apparatus 3. You may remember.

DESCRIPTION OF SYMBOLS 1 Image reconstruction system 2 Detection part 3 Image processing apparatus 10 Detection plate 11 Scintillation detector 12 Detection surface 14 Crystal pair 20 Object 21 Source 22 Detection area 23 Cancer cell 24 Two-dimensional image 24NZ Non-zero area 24ZA Zero area 25 Analysis model 251 Schematic model 252 Fine model 26 Three-dimensional initial image 26 NZ Non-zero region 26 ZA Zero region 27 Two-dimensional plane 28 High-density region 31 Projection data generation unit 32 Repetition calculation unit 33 Forward projection calculation unit 34 Back projection calculation unit 35 Density calculation Unit 36 display control unit 37 memory 38 data storage unit 39 parameter storage unit 40 control unit 41 drive mechanism 50 display 51 input operation unit 60 initial value setting unit 61 two-dimensional image generation unit 62 three-dimensional initial image generation unit AX rotation axis Bx voxel x voxel value y measured projection data

Claims (15)

Projection data indicating the amount of photons flying from the object to be measured to each of the plurality of detector elements, an analysis model in which a plurality of voxels are three-dimensionally arranged inside the detection region where the object to be measured is disposed, And reconstructing a three-dimensional image of the internal structure of the object to be measured by repeatedly executing a calculation process for estimating a voxel value for each voxel from the projection data using a statistical image reconstruction method. An iterative operation unit;
a parameter storage unit storing data indicating a detection probability C _{ij at} which the photon that has passed through the j th voxel is detected by the i th detector element;
An image processing apparatus comprising:
The iterative operation unit is
(A) a first product-sum operation process of the detection probability C _ij read from the parameter storage unit and the voxel value of the j-th voxel;
A second product-sum operation process of the detection probability C _ij and the projection data applied to the i-th detector element;
And performing the arithmetic processing including
(B) An image processing apparatus that skips execution of at least one of the first or second product-sum operation processing for the voxel corresponding to the detector element whose projection data is equal to or less than a predetermined value.   The image processing apparatus according to claim 1, wherein the iterative operation unit sets a voxel value of the voxel corresponding to a detector element whose projection data is equal to or less than the predetermined value to zero.   The image processing apparatus according to claim 2, further comprising: an initial value setting unit that gives zero as an initial value of the voxel value to the voxel corresponding to a detector element whose projection data is equal to or less than the predetermined value.   And when the projection data applied to a plurality of different detector elements corresponding to the voxel is less than or equal to the predetermined value, the execution of at least one of the first or second product-sum operation processing is skipped for the voxel. Item 4. The image processing device according to any one of Items 1 to 3. A zero region that is a set of voxels corresponding to the detector elements whose projection data is less than or equal to the predetermined value, and a non-zero that is a set of voxels corresponding to the detector elements whose projection data is greater than the predetermined value A two-dimensional image generation unit that generates a two-dimensional image of voxel values, including a region,
Generating a three-dimensional initial image including the two-dimensional image, and determining a voxel corresponding to the zero region in the three-dimensional initial image;
Further comprising
The said iterative operation part skips execution of at least one of said 1st or 2nd product-sum operation process about the voxel corresponding to the said zero area | region in the said three-dimensional initial image. Image processing device.   The image processing apparatus according to claim 5, wherein the three-dimensional initial image generation unit generates the three-dimensional initial image formed by rotating the two-dimensional image.   The image processing apparatus according to claim 5, wherein the three-dimensional initial image generation unit generates the three-dimensional initial image obtained by sweeping the two-dimensional image in a direction orthogonal to the two-dimensional plane. The three-dimensional initial image generation unit
A designation input unit for receiving input information for designating a type of the three-dimensional initial image to be generated;
Based on the input information received by the designation input unit, the three-dimensional initial image obtained by axially rotating the two-dimensional image, or the two-dimensional image is swept in a direction orthogonal to the two-dimensional plane. The image processing apparatus according to claim 5, wherein any one of a three-dimensional initial image is generated.   The image processing apparatus according to claim 5, wherein a plurality of voxel values in the non-zero region of the two-dimensional image are different from each other according to the value of the projection data. The image processing apparatus according to any one of claims 5 to 9, wherein the arithmetic processing is repeatedly performed by reducing the size of the voxel in a plurality of stages.
For the voxels after miniaturization included in the voxels corresponding to the zero region in the three-dimensional initial image before miniaturization, the iterative operation unit includes at least one of the first or second product-sum operation processing An image processing apparatus characterized by skipping execution. The image processing apparatus according to any one of claims 1 to 10, comprising a pair of scintillation detectors arranged opposite to each other as the detector element.
The image reconstruction system, wherein the scintillation detector detects a pair of photons emitted from the object to be measured on both sides in a predetermined direction.   The image reconstruction system according to claim 11, wherein the image reconstruction system includes a pair of opposing detection surfaces in which a plurality of scintillation detectors are respectively arranged, and an interval between the detection surfaces is variable. An image for reconstructing a three-dimensional image of the internal structure of the object to be measured using a statistical image reconstruction method based on projection data indicating the amount of photons flying from the object to be measured to each of the plurality of detector elements. A processing method,
(A) a model generation step of generating an analysis model by three-dimensionally arranging a plurality of voxels within a detection region in which the object to be measured is arranged;
(B) A forward projection step of performing a first product-sum operation process of a detection probability C _{ij in} which the photon that has passed through the j-th voxel is detected by the i-th detector element and a voxel value of the j-th voxel When,
(C) a back projection step of performing a second product-sum operation process of the detection probability C _ij and the projection data applied to the i-th detector element;
(D) an iterative operation step of repeatedly executing the forward projection step and the backward projection step to estimate a voxel value for each voxel;
Including
An image processing method characterized by skipping execution of at least one of the first or second product-sum operation processing for the voxel corresponding to the detector element whose projection data is equal to or less than a predetermined value.   The image processing method according to claim 13, further comprising an actual measurement step of measuring a light quantity of the photon in a state where the object to be measured is pressed and widened. An image for reconstructing a three-dimensional image of the internal structure of the object to be measured using a statistical image reconstruction method based on projection data indicating the amount of photons flying from the object to be measured to each of the plurality of detector elements. A program for causing an information processing apparatus to execute a reconfiguration process,
The image reconstruction process includes:
(A) a model generation process for generating an analysis model by three-dimensionally arranging a plurality of voxels within a detection region in which the object to be measured is arranged;
(B) A forward projection process in which a first product-sum operation process is performed between the detection probability C _{ij at} which the photon that has passed through the j-th voxel is detected by the i-th detector element and the voxel value of the j-th voxel. When,
(C) a back projection process for performing a second product-sum operation process on the detection probability C _ij and the projection data applied to the i-th detector element;
(D) an iterative calculation process of repeatedly executing the forward projection process and the backward projection process to estimate a voxel value for each voxel;
Including
A program that skips execution of at least one of the first and second sum-of-products operations for the voxel corresponding to the detector element whose projection data is less than or equal to a predetermined value.

Images