WO2016035706A1 - 画像化装置及び方法 - Google Patents

画像化装置及び方法 Download PDF

Info

Publication number
WO2016035706A1
WO2016035706A1 PCT/JP2015/074427 JP2015074427W WO2016035706A1 WO 2016035706 A1 WO2016035706 A1 WO 2016035706A1 JP 2015074427 W JP2015074427 W JP 2015074427W WO 2016035706 A1 WO2016035706 A1 WO 2016035706A1
Authority
WO
WIPO (PCT)
Prior art keywords
event
pixel
detector
compton
compton scattering
Prior art date
Application number
PCT/JP2015/074427
Other languages
English (en)
French (fr)
Inventor
敬弘 猪田
信治 本村
秀一 榎本
Original Assignee
国立研究開発法人理化学研究所
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 国立研究開発法人理化学研究所 filed Critical 国立研究開発法人理化学研究所
Priority to JP2016546612A priority Critical patent/JP6607576B2/ja
Priority to EP15837605.3A priority patent/EP3190431B1/en
Priority to US15/509,112 priority patent/US9784857B2/en
Publication of WO2016035706A1 publication Critical patent/WO2016035706A1/ja

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01TMEASUREMENT OF NUCLEAR OR X-RADIATION
    • G01T1/00Measuring X-radiation, gamma radiation, corpuscular radiation, or cosmic radiation
    • G01T1/29Measurement performed on radiation beams, e.g. position or section of the beam; Measurement of spatial distribution of radiation
    • G01T1/2914Measurement of spatial distribution of radiation
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01TMEASUREMENT OF NUCLEAR OR X-RADIATION
    • G01T1/00Measuring X-radiation, gamma radiation, corpuscular radiation, or cosmic radiation
    • G01T1/16Measuring radiation intensity
    • G01T1/161Applications in the field of nuclear medicine, e.g. in vivo counting
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01TMEASUREMENT OF NUCLEAR OR X-RADIATION
    • G01T1/00Measuring X-radiation, gamma radiation, corpuscular radiation, or cosmic radiation
    • G01T1/16Measuring radiation intensity
    • G01T1/161Applications in the field of nuclear medicine, e.g. in vivo counting
    • G01T1/164Scintigraphy
    • G01T1/1641Static instruments for imaging the distribution of radioactivity in one or two dimensions using one or several scintillating elements; Radio-isotope cameras
    • G01T1/1647Processing of scintigraphic data

Definitions

  • the present invention relates to an imaging apparatus and method.
  • a method using a Compton camera has been proposed as a method for imaging the spatial distribution of a gamma ray source (see, for example, Patent Document 1 below).
  • the Compton camera has a wider imaging field of view and a wider detection energy band than the existing nuclear medicine technologies positron emission tomography (PET) and single photon emission computed tomography (SPECT).
  • PET positron emission tomography
  • SPECT single photon emission computed tomography
  • the Compton camera is not only an improvement of the existing technology, but also enables simultaneous multi-molecule imaging that simultaneously images a wide variety of radionuclides, which has been difficult in the past. This is a technique capable of noninvasively visualizing the behavior of a plurality of radiodiagnostic drugs and biofunctional molecules, and can acquire information on life functions that are more advanced than conventional techniques. Therefore, the Compton camera is expected to contribute to early diagnosis in life science research and clinical medicine.
  • a typical Compton camera includes two radiation detectors capable of measuring an interaction position and energy as a front-stage detector and a rear-stage detector.
  • An imaging device using a Compton camera pays attention to a gamma ray detection event in which a gamma ray from a gamma ray source is Compton scattered by a front stage detector and then photoelectrically absorbed by a back stage detector.
  • the interaction information that is, the detection position and detection energy of the gamma ray at each detector
  • the spatial distribution of the gamma ray source is imaged.
  • a three-dimensional image showing the spatial distribution of the gamma ray source is called a distribution image.
  • the pixel value ⁇ j of the pixel j of the distribution image is obtained according to the following equation (1) (see, for example, Non-Patent Document 1 below).
  • the distribution of the estimated gamma ray source is brought closer to the true distribution by updating the pixel value by iterative calculation.
  • ⁇ j (l) and ⁇ j (l + 1) indicate pixel values of the pixel j obtained by the 1st and (l + 1) th iterative calculations, respectively.
  • s j is a detection sensitivity parameter for the pixel j, and is composed of the geometric efficiency of the previous detector viewed from the pixel j. Division by the detection sensitivity parameter s j means correction of the pixel value in consideration of detection sensitivity. The significance of other parameters in equation (1) will be described in detail later.
  • the detection sensitivity is determined depending only on the pixel j.
  • the response function changes greatly for each event, so it is considered desirable to optimize the correction of detection sensitivity for each event.
  • the correction of the pixel value with the appropriate detection sensitivity contributes to the generation of an appropriate distribution image.
  • an object of the present invention is to provide an imaging apparatus and method that contribute to the optimization of a distribution image showing the spatial distribution of a gamma ray source.
  • An imaging apparatus includes a post-stage detector, a pre-stage detector disposed between a gamma ray source and the post-stage detector, and gamma rays that are Compton scattered by the pre-stage detector.
  • An event detection means for detecting an event photoelectrically absorbed by the detector, and a distribution image showing a spatial distribution of the gamma ray source based on the measurement data of the interaction between each detector and the gamma ray in a plurality of events.
  • Computing means for generating an image in an image space containing the image, and the computing means includes a probability parameter indicating a probability that Compton-scattered gamma rays have arrived from the image space in each event. Individually set for each of a plurality of pixels in the distribution image based on the measurement data of the event, and using the set probability parameter And generating the distribution image.
  • the probability that the gamma ray has arrived from within the image space is Since it can be said that it is relatively low, the measurement data of the event should not be strongly reflected in the distribution image. And the state of the said intersection changes for every event. Therefore, a probability parameter indicating the probability may be set for each event. This reduces artifacts. Further, if an appropriate probability parameter is set for each pixel using the interaction measurement data, further reduction of artifacts is expected and generation of an appropriate distribution image is expected.
  • An imaging apparatus includes a post-stage detector, a pre-stage detector disposed between a gamma ray source and the post-stage detector, and gamma rays Compton scattered by the pre-stage detector.
  • An event detection means for detecting an event photoelectrically absorbed by the detector, and an operation for generating a distribution image showing a spatial distribution of the gamma ray source based on measurement data of interaction between each detector and gamma rays in a plurality of events.
  • An imaging device comprising: means for calculating a detection sensitivity parameter indicating Compton scattering detection sensitivity at each event based on the measurement data of each event for each of the plurality of pixels in the distribution image. The distribution image is generated individually using the set detection sensitivity parameter.
  • An imaging method is a method used for a Compton camera including a post-stage detector and a pre-stage detector disposed between a gamma ray source and the post-stage detector. Based on the measurement data of the interaction between each detector and the gamma ray at each event in which the gamma rays that are Compton scattered by the detector are photoelectrically absorbed by the subsequent detector, a distribution image showing the spatial distribution of the gamma ray source is obtained.
  • a probability parameter indicating a probability that Compton-scattered gamma rays have arrived from the image space in each event is set based on the measurement data of each event.
  • Each of the plurality of pixels is individually set, and the distribution image is generated using the set probability parameter.
  • An imaging method is a method used in a Compton camera including a rear-stage detector and a front-stage detector disposed between a gamma-ray source and the rear-stage detector. Generates a distribution image showing the spatial distribution of the gamma ray source based on the measurement data of the interaction between each gamma ray and each detector at each event in which the gamma rays scattered by the Compton scattering at the detector are photoelectrically absorbed by the latter detector.
  • the detection sensitivity parameter indicating the detection sensitivity of Compton scattering at each event is individually set and set for each of the plurality of pixels in the distribution image based on the measurement data at each event.
  • the distribution image is generated using the detection sensitivity parameter.
  • an imaging apparatus and method that contribute to the optimization of a distribution image showing the spatial distribution of a gamma ray source.
  • FIG. 4 is a diagram showing how measurement data of a plurality of events are obtained. These are the relationship diagrams between each detector and the image space.
  • A) And (b) is a figure which shows a mode that an image space is divided
  • FIG. 6 is a diagram illustrating an example of a state of intersection between a Compton conical surface and an image space of two events.
  • (A) And (b) is a figure which shows the example of the intersection state of the Compton conical surface area
  • (A) And (b) is a figure which shows the example of the intersection state of the Compton conical surface area
  • (A) And (b) is a comparison figure of a reference
  • FIG. 11 is a flowchart illustrating a distribution image generation operation based on measurement data of each event. These are the schematic top views of the cancer bearing mouse
  • reference numerals 11 and 12 denote a front-stage detector and a rear-stage detector that form the Compton camera 10, respectively.
  • Each of the detectors 11 and 12 is a radiation detector capable of measuring an interaction position and energy, and is formed of a semiconductor or a light emitting substance (scintillator).
  • Reference numeral 201 denotes a gamma ray source that emits gamma rays.
  • Reference numeral 200 denotes an imaging target including the gamma ray source 201.
  • the detectors 11 and 12 are typically arranged in parallel and spaced from each other. When viewed from the imaging target 200 that includes the gamma ray source 201, the front detector 11 and the rear detector 12 are arranged in this order. That is, the front detector 11 is disposed between the imaging target 200 including the gamma ray source 201 and the rear detector 12. In FIG. 1B, only one gamma ray source 201 is shown in the imaging target 200, but many gamma ray sources 201 can be scattered in the imaging target 200.
  • a gamma ray from the gamma ray source 201 is incident on the front-stage detector 11 and Compton scattered at a position R 1 in the front-stage detector 11, and the scattered gamma ray is incident on the rear-stage detector 12 and is positioned in the rear-stage detector 12. and it is photoelectrically absorbed by R 2.
  • Compton scattering and photoelectric absorption are a kind of interaction between gamma rays and a detector. Further, it is assumed that all energy of gamma rays after Compton scattering is absorbed at a position R 2 in the post-stage detector 12. At this time, the energy detected by the detectors 11 and 12 is represented by E 1 and E 2, respectively.
  • the energy E 1 is the energy given to the electron at the position R 1 by the gamma rays from the gamma ray source 201 (that is, the energy lost by the Compton scattering in the energy of the gamma rays from the gamma ray source 201).
  • the energy E 2 is the total energy of gamma rays after Compton scattering, and is the photoelectric absorption energy at the position R 2 . Then, the following equation (A1) is established according to the kinematics of Compton scattering.
  • ⁇ C is the Compton scattering angle
  • me is the static mass of electrons
  • c is the speed of light in vacuum
  • E 0 is the initial energy of gamma rays.
  • the Compton scattering angle ⁇ C is an angle formed by the flight direction (traveling direction) of gamma rays before Compton scattering and the flight direction (traveling direction) of gamma rays after Compton scattering.
  • the initial energy E 0 of gamma rays is the initial energy of gamma rays emitted from the gamma ray source 201 and may be known to the gamma ray source distribution imaging apparatus 1 described later.
  • a parameter is known to the imaging apparatus 1 means that the imaging apparatus 1 recognizes the value of the parameter in advance.
  • an antibody labeled with a predetermined gamma-ray emitting nucleus is administered to a living body as the imaging target 200.
  • the value of the energy of gamma rays emitted from the gamma ray emission nucleus as the gamma ray source 201 may be recognized in advance by the gamma ray source distribution imaging apparatus 1.
  • the gamma ray source 201 may be a positron emitting nucleus.
  • the gamma rays accompanying the annihilation of positrons emitted from the nuclei, which are positron emitting nuclei, by positron decay become the gamma rays from the gamma ray source 201.
  • the imaging device 1 creates a gamma-ray energy spectrum during imaging by the Compton camera 10. It is possible to recognize the initial energy E 0 of the gamma ray from its energy spectrum.
  • the conical surface 202 has the Compton scattering position R 1 as the apex, the Compton scattering angle ⁇ C as the half apex angle (half apex angle; open angle), and on the straight line of the flight path of gamma rays after Compton scattering. (That is, on a straight line passing through the positions R 1 and R 2 ) is a conical surface having a central axis (hereinafter also referred to as a Compton conical surface).
  • a position R 1 be a conical surface (i.e.
  • An event in which a set of “R 1 , R 2 , E 1 and E 2 ” is detected is called an event. That is, in one event, one set of “R 1 , R 2 , E 1 and E 2 ” is detected by one-time Compton scattering and total energy absorption of gamma rays after Compton scattering. Two Compton conical surfaces are required. In one event, we can only know that there is a gamma ray source somewhere on the drawn Compton cone, but we measure a lot of Compton scattering phenomena, and multiple Compton cones for multiple events are displayed in the image space. Drawing above suggests that there are gamma ray sources at the intersection of many Compton cones.
  • FIG. 2 is a perspective view showing the configuration of the Compton camera 10.
  • the Compton conical surface 211C is set for the gamma ray 211 detected in one event
  • the Compton conical surface 212C is set for the gamma ray 212 detected in another event.
  • a similar Compton conical surface is additionally set. Since the position where more Compton conical surfaces overlap corresponds to the position where many gamma rays are generated, for example, by setting the pixel value of each pixel in the image space according to the multiplicity of overlap of the Compton conical surfaces, The state of the source distribution can be imaged.
  • Each of the detectors 11 and 12 is a radiation detector capable of measuring the interaction position between itself and the gamma ray and the energy of the gamma ray.
  • movement of the detectors 11 and 12 are mainly demonstrated for the example of the structure in which each of the detectors 11 and 12 is formed using the semiconductor flat plate (planar type detector) which consists of high purity germanium.
  • the one surface and the other surface of the semiconductor flat plate are referred to as a front surface and a rear surface, respectively.
  • the normal direction of each of the front surface and the rear surface coincides with the thickness direction of the semiconductor flat plate.
  • the front surface is located closer to the imaging target 200 including the gamma ray source 201 than the rear surface.
  • the front and rear surfaces of the front detector 11 refer to the front and rear surfaces of the semiconductor flat plate forming the front detector 11, and the front and rear surfaces of the rear detector 12 are the front and rear surfaces of the semiconductor flat plate forming the rear detector 12. Point to.
  • a plurality of strip electrodes arranged in a stripe shape are arranged on each of the front surface and the rear surface.
  • the X axis, Y axis, and Z axis that are orthogonal to each other are defined.
  • a plurality of strip electrodes extending in the Y-axis direction are arranged along the X-axis direction on the front surface
  • a plurality of strip electrodes extending in the X-axis direction are arranged on the rear surface on the Y-axis. It is arranged along the direction.
  • the strip electrodes adjacent to each other are electrically insulated through spacing. The same applies to the post-stage detector 12.
  • each detector 11 and 12 on the front side uppermost n + layer of semiconductor plates uppermost rear side are in the p + layer, a region sandwiched between the n + layer and p + layer For example, a high-purity p-type semiconductor is used.
  • each detector is a detector using germanium, each detector is cooled to a liquid nitrogen temperature ( ⁇ 196 ° C.), for example. For this reason, in each semiconductor flat plate of the detectors 11 and 12, supply of conduction carriers due to thermal excitation of electrons from the valence band is suppressed.
  • the detectors 11 and 12 may be formed of a semiconductor that does not require cooling for measurement.
  • a reverse bias voltage is applied to each strip electrode on both sides so as to generate an electric field in the direction of the thickness of the semiconductor flat plate.
  • a region called a depletion layer in which supply of conductive carriers from the impurity level is suppressed is formed in the semiconductor, and the semiconductor is in a high resistance state.
  • the gamma rays enter the depletion layer, the gamma rays interact with electrons in the semiconductor, and these electrons receive energy from the gamma rays.
  • a large number of carrier charges are generated by excitation of electrons from the valence band along the path of the electrons that have received energy, and the generated carrier charges are extracted by the reverse bias voltage.
  • the current due to the carrier charge has a current value reflecting the energy of the electrons interacting with the gamma rays, and the charge amount obtained by integrating the current is proportional to the energy. For this reason, when Compton scattering occurs in the detector 11 as the above interaction, the detector 11 detects a current corresponding to the energy E 1 from which the gamma rays have been lost due to Compton scattering.
  • the energy E 1 is specified by the amount of charge corresponding to the integral value of the current detected by the detector 11.
  • gamma rays after Compton scattering are absorbed by the detector 12, a photoelectric effect occurs as the above interaction, and a current corresponding to the total energy of gamma rays after Compton scattering is detected by the detector 12.
  • the energy E 2 is specified by the amount of charge corresponding to the integral value of the current detected by the detector 12.
  • an interaction position is detected by specifying a combination of front and rear strip electrodes where a current signal due to the interaction between the gamma ray and the detector is strongly detected.
  • the interaction position detected by the upstream detector 11 is the Compton scattering position R 1
  • the interaction position detected by the downstream detector 12 is the photoelectric absorption position R 2 .
  • the block (x A , y A ) when a current signal is detected from the x A th strip electrode on the front surface and a current signal is detected from the y A th strip electrode on the rear surface, the block (x It can be detected that an interaction with the gamma ray has occurred at ( A , y A ) (the same applies to the detector 12).
  • the position of the block (x A , y A ) In the X axis direction, the position of the block (x A , y A ) is the same as the position of the x A th strip electrode on the front surface, and in the Y axis direction, the position of the block (x A , y A ) is The position is the same as the position of the y A th strip electrode.
  • x A is an arbitrary natural number equal to or less than the number of strip electrodes provided on the front surface
  • y A is an arbitrary natural number equal to or less than the number of strip electrodes provided on the rear surface.
  • the number of strip electrodes arranged on each surface of the front detector 11 is arbitrary (the same applies to the rear detector 12).
  • the detection resolution of the interaction position in the X-axis direction depends on the number of strip electrodes arranged on the front surface
  • the detection resolution of the interaction position in the Y-axis direction depends on the number of strip electrodes arranged on the rear surface.
  • the Z-axis component of the interaction position detected in each detector may be constant (center position of each detector in the Z-axis direction). However, in each of the detectors 11 and 12, the interaction position in the Z-axis direction (Z-axis component of the interaction position) is based on the time difference between the current signal from the front strip electrode and the current signal from the rear strip electrode. May be derived with high accuracy.
  • the semiconductor forming the detectors 11 and 12 may be a semiconductor other than germanium as long as it exhibits sensitivity to gamma rays.
  • silicon, cadmium telluride, cadmium zinc telluride, or diamond may be used as the semiconductor forming the detectors 11 and 12.
  • each of the detectors 11 and 12 may be any radiation detector that can measure the interaction position between itself and the gamma rays and the energy of the gamma rays. You may have such a structure.
  • each of the detectors 11 and 12 may be formed by a scintillation detector or a TPC (time projection chamber) using gas or liquid.
  • FIG. 3 is a schematic overall configuration diagram of the gamma ray source distribution imaging apparatus 1 according to the present embodiment.
  • the imaging apparatus 1 includes respective parts referred to by reference numerals 20, 30 and 41 to 44. However, it can be considered that the Compton camera 10 is not included in the components of the imaging apparatus 1.
  • the upstream detector 11 When an interaction with a gamma ray occurs inside the upstream detector 11, the upstream detector 11 outputs a signal 11 OUT corresponding to the content of the interaction.
  • the downstream detector 12 interacts with gamma rays inside itself, the downstream detector 12 outputs a signal 12 OUT corresponding to the content of the interaction.
  • the interaction position (for example, position R 1 ) with the gamma ray in the detector 11 by the signal 11 OUT and the interaction energy (for example, energy E 1) given to the detector 11 from the gamma ray by the interaction in the detector 11. ) Is identified.
  • the interaction position (for example, position R 2 ) with the gamma ray in the detector 12 by the signal 12 OUT and the interaction energy (for example, energy E 2) given from the gamma ray to the detector 12 by the interaction in the detector 12. ) Is identified.
  • the event measurement data acquisition unit 20 includes a simultaneous event extraction unit 21 and a total energy absorption determination unit 22, and event measurement data that is interaction data between the detectors 11 and 12 and gamma rays in each event is used. Generate and get.
  • FIG. 4 is a flowchart for acquiring event measurement data by the acquisition unit 20.
  • the simultaneous event extraction unit 21 outputs corresponding to events that occur simultaneously in the detectors 11 and 12 out of the output signal 11 OUT of the upstream detector 11 and the output signal 12 OUT of the downstream detector 12.
  • a set of signals is extracted as a simultaneous event detection signal.
  • the simultaneous event extracting unit 21 each time from the previous stage detector 11 signal 11 OUT is output, stored in the storage unit 44 a signal 11 OUT after having granted a time stamp indicating its output time to the signal 11 OUT with previously, each time from the subsequent detector 12 signal 12 OUT is outputted, allowed to store the signal 12 OUT after having granted a time stamp indicating its output time to the signal 12 OUT in the storage unit 44.
  • the simultaneous event extracting unit 21 reads the arbitrary signal 11 OUT stored in the storage unit 44, and stores the signal 12 OUT whose time difference from the read signal 11 OUT is within a predetermined allowable time difference. Search from within 44 (step S11).
  • step S11 When the signal 12 OUT whose time difference from the read signal 11 OUT is within a predetermined allowable time difference exists in the storage unit 44 (Y in step S11), the simultaneous detection of the read signal 11 OUT and the searched signal 12 OUT is detected. It extracts as a signal (step S12).
  • Step S15 When the signal 12 OUT whose time difference from the read signal 11 OUT is within the predetermined allowable time difference does not exist in the storage unit 44 (N in Step S11), the process proceeds to Step S15.
  • the second signal whose time difference from the first signal is within a predetermined allowable time difference is associated with the time stamp indicating a time whose time difference from the time indicated by the time stamp corresponding to the first signal is within the predetermined allowable time difference. The second signal is indicated.
  • the simultaneous event extraction unit 21 extracts the simultaneous event detection signal in real time from the output signals 11 OUT and 12 OUT of the detectors 11 and 12 without the storage of the signals 11 OUT and 12 OUT by the storage unit 44.
  • step S13 the total energy absorption determination unit 22 performs a total energy absorption determination process on the extracted simultaneous event detection signal. In the total energy absorption determination process, it is determined whether the signal 12 OUT in the simultaneous event detection signal indicates total energy absorption. Specifically, the total energy absorption determination unit 22 determines the simultaneous event determination formula “E 0 ⁇ E 0 ⁇ E for the interaction energies E 1 and E 2 indicated by the signals 11 OUT and 12 OUT in the simultaneous event detection signal.
  • ⁇ E 0 is introduced into the simultaneous event determination formula.
  • ⁇ E 0 has a positive predetermined value. For example, when E 0 is 511 keV (kiloelectron volts), ⁇ E 0 is set to 5 keV.
  • the total energy absorption determination unit 22 determines that the signal 12 OUT in the time event detection signal does not indicate total energy absorption (N in step S13), and discards the simultaneous event detection signal. Then, the process proceeds to step S15. This is because the Compton scattering angle cannot be correctly estimated from the measurement data if the simultaneous event determination formula does not hold even if Compton scattering occurs.
  • step S14 the event measurement data acquisition unit 20 recognizes that the extracted simultaneous event detection signal represents a detection signal in a valid event, and event measurement data based on the signals 11 OUT and 12 OUT in the valid event. (Interaction measurement data) is acquired, and the acquired event measurement data is stored in the storage unit 44. After step S14, the process proceeds to step S15.
  • One event measurement data includes information identifying interaction positions R 1 and R 2 and energies E 1 and E 2 in one event.
  • the acquisition unit 20 detects an event using the extraction unit 21 and the determination unit 22 and extracts event measurement data from the output signals 11 OUT and 12 OUT of the detectors 11 and 12.
  • step S15 the event measurement data acquisition unit 20 determines whether or not a predetermined measurement data acquisition end condition is met. If the measurement data acquisition end condition is satisfied, the event measurement data acquisition process ends. If the measurement data acquisition end condition is not satisfied, the process returns to step S11. For example, the measurement data acquisition end condition is satisfied when the number of event measurement data acquisition reaches a predetermined number (for example, several 10000). Alternatively, for example, when the process related to the extraction of the simultaneous event detection signal is completed for all the signals 11 OUT stored in the storage unit 44, the measurement data acquisition end condition is satisfied.
  • a predetermined number for example, several 10000.
  • the main arithmetic unit 30 includes a highly parallel calculation accelerator such as a CPU (Central Processing Unit) and a SIMD (single instruction multiple data) arithmetic unit.
  • the main calculation unit 30 sets a predetermined image space in the space including the imaging target 200, and based on event measurement data generated for a plurality of events, a distribution image (gamma ray source distribution image) indicating the density distribution of the gamma ray source As an image in the image space. A method for generating the distribution image will be described later.
  • the main calculation unit 30 may be formed using a dedicated calculation device for generating a reconstructed image including a distribution image.
  • the dedicated arithmetic device is a device having a hardware configuration optimum for image reconstruction, which is constructed using, for example, ASIC (application specific integrated circuit) or FPGA (field-programmable gate array).
  • the overall control unit 41 is realized by using a CPU (Central Processing Unit) or the like, and comprehensively controls the operation of each part in the imaging apparatus 1.
  • the display unit 42 includes a liquid crystal display panel and the like, and displays an arbitrary image including a distribution image under the control of the overall control unit 41.
  • the operation unit 43 includes a pointing device, a keyboard, and the like, and accepts arbitrary instructions and operations from the user of the imaging apparatus 1.
  • the storage unit 44 includes a ROM (Read Only Memory) and a RAM (Random Access Memory).
  • the storage unit 44 stores arbitrary programs as well as programs to be operated by the main arithmetic unit 30, the overall control unit 41, and the like.
  • the storage unit 44 may include a secondary storage device such as an HDD (Hard disk drive) or a flash memory.
  • FIG. 5 shows how event measurement data for a plurality of events are acquired one after another.
  • the event measurement data in the i-th event (that is, the i-th event) is also referred to as i-th event measurement data.
  • i is an arbitrary natural number.
  • the detected energies E 1 and E 2 at the detectors 11 and 12 in the measurement data of the i-th event are represented by E i1 and E i2, respectively.
  • the interaction position at the detector 11 ie, the detected Compton scattering position
  • R 1 and the interaction position (ie, the detected photoelectric absorption position) R 2 at the detector 12 in the measurement data of the i-th event In particular, R i1 and R i2 respectively.
  • the main calculation unit 30 calculates the Compton scattering angle ⁇ C and the Compton conical surface using the above formula (A1) based on the event measurement data (E i1 , R i1 , E i2 , R i2 ) in the event. Can be derived and set.
  • FIG. 6 shows an image space IS set and defined by the main arithmetic unit 30.
  • the image space IS may be set in a dedicated memory device that can be included in the storage unit 44 or the like.
  • the image space IS is a three-dimensional virtual space set in association with the real space, and includes at least the position where the imaging target 200 exists.
  • the image space IS has a finite size.
  • the shape of the image space IS may be arbitrary, but here, the image space IS is a rectangular parallelepiped space having two planes parallel to the XY plane, the YZ plane, and the ZX plane.
  • the XY plane is a plane parallel to the X axis and the Y axis
  • the YZ plane is a plane parallel to the Y axis and the Z axis
  • the ZX plane is a plane parallel to the Z axis and the X axis.
  • one of the planes in the image space IS parallel to the XY plane is adjacent to the front surface of the detector 11, and the center of the detector 11, the center of the detector 12, and the center of the image space IS are aligned. Deploy.
  • the main calculation unit 30 can generate a distribution image representing the density distribution state of the gamma ray source as a three-dimensional image in the image space IS.
  • FIG. 7A shows a state where the image space IS is equally divided into a plurality of parts in each of the X, Y, and Z axis directions.
  • the image space IS is formed from element areas equally divided into a plurality of parts in each of the X, Y, and Z axis directions, and each element area functions as a pixel of a distribution image (see FIG. 7B).
  • Each pixel forming the distribution image is a three-dimensional pixel having a size in each of the X, Y, and Z axis directions, and can be called a voxel.
  • each pixel forming the distribution image has a cubic shape.
  • a certain focused pixel is expressed as a pixel j or a pixel k using an arbitrary integer j or k indicating a pixel number.
  • the position of the pixel j is represented by the symbol r j (see FIG. 8).
  • FIG. 8 shows the relationship between the Compton conical surface of the i-th event, the pixel j, the Compton scattering position R i1, and the like.
  • the position r j of the pixel j represents the center position of the pixel j (the center coordinate of the pixel j in the image space IS).
  • a pixel refers to a pixel in the image space IS (and thus a pixel in the distribution image).
  • the main arithmetic unit 30 can generate a distribution image of the gamma ray source by an image reconstruction method based on a list-mode maximum-likelihood expectation-maximization (LM-ML-EM) method.
  • LM-ML-EM list-mode maximum-likelihood expectation-maximization
  • lambda j represented by lambda j (l) or lambda j (l + 1) indicates the pixel value of the pixel j in the distribution image (the j th pixel distribution image).
  • the pixel value of the pixel j represents an expected value of the number of radionuclides that have collapsed at the pixel j within the imaging time of the Compton camera 10. Since the number of decays of radionuclide per unit time follows a Poisson distribution, ⁇ j is an expected value of this Poisson distribution.
  • the radionuclide here is a gamma ray emitting nucleus or positron emitting nucleus as the gamma ray source 201.
  • the pixel value is updated by iterative calculation to bring the gamma ray source distribution closer to the estimated true distribution.
  • ⁇ j (l) indicates the pixel value of the pixel j obtained in the first iteration (the above expected value for pixel j), and ⁇ j (l + 1) is the (l + 1) th iteration.
  • the pixel value of the pixel j obtained by (the above expected value for the pixel j) is shown. Accordingly, the symbol “l” associated with ⁇ j indicates the number of iterations (thus an integer).
  • t ij is called a system response function, and represents the probability that the gamma ray emitted from the pixel j is measured as the i-th event.
  • the system response function t ij is set for pixels that intersect the Compton conical surface estimated from the event measurement data.
  • Y i represents the number of gamma rays detected by the Compton camera 10 in the i-th event.
  • Y i 1”.
  • the actual number of detected gamma rays (corresponding to “ ⁇ k t ik ⁇ k ” of the denominator) obtained by calculation (corresponding to the number of detected gamma rays) (Corresponding to “Y i ”).
  • the pixel value is updated by iterative calculation based on this ratio.
  • the pixel value is updated in the iterative calculation so that the calculated gamma ray detection number (denominator) approaches the actually measured gamma ray detection number (numerator).
  • the pixel value is updated so that the calculated gamma ray detection number becomes smaller than before the update, and the calculated gamma ray detection number becomes the actual gamma ray detection number. If it is smaller than the number, the pixel value is updated so that the number of detected gamma rays is larger than before the update.
  • the detection sensitivity parameter s j is a detection sensitivity parameter for pixel j, and division by s j means correction of the pixel value in consideration of detection sensitivity.
  • the detection sensitivity parameter s j mainly includes the geometric efficiency of the detector 11 viewed from the pixel j, the interaction probability between the gamma ray and the detectors 11 and 12, and the time from when the gamma ray is emitted until it is detected. Consists of factors such as attenuation. Guests can set the detection sensitivity parameter s j through a numerical calculation using these factors may be obtained detection sensitivity parameter s j by Monte Carlo simulation.
  • the geometric efficiency of the detector 11 viewed from the pixel j represents the size (solid angle) of the detector 11 viewed from the pixel j.
  • the probability that the gamma ray emitted from the pixel j hits the detector 11 depends on the geometric efficiency of the detector 11 viewed from the pixel j. The closer the distance between the pixel j and the detector 11 is, the larger the size of the detector 11 viewed from the pixel j is, and the probability that the gamma rays emitted from the pixel j hit the detector 11 increases.
  • the geometric efficiency of the detector 11 viewed from the pixel j is approximately proportional to the inverse square of the distance between the pixel j and the detector 11.
  • Formula (B1a) and Formula (B1b) are a part of Formula (B1).
  • i and k indicate an event number and a pixel number, respectively.
  • Expression (B1a) represents the total sum of “t ik ⁇ k (l) ” for all pixels in the i-th event.
  • i and j represent an event number and a pixel number, respectively.
  • Expression (B1b) represents the sum of all the fractional events having expression (B1a) as the denominator and “Y i t ij ” as the numerator. All events refer to all of the events used for obtaining the pixel value of each pixel of the distribution image.
  • the arrival direction of the gamma ray toward the upstream detector 11 is limited to the generatrix direction of the Compton conical surface.
  • the Compton scattering angle is determined by the energy at Compton scattering, it has uncertainty derived from a physical phenomenon during gamma ray detection.
  • the distribution function of the Compton scattering angle (the distribution function of the error distribution of the Compton scattering angle) showing this uncertainty is called a scattering angle distribution function.
  • the Compton conical surface also has an uncertainty according to the scattering angle distribution function, and it is estimated that the gamma ray source exists stochastically in the spread derived from the scattering angle distribution function.
  • a function f O ( ⁇ ) represented by the following formula (C1) is defined as a scattering angle distribution function.
  • FIG. 9 shows a graph of the scattering angle distribution function f O ( ⁇ ).
  • Z represents the atomic number of a semiconductor material that constitutes the detectors 11 and 12 and interacts with gamma rays (hereinafter referred to as detector material), and n l represents the detection thereof. This represents the number of electrons in the l-th electron shell in the organic material.
  • the summation operation symbol “ ⁇ ” in equation (C1) indicates the summation for all electron shells in the detector material.
  • ⁇ e represents a standard deviation obtained from the energy detection resolution when the detector 11 detects the energy E 1 .
  • ⁇ d, l represents the standard deviation of the detected energy due to Doppler broadening due to the electron momentum of the l th electron shell in the detector material.
  • parameters other than ⁇ C , ⁇ , ⁇ e and ⁇ d, l have predetermined values known to the imaging apparatus 1.
  • ⁇ C represents the Compton scattering angle obtained by the above formula (A1).
  • represents the Compton scattering angle when it is assumed that the gamma rays emitted from the target position NP are Compton scattered at the interaction position R 1 in the detector 11. That is, the angle formed by the straight line connecting the interaction positions R 1 and R 2 of the detectors 11 and 12 (the central axis of the Compton conical surface) and the straight line connecting the target position NP and the position R 1 corresponds to the angle ⁇ . To do.
  • the angle ( ⁇ C ) is the Compton scattering angle ⁇ obtained when the gamma rays emitted from the target position NP are Compton scattered at the interaction position R 1 in the detector 11 and the Compton obtained by the above formula (A1).
  • the error angle ⁇ with respect to the scattering angle ⁇ C is represented.
  • the intensity of the scattering angle distribution function f O ( ⁇ ) (that is, the value of the function f O ( ⁇ )) is maximum when the error angle ⁇ is zero, and is the absolute value of the error angle ⁇ . As it increases, it decreases toward zero.
  • a back projection operation is performed to project the conical gamma ray arrival direction estimated from the gamma ray measurement data onto the image space IS.
  • the conical gamma ray arrival direction estimated from the gamma ray measurement data includes an error based on the above-described uncertainty, it is necessary to project a cone having a thickness due to the error in the back projection operation. Since the cone of thickness to be considered in the backprojection is due uncertainty Compton scattering angle may be determined the thickness from the scattering angle distribution function f O ( ⁇ ).
  • a Compton conical surface having a thickness is called a Compton conical surface region.
  • the main calculation unit 30 can set the Compton conical surface area by giving the thickness to the Compton conical surface for each event as follows.
  • the integrated values when the function f O ( ⁇ ) is integrated at the angle ⁇ are referred to as INT A and INT B, respectively.
  • INT A and INT B are referred to as INT A and INT B, respectively.
  • a positive angle ⁇ REF at which the ratio “INT B / INT A ” is close to 1 and becomes a predetermined positive value (for example, 95%) less than 1 is obtained.
  • FIG. 11 shows three conical surfaces 311 to 313 considered for one event.
  • the conical surfaces 311 to 313 are all centered on the straight line of the flight path of the gamma ray after Compton scattering (that is, on the straight line passing through the positions R 1 and R 2 ) with the measured Compton scattering position R 1 as the vertex.
  • the conical surface 312 is a Compton conical surface having an angle ⁇ C as a half apex angle (half apex angle; opening angle), while the conical surfaces 311 and 313 each have an angle ( ⁇ C ⁇ REF ) as a half.
  • the main calculation unit 30 sets a region between the lower limit conical surface 311 and the upper limit conical surface 313 as the Compton conical surface region.
  • the Compton conical surface 312 is included in the region between the lower limit conical surface 311 and the upper conical surface 313, that is, the Compton conical surface region.
  • ( ⁇ REF ⁇ 2) corresponds to the above-described thickness expressed as an angle.
  • the thickness of the Compton conical surface at a portion where the distance from the position R 1 is d corresponds to the length of an arc having a radius d and a central angle ( ⁇ REF ⁇ 2).
  • a first improvement method relating to generation of a distribution image will be described.
  • the first improvement method and each improvement method to be described later are also based on the LM-ML-EM method as in the reference method described above.
  • a distribution image is generated based on the following formula (D1).
  • t ij ⁇ j represents an expected value measured by gamma rays emitted from the pixel j in which the i-th event has an accumulation amount (in other words, an abundance) of radionuclides corresponding to the pixel value ⁇ j .
  • v ij and s ij will be described later.
  • the significance of the sum ( ⁇ ) in formula (D1) will be supplementarily described.
  • the formula (D1a) and the formula (D1b) are a part of the formula (D1).
  • i and k indicate an event number and a pixel number, respectively.
  • Expression (D1a) represents the sum of “t ik ⁇ k (l) ” for all pixels in the i-th event.
  • i and j represent an event number and a pixel number, respectively.
  • Equation (D1b) represents the sum of all fractional events having the product of equation (D1a) and s ij as the denominator and “v ij t ij ” as the numerator. All events refer to all of the events used for obtaining the pixel value of each pixel of the distribution image.
  • a cone 360 is a cone having a Compton conical surface derived based on the measurement data of the first event as a side surface and the length of the generatrix as d A
  • the cone 370 is the measurement data of the second event. Is a cone having a Compton conical surface derived from the above as a side surface and a bus bar length d A. Near the bottom surface of the cone 360, the cone 360 is within the image space IS, while a portion of the cone 370 protrudes from the image space IS near the bottom surface of the cone 370. In the illustration of the cone 370 in FIG.
  • a portion where the cone 370 and the image space IS intersect (that is, a portion of the cone 370 located in the image space IS) is indicated by a solid line, and the cone 370 and the image are displayed.
  • a portion where the space IS does not intersect is indicated by a broken line (corresponding to a broken line 371).
  • a spherical surface having the Compton scattering position R 11 in the first event and having a radius d A (hereinafter referred to as a first spherical surface), and having a Compton scattering position R 21 in the second event and having a radius d A.
  • a second spherical surface Is assumed to be a spherical surface (hereinafter referred to as a second spherical surface).
  • reference numeral 380 shows a curved surface that is the intersection of the first spherical surface and the image space IS projected onto the XY plane.
  • 390 shows a curved surface that is the intersection of the second spherical surface and the image space IS projected on the XY plane.
  • FIG. 14A Although the shape of the projection of the intersection of the first spherical surface and the image space IS on the XY plane may not actually be a quadrangle, in FIG. 14A, it is assumed that the shape is a quadrangle for simplicity. is doing. The same applies to FIG. 14B.
  • the thickness of the Compton conical surface based on the measurement data of the first event is given a thickness based on the scattering angle distribution function, so that the Compton conical surface region of the first event is set.
  • a hatched area 381 is an intersection area between the Compton conical surface area of the first event and the curved surface 380 (that is, an intersection area of the Compton conical surface area of the first event, the first spherical surface, and the image space IS).
  • a hatched area 391 is an intersection area between the Compton conical surface area of the second event and the curved surface 390 (that is, an intersection area of the Compton conical surface area of the second event, the second spherical surface, and the image space IS).
  • a dot region 392 represents a region that does not intersect the image space IS among the intersecting regions of the second event Compton conical surface region and the second spherical surface.
  • the ratio of the area intersecting the image space IS in the intersecting area between the Compton conical surface area and the first spherical surface of the first event is 100%.
  • a distance d A from Compton scattering position R 11 in and first event pixels belonging to the region 381 is a 100 and uniform intensity of the scattering angle distribution function in the region 381
  • the probability is approximately 1/100, and the system response function t based on such estimation is assumed. 1j is set.
  • the ratio of the region intersecting the image space IS in the intersecting region between the Compton conical surface region and the second spherical surface of the second event is 50%. Therefore, the number of pixels belonging to the region 391 is 50, the gamma ray source in the second event is located at a distance d A from the Compton scattering position R 21 , and the intensity of the scattering angle distribution function in the region 391 is uniform. Assuming that the gamma ray source of the second event exists in each pixel belonging to the region 391, the probability that the second event gamma ray source exists is approximately “1/50”, and the system response function t based on such estimation is assumed. 2j is set.
  • the regions 391 and 392 are present in the image space IS, it should be estimated that gamma rays have arrived from approximately 100 pixels in each event, and therefore belong to the region 391.
  • the probability that a second event gamma source is present in each pixel should be estimated to be approximately "1/100" equally.
  • the gamma ray source of the second event is assigned to each pixel belonging to the region 391. It is estimated that there is a probability of “1/50”.
  • a cone 410 is a cone having a Compton conical surface derived from the measurement data of the i-th event as a side surface and the length of the generatrix as d A
  • the cone 420 is the measurement data of the i-th event.
  • the events corresponding to the cones 410 and 420 are common, and since “0 ⁇ d A ⁇ d B ”, the conical surface of the cone 410 is part of the conical surface of the cone 420.
  • the cone 410 is within the image space IS, while a portion of the cone 420 protrudes from the image space IS near the bottom surface of the cone 420.
  • a portion where the cone 420 intersects the image space IS is indicated by a solid line, and the cone 420 and the image are displayed.
  • a portion where the space IS does not intersect is indicated by a broken line (corresponding to the broken line 421).
  • a sphere having a Compton scattering position R i1 in the i-th event and having a radius d A (hereinafter referred to as a spherical surface A), and a Compton scattering position R i1 in the i-th event and having a radius d B in the center.
  • a spherical surface (hereinafter referred to as spherical surface B) is assumed.
  • reference numeral 440 indicates a curved surface that is the intersection of the spherical surface A and the image space IS projected onto the XY plane.
  • reference numeral 450 denotes a curved surface that is the intersection of the spherical surface B and the image space IS projected onto the XY plane.
  • the shape of the projected portion of the intersection of the spherical surface A and the image space IS on the XY plane may not actually be a quadrangle, in FIG. 16A, it is assumed that the shape is a quadrangle for simplicity. ing. The same applies to FIG.
  • the thickness of the Compton conical surface based on the measurement data of the i-th event is given a thickness based on the scattering angle distribution function, whereby the Compton conical surface region of the i-th event is set.
  • a hatched area 441 is a crossing area between the Compton conical surface area of the i-th event and the curved surface 440 (that is, a crossing area of the Compton conical surface area of the i-th event, the spherical surface A, and the image space IS). Represents.
  • FIG. 16A a hatched area 441 is a crossing area between the Compton conical surface area of the i-th event and the curved surface 440 (that is, a crossing area of the Compton conical surface area of the i-th event, the spherical surface A, and the image space IS).
  • a hatched area 451 is a crossing area between the Compton conical surface area of the i-th event and the curved surface 450 (that is, a crossing area of the Compton conical surface area of the i-th event, the spherical surface B, and the image space IS).
  • a dot region 452 represents a region that does not intersect the image space IS among the intersecting regions of the Compton conical surface region of the i-th event and the spherical surface B.
  • the ratio of the area intersecting the image space IS in the intersection area between the Compton conical surface area of the i-th event and the spherical surface A is 100%. Therefore, the number of pixels belonging to the region 441 is N A , the gamma ray source in the i-th event is located at a distance d A from the Compton scattering position R i1 , and the intensity of the scattering angle distribution function in the region 441 is uniform.
  • the ratio of the area intersecting the image space IS among the intersecting areas of the Compton conical surface area and the spherical surface B of the i-th event is 50%. Therefore, the number of pixels belonging to the region 451 is N B , the gamma ray source at the i-th event is at a position away from the Compton scattering position R i1 by the distance d B , and the intensity of the scattering angle distribution function in the region 451 is uniform.
  • the probability parameter v ij (refer to the above formula (D1)) introduced in the first improvement method serves as a regularization factor for suppressing such a phenomenon (overfitting).
  • v ij is defined by the following formula (D2).
  • R i1 is an interaction position (ie, Compton scattering position) between the gamma ray detected in the i-th event and the pre-stage detector 11 (see FIG. 8).
  • r j indicates the position (center position) of the pixel j.
  • r represents an arbitrary position (coordinates) in the image space IS.
  • V k represents the volume of the pixel k.
  • “R ⁇ V k ” represents a position r included in the pixel k (in other words, coordinates r included in the volume of the pixel k).
  • v ij is p for all pixels k intersecting a sphere centered on the Compton scattering position R i1 of the i-th event and having a radius between the position (r j ) of pixel j and the distance between Compton scattering positions R i1. Represents the sum of ik .
  • the probability parameter v ij indicates the probability that gamma rays that have caused Compton scattering at the i-th event have arrived from within the image space IS.
  • the probability parameter is set for each pixel, and v ij is a probability parameter set for the pixel j. More specifically, v ij is gamma rays caused the Compton scattering and radius the distance between the position (r j) and Compton scattering position R i1 of and the pixel j around the Compton scattering position R i1 at the i event When it is assumed that the position has arrived from a position on the spherical surface, the probability that the gamma ray has arrived from within the image space IS is shown.
  • p ik is the gamma ray passing through the pixel k among the gamma-ray arrival directions estimated for the i-th event. Represents the percentage of the arrival direction. That is, p ik represents the probability density of the gamma ray arrival direction passing through the pixel k when the gamma ray arrival direction of the i-th event is regarded as a random variable.
  • v ij is the sum of p ik for all the pixels k including the position r satisfying “
  • ”. ik represents the contribution of pixel k to v ij .
  • p ik is based on the uncertainty of the Compton scattering angle and is defined by the equation (D3).
  • p ij is defined using variable j instead of variable k in p ik .
  • the angle ⁇ i (r) in the formula (D3) is defined by the formula (D3a).
  • f O, i represents the scattering angle distribution function (the distribution function of the error distribution of the Compton scattering angle) at the i-th event.
  • f O, i is regarded as a function of the angle ⁇ i (r) depending on the position r.
  • the position r represents an arbitrary position in the image space IS as described above, and may be considered to correspond to the target position NP in FIG.
  • ⁇ i (r) represents the Compton scattering angle of the gamma ray when it is assumed that the i-th event is caused by the gamma ray emitted from the position r.
  • R ii and R i2 are the interaction position of the gamma ray detected at the i-th event and the detector 11 (ie, the detected Compton scattering position) and the gamma ray detected at the i-th event and the detector 12 respectively.
  • Interaction position ie, detected photoelectric absorption position
  • V j represents the volume of pixel j.
  • the position r is used as an integration variable, and “f O, i [ ⁇ i (r)] / ⁇
  • ⁇ C in the above formula (C1) is calculated using the formula (A1) for the i-th event. Compton scattering angle.
  • ⁇ in the above equation (C1) represents ⁇ i (r), and positions R i1 and R i2 Is equivalent to the angle formed by the straight line connecting the two (the central axis of the Compton conical surface) and the straight line connecting the positions R i1 and r.
  • f O, i [ ⁇ i (r)] represents the probability density of the error distribution of the Compton scattering angle at the position r, and p ij is obtained by integrating this over the volume V j of the pixel j. .
  • simple integration does not take into account the geometric efficiency of the detector 11 as viewed from the pixel j.
  • f O, i is a function of angle, it is necessary to consider a solid angle due to a difference in angle. Therefore, in the equation (D3), f O, i [ ⁇ i (r)] is expressed as the reciprocal of the geometric efficiency approximation “
  • the value of the function f O, i for pixels that do not belong to the Compton conical surface region is considered to be sufficiently small. Therefore, the right side of the equation (D2) is for all the pixels k that include the position r satisfying “
  • the main arithmetic unit 30 uses the probability parameter v ij indicating the probability that the Compton-scattered gamma rays arrived from the image space IS in each event as the measurement data (R i1 , R i2 , E i1 , E i2 ) are set for each pixel of the distribution image by calculation, and a distribution image is generated using the probability parameter v ij set for each pixel (that is, the pixel value of each pixel is obtained). ).
  • the pixels (target pixels) for which the probability parameter v ij is to be set based on the measurement data at each event need not be all the pixels of the distribution image (however, they may be all the pixels of the distribution image). .
  • the main calculation unit 30 individually calculates the probability parameter v ij for each of a plurality of target pixels in the distribution image by calculation based on measurement data at each event. It can be said that a distribution image is generated using the set probability parameter v ij (that is, a pixel value of each pixel is obtained).
  • the main calculation unit 30 sets the probability parameter v ij for each event and for each target pixel. That is, the main calculation unit 30 pays attention to each of a combination of a plurality of events and a plurality of target pixels, and sets the probability parameter v ij for each combination.
  • the main processing unit 30 When attention is paid to the i-th event and one pixel j (when attention is paid to the combination of the i-th event and pixel j), the main processing unit 30 is centered on the detected Compton scattering position R i1 and the position of the focus pixel j (R j ) and a sphere having a radius (eg, d A or d B in FIG. 15) between the Compton scattering position R i1 , the Compton conical region of the i-th event of interest, and the intersection of the image space IS Based on the state, a probability parameter v ij is set for the pixel of interest j of the i-th event that has been noticed.
  • the hatched area 441 in FIG. 16A or the hatched area 451 in FIG. 16B corresponds to the intersecting area of the spherical surface, the Compton conical surface area, and the image space IS, and FIGS. These crossing states are different from each other.
  • the main calculation unit 30 when focusing on the i-th event and one pixel j, the main calculation unit 30 is centered on the detected Compton scattering position R i1 and the position (r j ) and Compton scattering position of the target pixel j.
  • each pixel arranged in the intersection region is extracted, and an index (p ij ) corresponding to the intensity of the scattering angle distribution function (f O, i ) is derived for each extracted pixel (formula (D3 )reference). Then, based on the sum of the indices derived for each extracted pixel (see Expression (D2)), the probability parameter v ij for the pixel of interest j of the i-th event of interest is set.
  • v ij is considered to be an integral of the scattering angle uncertainty with the integral range being the range of the total scattering angle
  • the ideal value of v ij is 1.
  • v ij is smaller than 1 when a part of the Compton conical surface region is outside the image space IS. Therefore, an event in which the Compton conical surface region slightly intersects the image space IS has a reduced contribution to the reconstructed image (distributed image) due to the introduction of v ij , thereby reducing artifacts.
  • the intersection state is different if the target pixel is different in one event. In the first improvement method, this is taken into consideration, and the optimum v ij is set for each pixel in each event, so that the artifact reduction effect is enhanced.
  • the factor that constitutes the detection sensitivity has the greatest influence on the geometric efficiency of the detector 11 viewed from the individual pixel j.
  • the geometric efficiency is approximately the pixel j and the detector 11. It is proportional to the inverse square of the distance between.
  • the detection sensitivity (s j ) depends only on the position of the pixel j, and thus depends on the difference in Compton scattering position (interaction position in the detector 11). Differences in geometric efficiency are averaged.
  • the distance (511, 521) between the Compton scattering position and the pixel j differs between when the Compton scattering position is the position 510 and when it is the position 520.
  • the geometric efficiency Should be different, but assuming that the detection sensitivity (s j ) depends only on the position of pixel j, the difference in geometric efficiency due to the difference in Compton scattering position will be averaged. . The closer to the detector 11, the greater the influence of this averaging (because the difference between the distances 511 and 521 is likely to occur), and the accuracy of sensitivity correction using the detection sensitivity parameter is reduced.
  • the detection sensitivity parameter for each pixel is set in consideration of the Compton scattering position R i1 for each event.
  • s ij is a detection sensitivity parameter indicating the detection sensitivity of the detector 11 with respect to Compton scattering of the i-th event.
  • the detection sensitivity parameter is set for each pixel
  • s ij is a detection sensitivity parameter set for the pixel j.
  • division by s jj means correction of the pixel value in consideration of detection sensitivity.
  • FIGS. 19A and 19B show a comparison between the reference method and the first improved method.
  • the detection sensitivity parameter s ij is defined by the following equation (D4).
  • ⁇ t (E 0 ) represents the total interaction cross section of the gamma ray having the energy of E 0 and the detector material
  • d i, j1 is the detection of the gamma ray from the pixel j at the i-th event.
  • the range in the vessel 11 is represented (see FIG. 20). That is, d i, j1 represents the distance that the gamma ray from the pixel j in the i-th event has passed through the detector 11 until Compton scattering occurs.
  • exp [ ⁇ t (E 0 ) d i, j1 ] means that the gamma ray from the pixel j does not interact with the detector substance in the detector 11 until the Compton scattering position R i1 in the i th event. It represents the probability of arrival (in other words, the attenuation that occurs in the detector 11 until the gamma ray from the pixel j in the i-th event is Compton scattered).
  • the value of ⁇ t (E 0 ) is uniquely determined based on the initial energy E 0 of gamma rays emitted from the gamma ray source 201 and the physical properties of the detector material, and is known to the imaging apparatus 1.
  • the range d i, j1 is determined from the position r j of the pixel j, the Compton scattering position R i1 of the i-th event, and the shape of the detector 11.
  • the detector 11 can detect the interaction position in the X, Y, and Z axis directions in m, n, and o stages, respectively.
  • M, n and o are any integers of 2 or more).
  • the detector 11 can be considered as a collection of (m ⁇ n ⁇ o) independent detector elements.
  • each detector element has detection sensitivity for all pixels in the image space IS.
  • R i1 may be considered to represent the center position of the detector element including the Compton scattering position at the i-th event.
  • ⁇ 2 ” represents an approximation of the geometric efficiency of the detector element centered on the position R i1 as viewed from the pixel j.
  • the main calculation unit 30 performs calculation based on the detection sensitivity parameter s ij indicating the detection sensitivity of Compton scattering at each event based on the measurement data (including R i1 ) at each event. Is set for each pixel of the distribution image, and the distribution image is generated using the detection sensitivity parameter s ij set for each pixel (that is, the pixel value of each pixel is obtained).
  • the pixels (target pixels) for which the detection sensitivity parameter s ij should be set based on the measurement data at each event need not be all pixels of the distribution image (however, they may be all pixels of the distribution image). ). In each event, it is sufficient to include each pixel in the Compton conical surface region as a target pixel for which the detection sensitivity parameter s ij is to be set. That is, in each event, it is sufficient to calculate s ij based on the equation (D4) only for each pixel in the Compton conical surface area. In each event, a constant value may be substituted for s ij for each pixel outside the Compton conical surface area.
  • the main calculation unit 30 applies the detection sensitivity parameter s ij to each of a plurality of target pixels in the distribution image by calculation based on measurement data at each event. It can also be said that the distribution image is generated (that is, the pixel value of each pixel is obtained) using the set detection sensitivity parameter s ij individually.
  • the main calculation unit 30 sets the detection sensitivity parameter s ij for each event and for each target pixel. That is, the main calculation unit 30 sets the detection sensitivity parameter s ij for each combination by paying attention to each combination of a plurality of events and a plurality of target pixels.
  • the main calculation unit 30 sets the detected Compton scattering position R i1 and the position r j of the noticed pixel j to each other. Accordingly, the detection sensitivity parameter (s ij ) is set. More specifically, when attention is paid to the i-th event and one pixel j, the main calculation unit 30 calculates gamma rays before Compton scattering in the detector 11 based on the Compton scattering position R i1 and the position r j of the target pixel j.
  • Sensitivity parameter for the target pixel j of the i-th event of interest based on the range d i, j1 of the current and the distance
  • the distance between the gamma ray source and the detector 11 is relatively small, the gamma rays from the gamma ray source are likely to hit the detector 11 and thus are easily detected by the detector 11, and the distance between the gamma ray source and the detector 11 is relatively long.
  • the gamma ray from the gamma ray source is hard to hit the detector 11, so it is difficult for the detector 11 to detect it. If a distribution image is generated without considering such ease of detection / hardness of detection, the presence of a gamma ray source that is far from the detector 11 is less likely to be reflected in the distribution image.
  • the detection sensitivity parameter (s ij ) is set in consideration of not only the position of each pixel but also the interaction position R i1 in the detector 11, thereby detecting without considering the interaction position R i1.
  • the detection sensitivity parameter (s ij )
  • the detection sensitivity parameter (s ij ) is set in consideration of not only the position of each pixel but also the interaction position R i1 in the detector 11, thereby detecting without considering the interaction position R i1.
  • the sensitivity parameter (FIG. 19A)
  • more appropriate sensitivity correction can be performed.
  • the distribution image more accurately represents the distribution of the gamma ray source.
  • the gamma ray attenuation (exp [ ⁇ t (E 0 ) d i, j1 ]) in the detector 11 is also considered using the interaction position R i1 , the validity of the distribution image is considered. (The accuracy of estimating the gamma source distribution increases.)
  • d ⁇ c, ij / d ⁇ represents the probability that the gamma ray from the pixel j causes Compton scattering in the i-th event. Since this probability varies depending on the Compton scattering angle, “d ⁇ c, ij / d ⁇ ” depends on the pixel j.
  • E 0 ′ represents the energy of gamma rays after Compton scattering. Therefore, E 0 ′ in the equation (D5) corresponds to a value obtained by subtracting the detected energy E 1 at the i-th event from the initial energy E 0 of the gamma ray.
  • ⁇ t (E 0 ′) represents the total interaction cross section between the gamma ray having the energy of E 0 ′ and the detector substance
  • d i, 12 represents the gamma ray after Compton scattering at the i-th event as the detector.
  • 12 represents the range of the gamma ray detectors 11 and 12 until photoelectric absorption is performed.
  • d i, 12 represents the distance that the gamma rays after Compton scattering at the i-th event have passed through the detectors 11 and 12 until they are photoelectrically absorbed by the detector 12. Therefore, “exp [ ⁇ t (E 0 ′) d i, 12 ]” interacts with the detector substance in the detectors 11 and 12 until the gamma ray after Compton scattering at the i-th event reaches the photoelectric absorption position R i2. (In other words, attenuation that occurs in the detectors 11 and 12 before the gamma ray after Compton scattering in the i-th event is photoelectrically absorbed by the detector 12 at the position R i2 ).
  • ⁇ 2 represents an approximation of the geometric efficiency (solid angle) of the photoelectric absorption position R i2 viewed from the Compton scattering position R i1 .
  • ⁇ p (E 0 ′) represents the photoelectric absorption cross section between the gamma ray having the energy of E 0 ′ and the detector material, that is, the probability that the gamma ray after Compton scattering causes photoelectric absorption.
  • the main calculation unit 30 can derive and set the system response function t ij for each event and for each pixel based on the measurement data (R i1 , R i2 , E i1 , E i2 ) for each event.
  • the measurement data R i1 , R i2 , E i1 , E i2
  • D5 the value of the physical quantity that does not depend on the measurement data of the i-th event is known to the imaging apparatus 1.
  • FIG. 22 is a flowchart for generating a distribution image based on the measurement data of each event.
  • a procedure for generating a distribution image based on measurement data of each event will be described.
  • the main calculation unit 30 substitutes a predetermined common initial value (for example, 1) for the pixel value ⁇ k (0) of all the pixels of the distribution image, and sets a variable l indicating the number of iterations to zero. substitute.
  • the main arithmetic unit 30 assigns 1 to the variable i in the subsequent step S22, and then sequentially executes the processes of steps S23 to S28.
  • the main operation unit 30 reads the measurement data (R i1 , R i2 , E i1 , E i2 ) of the i-th event from the storage unit 44 in step S23, and measures the i-th event in step S24. Based on the data, the Compton scattering angle ⁇ C is calculated according to the above formula (A1) and the Compton conical surface is set. In subsequent step S25, the main computing unit 30 sets the Compton conical surface area of the i-th event by giving the thickness of the Compton conical surface of the i-th event based on the scattering angle distribution function f O ( ⁇ ).
  • step S26 the main calculation unit 30 regards each pixel in the Compton conic surface area of the i-th event as a target pixel, and for each target pixel, v ij , Calculate and set s ij and t ij .
  • the system response function t ij may be set to zero for each pixel (non-target pixel) outside the Compton conical surface area.
  • step S27 subsequent to step S26, the main arithmetic unit 30 adds v ij , s ij and t ij obtained for the i-th event to the sum calculation result of the equation (D1b) obtained up to the (i ⁇ 1) -th event.
  • step S28 it is confirmed whether or not the measurement data for all events has been read from the storage unit 44.
  • step S28 If there is measurement data that has not yet been read (N in step S28), 1 is added to the variable i in step S29, and the process returns to step S23, and the processes in steps S23 to S28 are repeated. If measurement data for all events has been read from the storage unit 44 (Y in step S28), the process proceeds to step S30.
  • step S30 the main calculation unit 30 performs a pixel value update process.
  • the updated pixel value ⁇ k (l + 1) is obtained according to the equation (D1). That is, the pixel value ⁇ k (l + 1) of each pixel of the distribution image is obtained by multiplying the sum total calculation result of the latest formula (D1b) obtained in step S27 by ⁇ k (l) .
  • the main arithmetic unit 30 confirms whether or not a predetermined update end condition is satisfied. If the update end condition is satisfied, the process proceeds to step S33, the distribution image having the latest pixel value ⁇ k is stored in the storage unit 44, and the operation of FIG. If the update end condition is not satisfied, 1 is added to the variable l in step S32, and then the process returns to step S22, and the processes after step S22 are repeated.
  • the update end condition is satisfied, for example, when the number of executions of the update process in step S30 reaches a predetermined number.
  • the update end condition may be satisfied when a predetermined end instruction operation is input to the operation unit 43 during the repeated execution of the processes of steps S22 to S32.
  • the amount of change in each pixel value ⁇ k by the update process in step S30 is sufficiently small (for example, the sum of the pixel values ⁇ k (l + 1) in all the pixels of the distribution image and the pixel value ⁇ k When the difference from the sum of (l) is equal to or less than a predetermined value), the update end condition may be satisfied.
  • Second Improvement Method Regarding Distribution Image Generation A second improvement method relating to the generation of the distribution image will be described.
  • a distribution image is generated based on the following formula (E1).
  • equation (E1) v ij is introduced based on the reference method, but the introduction of s ij is postponed.
  • the second improvement method is the same as the first improvement method except that s j is incorporated in the calculation formula of the pixel value ⁇ j instead of s ij .
  • the detection sensitivity parameter s j is obtained for each pixel, for example, according to the following formula (E2).
  • ⁇ t (E 0 ) represents the total interaction cross section between the gamma ray having energy of E 0 and the detector material
  • r j represents the position (center position) of the pixel j.
  • p is a detector element in the first to (m ⁇ n ⁇ o) th detector elements forming the detector 11 (see FIGS. 21A and 21B).
  • R p represents the center position of the p-th detector element.
  • d p, j1 is the length of the path located in the detector 11 in the straight path from the position r j to the position R p . That is, d p, j1 is similar to d i, j1 in FIG.
  • Equation (E2) Equation (E2) interacts with the detector substance in the detector 11 until the gamma ray from the pixel j reaches the p-th detector element. Represents the probability of reaching without.
  • ⁇ 2 represents an approximation of the geometric efficiency of the p th detector element as viewed from pixel j.
  • s j in the expression (E2) is calculated as the sum of detection sensitivity (gamma ray detection sensitivity) set for all detector elements. Is done. Since the number of detector elements is a fixed value, the sum of the detection sensitivities can be regarded as an average value of the detection sensitivities of all the detector elements.
  • the size of each detector element is the product (for example, 3 ⁇ 3 ⁇ 1 [mm 3 ]) of the X, Y, and Z axis components of the position detection accuracy of the detector 11.
  • the second improved method an advantageous effect is achieved by introducing v ij in comparison with the reference method.
  • the first improved method that also introduces s ij is more suitable for generating a more accurate distribution image than the second improved method.
  • a third improved method relating to the generation of the distribution image will be described.
  • a distribution image is generated based on the following formula (F1) or formula (F2).
  • s ij is introduced based on the reference method, but the introduction of v ij is postponed.
  • the third improvement method is the same as the first improvement method, except that Y i or v i is incorporated in the calculation formula of the pixel value ⁇ j instead of v ij .
  • Y i may be “1” in all events.
  • the probability parameter v i is set for each event, and one probability parameter v i is commonly applied to each pixel.
  • the probability parameter v i is calculated according to the following formula (F3).
  • the right side of the equation (F3) represents the sum of p ij for all pixels in the Compton conic surface area of the i-th event.
  • the third improved method an advantageous effect is achieved by introducing s ij in comparison with the reference method.
  • the first improved method that also introduces v ij is more suitable for generating a more accurate distribution image than the third improved method.
  • Using v i instead of v ij also reduces artifacts in comparison with the reference method, but as can be seen from the above description with reference to FIGS. 15 and 16 (a) and (b),
  • the first improved method of setting the optimal v ij for each pixel in each event has a higher artifact reduction effect.
  • FIG. 23 is a schematic plan view of a tumor-bearing mouse used in the experiment.
  • An imaging target 200 is obtained by administering a 64 Cu-labeled anti-HER2 antibody to a tumor bearing mouse transplanted with three tumor tissues A431, 4T1, and C6.
  • the 64 Cu-labeled anti-HER2 antibody is an antibody that contains 64 Cu and accumulates in a tumor containing HER2.
  • HER2 is most highly expressed in the tumor tissue A431, and therefore high antibody accumulation is expected against the tumor tissue A431.
  • 64 Cu in the tumor-bearing mouse functions as the gamma ray source 201
  • 511 keV (kiloelectron volt) gamma rays accompanying the annihilation of positrons emitted from 64 Cu function as the gamma rays from the gamma ray source 201.
  • the imaging time using the Compton camera 10 is 9 hours, and the measured number of events is 1.9 ⁇ 10 6 .
  • the three-dimensional distribution images obtained as a result of the experiment to the reference method, the second improvement method, and the first improvement method are referred to as three-dimensional distribution images 610, 620, and 630, respectively (see FIG. 24).
  • the number of iterations of the iterative calculation in the experiment is uniformly 60 times. That is, under the conditions of the experiment, the three-dimensional distribution images 610, 620, and 630 have ⁇ j (60) in the expressions (B1), (E1), and (D1) as the pixel value of the pixel j, respectively. It is a three-dimensional distribution image.
  • the MIP images 611, 612, and 613 are planar images obtained by projecting the three-dimensional distribution image 610 onto the XY plane, the ZX plane, and the YZ plane, respectively, using a maximum value projection method (maximum intensity projection).
  • the MIP images 621, 622, and 623 are planar images obtained by projecting the three-dimensional distribution image 620 onto the XY plane, the ZX plane, and the YZ plane, respectively, using the maximum value projection method.
  • the MIP images 631, 632, and 633 are planar images obtained by projecting the three-dimensional distribution image 630 onto the XY plane, the ZX plane, and the YZ plane, respectively, using the maximum value projection method. That is, for example, in the three-dimensional distribution image 610, the unit operation for projecting the pixel having the maximum pixel value on the XY plane from the group of pixels having the same position in the X and Y axis directions is performed on the entire three-dimensional distribution image 610. As a result, the MIP image 611 is formed on the XY plane that is the projection plane. The same applies to the MIP images 612 and 613, and the same applies to the MIP images 621 to 623 and 631 to 633.
  • VOI volume of interest
  • Regions 620 [A431], 620 [4T1], and 620 [C6] in which tumor tissues A431, 4T1, and C6 exist are set in the image 620, and tumor tissues A431, 4T1, and C6 exist in the three-dimensional distribution image 630 Areas 630 [A431], 630 [4T1], and 630 [C6] are set.
  • the size of the region set for the tumor tissue A431 (that is, the regions 610 [A431], 620 [A431], and 630 [A431]) is the same between the images 610, 620, and 630. The same applies to the regions set for the tumor tissues 4T1 and C6.
  • FIG. 27 shows the measurement results of the radioactivity of 64 Cu from each tissue of the tumor-bearing mouse after the imaging experiment.
  • FIG. 28 in the bar graph composed of the bars 651 to 653, after normalizing the measured radioactivity from the tumor tissue A431 shown by the bar 651 as 1, the measured radioactivity from the tumor tissues 4T1 and C6 is shown. Represented by bars 652 and 653, respectively.
  • the relative intensities indicated by the bars 671 to 673 in FIG. 28 represent the total values of the pixel values in the regions 630 [A431], 630 [4T1], and 630 [C6] based on the first improved method, respectively. However, the relative intensity indicated by the bars 671 to 673 is normalized by setting the relative intensity corresponding to the region 630 [A431] to 1.0.
  • the relative intensities indicated by the bars 681 to 683 in FIG. 28 represent the total values of the pixel values in the regions 620 [A431], 620 [4T1], and 620 [C6] based on the second improved method, respectively.
  • the relative intensity indicated by the bars 681 to 683 is normalized by setting the relative intensity corresponding to the region 620 [A431] to 1.0.
  • the relative intensities indicated by the bars 691 to 693 in FIG. 28 represent the total values of the pixel values in the regions 610 [A431], 610 [4T1], and 610 [C6] based on the reference method, respectively.
  • the relative intensity indicated by the bars 691 to 693 is normalized by setting the relative intensity corresponding to the region 610 [A431] to 1.0.
  • the relative intensities of the tumor tissues 4T1 and C6 viewed from the relative intensity (1.0) of the tumor tissue A431 are the measured radioactivity (intensities of the bars 652 and 653 viewed from the intensity of the bar 651). It can be said that the closer to, the better the quantitativeness of the image.
  • the strength of the bars 672 and 673 viewed from the strength of the bar 671 and the strength of the bars 682 and 683 viewed from the strength of the bar 681 are both
  • the strength of the bars 652 and 653 viewed from the strength of the bar 651 is close. That is, in comparison with the reference method, the first or second improvement method improves the quantitativeness of the image.
  • the target device that is the imaging device 1 itself or the target device that is a part of the imaging device 1 may be configured by hardware such as an integrated circuit or a combination of hardware and software. it can.
  • Arbitrary specific functions that are all or part of the functions realized in the target device may be described as a program, and the program may be stored in a memory (flash memory or the like) that can be mounted on the target device. . Then, the specific function may be realized by executing the program on a program execution device (for example, a microcomputer that can be mounted on the target device).
  • the program can be stored and fixed on an arbitrary recording medium.
  • the recording medium for storing and fixing the program may be mounted or connected to a device (such as a server device) different from the target device.
  • the imaging apparatus 1 includes event detection means for detecting an event and calculation means for generating a distribution image.
  • the event detection means is mainly realized by, for example, the event measurement data acquisition unit 20, and the calculation means is mainly realized by, for example, the main calculation unit 30 (see FIG. 3).
  • the computing means (30) is a conical surface setting means (S21; FIG. 22) for setting a Compton conical surface for each event, a conical surface area setting means (S22; FIG. 22) for setting a Compton conical surface region for each event, Probability parameter setting means (S23; FIG.
  • the sensitivity parameter setting means (S23; FIG. 22) for setting a probability parameter v ij for each pixel (target pixel) and detection for setting a detection sensitivity parameter s ij for each event and for each pixel (target pixel) It can be said that the sensitivity parameter setting means (S23; FIG. 22) is included.

Landscapes

  • Health & Medical Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Molecular Biology (AREA)
  • High Energy & Nuclear Physics (AREA)
  • General Physics & Mathematics (AREA)
  • Biomedical Technology (AREA)
  • Optics & Photonics (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Medical Informatics (AREA)
  • General Health & Medical Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Nuclear Medicine (AREA)
  • Measurement Of Radiation (AREA)

Abstract

 ガンマ線源からのガンマ線が前段検出器にてコンプトン散乱された後、後段検出器にて光電吸収されるガンマ線検出事象(イベント)に注目し、各検出器とガンマ線との相互作用の測定データを元に、ガンマ線源の空間分布を所定の画像空間内において画像化する。この際、コンプトン散乱を起こしたガンマ線が画像空間内から到来した確率を示す確率パラメータ(vij)、及び、ガンマ線の検出感度を示す検出感度パラメータ(sij)を、各イベントの測定データに基づきイベントごとに且つ画素ごとに設定し、それらのパラメータを用いて各画素の画素値(λj)を求める。

Description

画像化装置及び方法
 本発明は、画像化装置及び方法に関する。
 ガンマ線源の空間分布を画像化する方法としてコンプトンカメラを利用する方法が提案されている(例えば下記特許文献1参照)。コンプトンカメラは既存の核医学技術であるpositron emission tomography(PET)及びsingle photon emission computed tomography(SPECT) よりも広範囲な撮像視野と幅広い検出エネルギ帯域を有する。コンプトンカメラは単に既存技術の改善だけでなく、従来は困難であった、幅広い種類の放射性核種を同時に撮像する複数分子同時イメージングを可能とする。これは複数の放射性診断薬や生体機能分子の挙動を、非侵襲的に可視化することができる技術であり、従来技術よりも高度な生命機能の情報を取得することができる。そのため、コンプトンカメラは、生命科学研究及び臨床医療における早期診断などへの貢献が期待されている。
 典型的なコンプトンカメラは、相互作用位置及びエネルギを測定可能な2つの放射線検出器を前段検出器及び後段検出器として備える。コンプトンカメラを利用した画像化装置は、ガンマ線源からのガンマ線が前段検出器にてコンプトン散乱された後、後段検出器にて光電吸収されるガンマ線検出事象に注目し、各検出器とガンマ線との相互作用情報(即ち各検出器におけるガンマ線の検出位置及び検出エネルギ)を元に、ガンマ線源の空間分布を画像化する。
 ガンマ線源の空間分布を画像化する方法として、list-mode maximum-likelihood expectation-maximization(LM-ML-EM)法を元にした画像再構成法が提案されている。これは、観測されたガンマ線の測定データから空間中のガンマ線源分布を最尤推定により求めることが困難な場合でも、反復計算により尤度の期待値を最大化するようにガンマ線源分布を求めることを可能にする手法である。
 ガンマ線源の空間分布を示す三次元画像を分布画像と呼ぶ。LM-ML-EM法を利用した典型的な画像再構成法では、下記式(1)に従い、分布画像の画素jの画素値λjを求める(例えば下記非特許文献1参照)。LM-ML-EM法では、反復計算によって画素値を更新することで、推定されるガンマ線源の分布を真の分布へと近づけてゆく。λj (l)、λj (l+1)は、夫々、l回目、(l+1)回目の反復計算にて得られる画素jの画素値を示す。sjは画素jに対する検出感度パラメータであり、画素jから見た前段検出器の幾何学的効率などから構成される。検出感度パラメータsjによる除算は検出感度を考慮した画素値の補正を意味する。式(1)における他のパラメータの意義については後に詳説される。
Figure JPOXMLDOC01-appb-M000002
特許第4486623号公報
S. J. Wilderman, N. H. Clinthorne, J. A. Fessler, and W. L. Rogers, "List-mode maximum likelihood reconstruction of Compton scatter camera images in nuclear medicine",IEEE Nucl. Sci. Symp., 1998, vol. 3, pp. 1716-1720.
 上記式(1)に基づく方法では、“Σkikλk”が極めて小さい場合などにおいて、対応する画素値が非常に大きな値に更新されるオーバーフィッティングが生じる。これは、画像の辺縁部でアーチファクト(虚像)が生じる原因をとなる。アーチファクトの低減が有益であることは言うまでもない。アーチファクトの低減は、適正な分布画像の生成に寄与する。
 また、式(1)においては検出感度が画素jにのみ依存して決定されている。しかし、コンプトンカメラでは個々のイベントごとに応答関数が大きく変化するため、検出感度の補正も個々のイベントごとに最適化することが望ましいと考えられる。適正な検出感度による画素値の補正は、適正な分布画像の生成に寄与する。
 そこで本発明は、ガンマ線源の空間分布を示す分布画像の適正化に寄与する画像化装置及び方法を提供することを目的とする。
 本発明の一側面に係る画像化装置は、後段検出器と、ガンマ線源と前記後段検出器との間に配置された前段検出器と、前記前段検出器でコンプトン散乱されたガンマ線が前記後段検出器にて光電吸収されるイベントを検出するイベント検出手段と、複数のイベントにおける各検出器とガンマ線との相互作用の測定データに基づき、前記ガンマ線源の空間分布を示す分布画像を、前記ガンマ線源を内包する画像空間内の画像として生成する演算手段と、を備えた画像化装置において、前記演算手段は、各イベントでコンプトン散乱したガンマ線が前記画像空間内から到来した確率を示す確率パラメータを各イベントの前記測定データに基づき前記分布画像内の複数の画素の夫々に対して個別に設定し、設定した前記確率パラメータを用いて前記分布画像を生成することを特徴とする。
 例えば、或るイベントの測定データに基づいて円錐状に推定されるガンマ線到来方向が画像空間に対してかすめるように交差していた場合(図12参照)、ガンマ線が画像空間内から到来した確率は比較的低いと言えるため、当該イベントの測定データは分布画像に強く反映されるべきではない。そして、上記交差の状態はイベントごとに異なる。故に、上記確率を示す確率パラメータをイベントごとに設定すると良い。これによりアーチファクトの低減が図られる。また、相互作用測定データを用いて適正な確率パラメータを画素ごとに設定するようにすれば、アーチファクトの更なる低減が期待され、適正な分布画像の生成が期待される。
 本発明の他の側面に係る画像化装置は、後段検出器と、ガンマ線源と前記後段検出器との間に配置された前段検出器と、前記前段検出器でコンプトン散乱されたガンマ線が前記後段検出器にて光電吸収されるイベントを検出するイベント検出手段と、複数のイベントにおける各検出器とガンマ線との相互作用の測定データに基づき、前記ガンマ線源の空間分布を示す分布画像を生成する演算手段と、を備えた画像化装置において、前記演算手段は、各イベントでのコンプトン散乱の検出感度を示す検出感度パラメータを各イベントの前記測定データに基づき前記分布画像内の複数の画素の夫々に対して個別に設定し、設定した前記検出感度パラメータを用いて前記分布画像を生成することを特徴とする。
 各画素に対する検出感度パラメータを各イベントの相互作用測定データを考慮して設定することにより、相互作用の内容に応じた適正な検出感度パラメータを設定することが可能となる。結果、適正な分布画像の生成が期待される。
 本発明の一側面に係る画像化方法は、後段検出器及びガンマ線源と前記後段検出器との間に配置された前段検出器を備えたコンプトンカメラに利用される方法であって、前記前段検出器でコンプトン散乱されたガンマ線が前記後段検出器にて光電吸収される各イベントでの各検出器とガンマ線との相互作用の測定データに基づき、前記ガンマ線源の空間分布を示す分布画像を、前記ガンマ線源を内包する画像空間内の画像として生成する画像化方法において、各イベントでコンプトン散乱したガンマ線が前記画像空間内から到来した確率を示す確率パラメータを各イベントの前記測定データに基づき前記分布画像内の複数の画素の夫々に対して個別に設定し、設定した前記確率パラメータを用いて前記分布画像を生成することを特徴とする。
 本発明の他の側面に係る画像化方法は、後段検出器及びガンマ線源と前記後段検出器との間に配置された前段検出器を備えたコンプトンカメラに利用される方法であって、前記前段検出器でコンプトン散乱されたガンマ線が前記後段検出器にて光電吸収される各イベントでの各検出器とガンマ線との相互作用の測定データに基づき、前記ガンマ線源の空間分布を示す分布画像を生成する画像化方法において、各イベントでのコンプトン散乱の検出感度を示す検出感度パラメータを各イベントでの前記測定データに基づき前記分布画像内の複数の画素の夫々に対して個別に設定し、設定した前記検出感度パラメータを用いて前記分布画像を生成することを特徴とする。
 本発明によれば、ガンマ線源の空間分布を示す分布画像の適正化に寄与する画像化装置及び方法を提供することが可能である。
(a)及び(b)は、コンプトン散乱を利用したガンマ線源の分布推定方法を説明するための図である。 は、本発明の実施形態に係るコンプトンカメラの構成を示す斜視図である。 は、本発明の実施形態に係るガンマ線源分布画像化装置の概略的な全体構成図である。 は、本発明の実施形態に係るイベント測定データの取得動作フローチャートである。 は、複数のイベントの測定データが得られる様子を示した図である。 は、各検出器と画像空間との関係図である。 (a)及び(b)は、画像空間が3方向において分割される様子を示す図、及び、画像空間を形成する画素を示す図である。 は、コンプトン散乱に関与する位置、角度及び円錐面等を示した図である。 は、コンプトン散乱角の不確定性による散乱角分布関数を示す図である。 は、コンプトン散乱角の不確定性を説明するための図である。 は、コンプトン散乱角の不確定性に依存して考慮されるべき円錐面の厚みを説明するための図である。 は、コンプトン円錐面が画像空間に対して僅かしか交差しない状況を説明するための図である。 は、2つのイベントのコンプトン円錐面と画像空間との交差状態の例を示す図である。 (a)及び(b)は、2つのイベントのコンプトン円錐面領域と画像空間との交差状態の例を示す図である。 は、1つのイベントのコンプトン円錐面と画像空間との交差状態の例を示す図である。 (a)及び(b)は、1つのイベントのコンプトン円錐面領域と画像空間との交差状態の例を示す図である。 は、確率パラメータ(vij)の構成要素(pij)の意義を説明するための図である。 は、1つの画素と2つのコンプトン散乱位置との距離を示す図である。 (a)及び(b)は、基準方法と第1改良方法の比較図である。 は、ガンマ線の検出器内の飛程を説明するための図である。 (a)及び(b)は、検出器を複数の検出器要素の集合体と考えることができることを説明するための図である。 は、各イベントの測定データに基づく分布画像の生成動作フローチャートである。 は、実験に利用した担癌マウスの概略平面図である。 は、実験で得られた3つの三次元分布画像を示す図である。 は、実験で得られた3つの三次元分布画像の夫々から切り出された二次元分布画像(計9枚)を示す図である。 は、VOI評価に用いる領域設定の内容を説明するための図である。 は、撮像実験後における担癌マウスの各組織からの放射能の計測結果を示す図である。 は、各組織からの放射能の計測結果及びVOI評価の結果を示すグラフである。
 以下、本発明の実施形態の例を、図面を参照して具体的に説明する。参照される各図において、同一の部分には同一の符号を付し、同一の部分に関する重複する説明を原則として省略する。尚、本明細書では、記述の簡略化上、情報、信号、物理量又は部材等を参照する記号又は符号を記すことによって、該記号又は符号に対応する情報、信号、物理量又は部材等の名称を省略又は略記することがある。
[コンプトン散乱を利用したガンマ線源の分布推定方法の原理]
 図1(a)及び(b)を参照して、コンプトン散乱を利用したガンマ線源の分布推定方法を説明する。図1(a)及び(b)において、11及び12は、夫々、コンプトンカメラ10を形成する前段検出器及び後段検出器を表している。検出器11及び12の各々は相互作用位置及びエネルギを測定可能な放射線検出器であり、半導体や発光物質(シンチレーター)などで形成される。201は、ガンマ線を放射するガンマ線源を表している。200は、ガンマ線源201を内包する撮像対象を表している。検出器11及び12は、典型的には互いに離間した状態で平行に配置されている。ガンマ線源201を内包する撮像対象200から見て、前段検出器11、後段検出器12の順に並んでいる。即ち、ガンマ線源201を内包する撮像対象200と後段検出器12との間に前段検出器11が配置される。尚、図1(b)には、撮像対象200内にガンマ線源201が1つしか示されていないが、多数のガンマ線源201が撮像対象200内に散在しうる。
 ガンマ線源201からのガンマ線が前段検出器11に入射して前段検出器11内の位置R1にてコンプトン散乱し、散乱されたガンマ線が後段検出器12に入射して後段検出器12内の位置R2にて光電吸収されたとする。コンプトン散乱、光電吸収は、ガンマ線と検出器との相互作用の一種である。また、後段検出器12内の位置R2では、コンプトン散乱後のガンマ線の全エネルギが吸収されたとする。このとき、検出器11、12にて検出されるエネルギを、夫々、E1、E2にて表す。エネルギE1は、ガンマ線源201からのガンマ線がコンプトン散乱により位置R1における電子に与えたエネルギ(即ち、ガンマ線源201からのガンマ線のエネルギの内、コンプトン散乱により失われたエネルギ)である。エネルギE2は、コンプトン散乱後のガンマ線の全エネルギであって、位置R2における光電吸収エネルギである。そうすると、コンプトン散乱の運動学に従い、以下の式(A1)が成立する。
Figure JPOXMLDOC01-appb-M000003
 ここで、θCはコンプトン散乱角、meは電子の静止質量、cは真空中の光速、E0はガンマ線の初期エネルギを表す。コンプトン散乱角θCは、コンプトン散乱前のガンマ線の飛行方向(進行方向)とコンプトン散乱後のガンマ線の飛行方向(進行方向)との成す角度である。ガンマ線の初期エネルギE0は、ガンマ線源201が放射するガンマ線の初期エネルギであり、後述のガンマ線源分布画像化装置1にとって既知であって良い。或るパラメータが画像化装置1にとって既知であるとは、画像化装置1が当該パラメータの値を予め認識していることを指す。例えば、撮像対象200としての生体に、所定のガンマ線放出核により標識された抗体が投与される。ガンマ線源201としてのガンマ線放出核から放射されるガンマ線のエネルギの値を、ガンマ線源分布画像化装置1に予め認識させておけばよい。ガンマ線源201は陽電子放出核でも構わない。この場合、陽電子放出核である原子核から陽電子崩壊によって放出される陽電子の対消滅に伴うガンマ線が、ガンマ線源201からのガンマ線となる。尚、投与されたガンマ線放出核種が不明であって初期エネルギE0が画像化装置1にとって既知でなかったとしても、画像化装置1は、コンプトンカメラ10による撮像中にガンマ線のエネルギスペクトルを作成することが可能であり、そのエネルギスペクトルからガンマ線の初期エネルギE0を認識することが可能である。
 円錐面202は、コンプトン散乱位置R1を頂点とし、且つ、コンプトン散乱角θCを半頂角(頂角の半分;開角)とし、且つ、コンプトン散乱後のガンマ線の飛行経路の直線上に(即ち、位置R1及びR2を通る直線上)中心軸を持った円錐面(以下、コンプトン円錐面と呼ぶこともある)である。但し、円錐面202は、位置R1を頂点として位置R2から位置R1に向かう向きに半頂角の原点をとった円錐面(換言すれば、位置R1を頂点とし且つ位置R2から位置R1に向かう向きに位置R1より伸びる線分を中心軸とする円錐面)であるとする。コンプトン散乱位置R1及びコンプトン散乱角θCが求められたとき、前段検出器11に向かうガンマ線の到来方向は円錐面202の母線方向に限定される。
 1組の“R1、R2、E1及びE2”が検出される事象をイベントと呼ぶ。つまり、1回のイベントでは、1回のコンプトン散乱及びコンプトン散乱後のガンマ線の全エネルギ吸収により、1組の“R1、R2、E1及びE2”が検出され、その検出結果より1つのコンプトン円錐面が求められる。1回のイベントでは、描かれる1つのコンプトン円錐面上のどこかにガンマ線源があることしか分からないが、コンプトン散乱現象を多数測定し、複数回のイベント分の複数のコンプトン円錐面を画像空間上に描くことで、多くのコンプトン円錐面が交差した位置にガンマ線源があること示唆される。
 図2は、コンプトンカメラ10の構成を示す斜視図である。図2に示す例では、1つのイベントにおいて検出されたガンマ線211に対してコンプトン円錐面211Cが設定され、他のイベントにおいて検出されたガンマ線212に対してコンプトン円錐面212Cが設定されている。図2には描かれていないが、コンプトン散乱の事象が検出されるたびに、同様のコンプトン円錐面が追加設定される。コンプトン円錐面がより多く重なる位置はガンマ線が多数発生する位置に対応するため、例えば、画像空間内の各画素の画素値をコンプトン円錐面の重なりの多重度に応じたものとすることにより、ガンマ線源の分布の様子を画像化することが可能となる。
 検出器11及び12の夫々は、自身とガンマ線との相互作用位置及びガンマ線のエネルギを測定可能な放射線検出器である。以下では主として、検出器11及び12の夫々が高純度ゲルマニウムから成る半導体平板(プレナー型検出器)を用いて形成される構成を例にとり、検出器11及び12の構成及び動作を説明する。半導体平板の一方の面及び他方の面を、ここでは、夫々、前面及び後面と呼ぶ。前面及び後面の夫々の法線方向は、半導体平板の厚み方向と一致する。検出器11及び12の夫々において、前面は、後面よりも、ガンマ線源201を内包する撮像対象200の近くに位置する。前段検出器11の前面及び後面は、前段検出器11を形成する半導体平板の前面及び後面を指し、後段検出器12の前面及び後面は、後段検出器12を形成する半導体平板の前面及び後面を指す。
 半導体平板において、前面及び後面の夫々にはストライプ状に配列された複数のストリップ電極が配置されている。説明の明確化及び具体化のため、互いに直交するX軸、Y軸及びZ軸を定義する。前段検出器11及び後段検出器12の夫々において、前面では、Y軸方向に伸びる複数のストリップ電極がX軸方向に沿って並べられ、後面では、X軸方向に伸びる複数のストリップ電極がY軸方向に沿って並べられている。前段検出器11の各面に配置された複数のストリップ電極において、互いに隣接するストリップ電極間はスペーシングを介して電気的に絶縁されている。後段検出器12についても同様である。
 検出器11及び12の夫々において、半導体平板の前面側の最上層がn+層に、後面側の最上層がp+層にされており、n+層とp+層とに挟まれる領域は、例えば高純度のp型半導体とされる。各検出器がゲルマニウムを利用する検出器である場合、各検出器は例えば液体窒素温度(-196℃)に冷却されている。このため、検出器11及び12の各半導体平板では、価電子帯からの電子の熱励起による伝導キャリアの供給が抑制される。尚、測定に際し冷却を必要としない半導体にて検出器11及び12が形成されることもある。
 検出器11及び12の夫々において、両面の各ストリップ電極には、半導体平板の厚みの方向に電界を生じさせるように、両面間に逆バイアス電圧が印加されている。このとき、半導体中に不純物準位からの伝導キャリアの供給が抑制された空乏層と呼ばれる領域が形成され、半導体は高抵抗状態となる。その空乏層にガンマ線が入射すると、ガンマ線と半導体中の電子とが相互作用を起こし、この電子はガンマ線からエネルギを受け取る。そして、エネルギを受け取った電子のパスに沿って価電子帯からの電子の励起によるキャリア電荷が多数生成され、生成されたキャリア電荷が逆バイアス電圧よって取り出される。ここで、キャリア電荷による電流は、ガンマ線と相互作用した電子のエネルギを反映する電流値を持ち、その電流を積分した電荷量が当該エネルギに比例する。このため、検出器11内にて上記相互作用としてコンプトン散乱が生じた場合には、コンプトン散乱によりガンマ線が失ったエネルギE1に応じた電流が検出器11にて検出される。検出器11にて検出された電流の積分値に相当する電荷量によりエネルギE1が特定される。コンプトン散乱後のガンマ線が検出器12にて全エネルギ吸収される場合には上記相互作用として光電効果が起こり、コンプトン散乱後のガンマ線の全エネルギに応じた電流が検出器12にて検出される。検出器12にて検出された電流の積分値に相当する電荷量によりエネルギE2が特定される。
 各検出器において、ガンマ線と検出器との相互作用による電流信号が強く検出される前面側及び後面側のストリップ電極の組み合わせを特定することで、相互作用位置が検出される。各イベントにおいて、前段検出器11にて検出される相互作用位置はコンプトン散乱位置R1であり、後段検出器12にて検出される相互作用位置は光電吸収位置R2であある。
 例えば、検出器11において、前面の第xA番目のストリップ電極から電流信号が検出されると共に後面の第yA番目のストリップ電極から電流信号が検出されたとき、検出器11内のブロック(xA,yA)にてガンマ線との相互作用が生じたと検出できる(検出器12についても同様)。X軸方向において、ブロック(xA,yA)の位置は前面の第xA番目のストリップ電極の位置と同じであり、Y軸方向において、ブロック(xA,yA)の位置は後面の第yA番目のストリップ電極の位置と同じである。xAは、前面に設けられたストリップ電極の本数以下の任意の自然数であり、yAは、後面に設けられたストリップ電極の本数以下の任意の自然数である。
 前段検出器11の各面に配置される複数のストリップ電極の本数は任意である(後段検出器12も同様)。各検出器において、X軸方向における相互作用位置の検出分解能は前面に配置されるストリップ電極の本数に依存し、Y軸方向における相互作用位置の検出分解能は後面に配置されるストリップ電極の本数に依存する。
 各検出器において検出される相互作用位置のZ軸成分は一定(Z軸方向における各検出器の中心位置)であって良い。但し、検出器11及び12の夫々において、前面のストリップ電極からの電流信号と後面のストリップ電極からの電流信号との時間差に基づき、Z軸方向における相互作用位置(相互作用位置のZ軸成分)を高い精度で導出するようにしても良い。
 尚、検出器11及び12を形成する半導体は、ガンマ線に対して感度を示す限り、ゲルマニウム以外の半導体でも良い。例えば、検出器11及び12を形成する半導体として、シリコン、テルル化カドミウム、テルル化カドミウム亜鉛又はダイアモンドを利用しても良い。検出器11及び12が半導体で形成される場合について説明したが、検出器11及び12の夫々は、自身とガンマ線との相互作用位置及びガンマ線のエネルギを測定可能な放射線検出器であれば、どのような構成を有していても良い。例えば、シンチレーション検出器、又は、ガス若しくは液体を用いたTPC(time projection chamber)などによって、検出器11及び12の夫々を形成しても良い。
[ガンマ線源分布画像化装置の全体構成及び基本動作]
 図3は、本実施形態に係るガンマ線源分布画像化装置1の概略的な全体構成図である。画像化装置1は、コンプトンカメラ10に加え、符号20、30及び41~44によって参照される各部位を備える。但し、コンプトンカメラ10は画像化装置1の構成要素に含まれない、と考えることも可能である。
 前段検出器11は、自身の内部でガンマ線との相互作用が生じたとき、その相互作用の内容に応じた信号11OUTを出力する。後段検出器12は、自身の内部でガンマ線との相互作用が生じたとき、その相互作用の内容に応じた信号12OUTを出力する。信号11OUTによって、検出器11内におけるガンマ線との相互作用位置(例えば位置R1)及び検出器11内での相互作用にてガンマ線から検出器11に与えられた相互作用エネルギ(例えばエネルギE1)が特定される。信号12OUTによって、検出器12内におけるガンマ線との相互作用位置(例えば位置R2)及び検出器12内での相互作用にてガンマ線から検出器12に与えられた相互作用エネルギ(例えばエネルギE2)が特定される。
 イベント測定データ取得部20は、同時事象抽出部21及び全エネルギ吸収判定部22を有し、それらを用いて、各イベントにおける検出器11及び12とガンマ線との相互作用データであるイベント測定データを生成及び取得する。
 図4は、取得部20によるイベント測定データを取得するためのフローチャートである。ステップS11及びS12において、同時事象抽出部21は、前段検出器11の出力信号11OUT及び後段検出器12の出力信号12OUTの内、検出器11及び12にて同時に生じた事象に対応する出力信号の組を同時事象検出信号として抽出する。
 例えば、同時事象抽出部21は、前段検出器11から信号11OUTが出力されるたびに、その出力時刻を示すタイムスタンプを信号11OUTに付与した上で信号11OUTを記憶部44に記憶させておくと共に、後段検出器12から信号12OUTが出力されるたびに、その出力時刻を示すタイムスタンプを信号12OUTに付与した上で信号12OUTを記憶部44に記憶させておく。そして、任意のタイミングにおいて、同時事象抽出部21は、記憶部44に記憶された任意の信号11OUTを読み出し、読み出した信号11OUTとの時間差が所定の許容時間差以内の信号12OUTを記憶部44内から探索する(ステップS11)。読み出した信号11OUTとの時間差が所定の許容時間差以内の信号12OUTが記憶部44内に存在する場合(ステップS11のY)、読み出した信号11OUTと探索された信号12OUTを同時事象検出信号として抽出する(ステップS12)。
 読み出した信号11OUTとの時間差が所定の許容時間差以内の信号12OUTが記憶部44内に存在しない場合(ステップS11のN)、ステップS15に進む。尚、第1信号との時間差が所定の許容時間差以内の第2信号とは、第1信号に対応するタイムスタンプが示す時刻との時間差が所定の許容時間差以内の時刻を示すタイムスタンプに対応付けられた第2信号を指す。
 尚、同時事象抽出部21は、記憶部44による信号11OUT及び12OUTの記憶を介することなく、検出器11及び12の出力信号11OUT及び12OUTから同時事象検出信号をリアルタイムで抽出するようにしても良い。
 ステップS12にて同時事象検出信号が抽出されるとステップS13に進む。ステップS13において、全エネルギ吸収判定部22は、抽出された同時事象検出信号に対して全エネルギ吸収判定処理を行う。全エネルギ吸収判定処理では、同時事象検出信号における信号12OUTが全エネルギ吸収を示しているかを判定する。即ち具体的には、全エネルギ吸収判定部22は、同時事象検出信号中の信号11OUT及び12OUTにより示される相互作用エネルギE1及びE2について同時事象判定式“E0-ΔE0≦E1+E2≦E0+ΔE0”が成立するか否かを判定し、同時事象判定式が成立する場合に(ステップS13のY)、同時事象検出信号における信号12OUTが全エネルギ吸収を示していると判定してステップS14に進む。エネルギ検出の誤差を考慮し、同時事象判定式にΔE0が導入されている。ΔE0は正の所定値を持つ。例えば、E0が511keV(キロエレクトロンボルト)である場合、ΔE0は5keVとされる。
 同時事象判定式が成立しない場合、全エネルギ吸収判定部22は、時事象検出信号における信号12OUTが全エネルギ吸収を示していないと判定し(ステップS13のN)、その同時事象検出信号を破棄してステップS15に進む。コンプトン散乱が生じていても同時事象判定式が成立しない場合には、測定データからコンプトン散乱角を正しく推定することができないからである。
 ステップS14において、イベント測定データ取得部20は、抽出された同時事象検出信号が有効なイベントにおける検出信号を表していると認識し、その有効なイベントにおける信号11OUT及び12OUTに基づきイベント測定データ(相互作用測定データ)を取得すると共に、取得したイベント測定データを記憶部44に記憶させる。ステップS14の後、ステップS15に進む。1つのイベント測定データは、1つのイベントにおける相互作用位置R1及びR2並びにエネルギE1及びE2を特定する情報を含んでいる。
 従って、或る1つのガンマ線が検出器11にてコンプトン散乱し、そのコンプトン散乱されたガンマ線のエネルギが検出器12にて全吸収された場合、検出器11における相互作用位置R1及び検出エネルギE1を特定する情報を含んだ信号11OUT並びに検出器12における相互作用位置R2及び検出エネルギE2を特定する情報を含んだ信号12OUTが同時事象検出信号として抽出され、抽出された同時事象検出信号から1つのイベント測定データが生成及び取得されることになる。換言すれば、取得部20は、抽出部21及び判定部22を用いてイベントを検出し、検出器11及び12の出力信号11OUT及び12OUTからイベント測定データを抽出する。
 ステップS15において、イベント測定データ取得部20は、所定の測定データ取得終了条件の成否を判定する。測定データ取得終了条件が成立する場合にはイベント測定データの取得処理を終えるが、測定データ取得終了条件が成立しない場合にはステップS11に戻る。例えば、イベント測定データの取得数が所定数(例えば数10000)に達した時点で測定データ取得終了条件が成立する。或いは例えば、記憶部44に記憶された全ての信号11OUTについて同時事象検出信号の抽出に関する処理を終えると、測定データ取得終了条件が成立する。
 図3を再び参照する。主演算部30は、CPU(Central Processing Unit)及びSIMD(single instruction multiple data)演算器等の高並列の計算アクセラレータなどから成る。主演算部30は、撮像対象200を内包する空間に所定の画像空間を設定し、複数のイベントについて生成されたイベント測定データに基づき、ガンマ線源の密度分布を示す分布画像(ガンマ線源分布画像)を画像空間内の画像として生成する。分布画像の生成方法については後述する。分布画像を含む再構成画像を生成するための専用の演算装置を用いて主演算部30を形成しても良い。専用の演算装置は、例えば、ASIC (application specific integrated circuit)やFPGA(field-programmable gate array)等を用いて構築された、画像再構成に最適なハードウェア構成を有する装置である。
 統括制御部41は、CPU(Central Processing Unit)等を用いて実現され、画像化装置1内の各部位の動作を統括的に制御する。表示部42は、液晶ディスプレイパネル等から成り、統括制御部41の制御の下、分布画像を含む任意の画像を表示する。操作部43は、ポインティングデバイス、キーボート等から成り、画像化装置1のユーザからの任意の指示及び操作を受け付ける。記憶部44は、ROM(Read Only Memory)及びRAM(Random Access Memory)から成り、主演算部30及び統括制御部41等で動作させるプログラムを記憶する他、任意のデータを記憶できる。記憶部44はHDD(Hard disk drive)やフラッシュメモリなどの二次記憶装置を含んでいても良い。
 図5に、複数のイベントにおけるイベント測定データが次々と取得される様子を示す。第iイベント(即ち、第i回目のイベント)におけるイベント測定データを、第iイベントの測定データとも呼ぶ。iは任意の自然数である。第iイベントの測定データにおける検出器11及び12での検出エネルギE1及びE2を、特に夫々Ei1、Ei2にて表す。第iイベントの測定データにおける検出器11での相互作用位置(即ち、検出されたコンプトン散乱位置)R1及び検出器12での相互作用位置(即ち、検出された光電吸収位置)R2を、特に夫々Ri1、Ri2にて表す。
 主演算部30は、イベントごとに、当該イベントにおけるイベント測定データ(Ei1、Ri1、Ei2、Ri2)に基づき、上記式(A1)を用いてコンプトン散乱角θC及びコンプトン円錐面を導出及び設定できる。
 図6に、主演算部30にて設定及び定義される画像空間ISを示す。尚、画像空間ISは、記憶部44などに内包されうる専用のメモリ装置に設定されても良い。画像空間ISは、実空間に対応付けて設定される三次元仮想空間であり、少なくとも撮像対象200の存在位置を内包している。当然、画像空間ISは有限の大きさを有する。画像空間ISの形状は任意であって良いが、ここでは、画像空間ISはXY面、YZ面及びZX面に平行な面を2面ずつ有する直方体状の空間であるとする。XY面はX軸及びY軸に平行な平面であり、YZ面はY軸及びZ軸に平行な平面であり、ZX面はZ軸及びX軸に平行な平面である。典型的には例えば、画像空間ISにおけるXY面に平行な面の1つを検出部11の前面に隣接させ、検出器11の中心、検出器12の中心及び画像空間ISの中心を一直線上に配置する。
 主演算部30は、画像空間IS内に、ガンマ線源の密度分布状態を表す分布画像を三次元画像として生成できる。図7(a)に、画像空間ISがX、Y及びZ軸方向の夫々において複数に等分割される様子を示す。画像空間ISは、X、Y及びZ軸方向の夫々において複数に等分割された要素領域から形成され、各要素領域は分布画像の画素(図7(b)参照)として機能する。分布画像を形成する各画素は、X、Y及びZ軸方向の夫々に大きさを持つ三次元画素であり、ボクセルと呼ばれうる。ここでは、分布画像を形成する各画素は立方体形状を持っているものとする。分布画像を形成する画素の内、或る注目した画素を、画素番号を示す任意の整数j又はkを用いて画素j又は画素kと表現する。また、画素jの位置を記号rjによって表す(図8参照)。図8には、第iイベントのコンプトン円錐面と画素jとコンプトン散乱位置Ri1などの関係が示される。画素jの位置rjは、画素jの中心位置(画像空間ISにおける画素jの中心座標)を表す。以下の説明において、特に記述なき限り、画素とは画像空間IS内の画素(従って分布画像内の画素)を指す。
<<分布画像の生成に関する基準方法>>
 主演算部30は、list-mode maximum-likelihood expectation-maximization(LM-ML-EM)法を元にした画像再構成法により、ガンマ線源の分布画像を生成することができる。これは、観測されたガンマ線の測定データが得られる確率が、統計的に最も高くなるように、空間中のガンマ線源分布を推定する手法である。
 LM-ML-EM法を利用した、分布画像の生成に関する基準方法を説明する。当該基準方法では、下記式(B1)に基づき分布画像が再構成される。
Figure JPOXMLDOC01-appb-M000004
 λj (l)又はλj (l+1)によって表されるλjは、分布画像の画素j(分布画像の第j番目の画素)の画素値を示す。画素jの画素値は、コンプトンカメラ10による撮像時間内に画素jにおいて崩壊した放射性核種の数の期待値を表す。放射性核種の単位時間あたりの崩壊数はポアソン分布に従うため、λjは、このポアソン分布の期待値となる。ここにおける放射性核種はガンマ線源201としてのガンマ線放出核又は陽電子放出核である。LM-ML-EM法では、反復計算によって画素値を更新することで、ガンマ線源の分布を推定される真の分布へと近づけてゆく。λj (l)は、第l回目の反復計算にて得られる画素jの画素値(画素jについて上記期待値)を示し、λj (l+1)は、第(l+1)回目の反復計算にて得られる画素jの画素値(画素jについて上記期待値)を示す。従って、λjに付随する記号“l”は反復計算の回数(従って整数)を示している。
 tijは、システム応答関数と呼ばれ、画素jから放出されたガンマ線が第iイベントとして測定される確率を表す。システム応答関数tijは、イベント測定データから推定されるコンプトン円錐面と交差する画素に対して設定される。Yiは、第iイベントにおいてコンプトンカメラ10により検出されるガンマ線の数を表す。ここでは、1つのガンマ線の検出事象を1つのイベントと捉えているため、“Yi=1”である。
 tijλjは、第iイベントが画素値λjに相当する放射性核種の集積量がある画素jから放出されたガンマ線によって測定される期待値を表す。このため、第iイベントに関し、全画素についての“tijλj”の総和は、理想的には、実際のガンマ線の検出数である“Yi=1”と一致する。よって、式(B1)の右辺における“i”の総和演算の内部では、計算によって得られたガンマ線検出数(分母の“Σkikλk”に相当)に対する実測のガンマ線検出数(分子の“Yi”に相当)の比をとっている。反復計算による画素値の更新は、この比に基づいて行われる。計算によるガンマ線検出数(分母)が実測によるガンマ線検出数(分子)に近づいてゆくように、反復計算の中で画素値が更新される。即ち例えば、計算によるガンマ線検出数が実測によるガンマ線検出数より大きい場合には、計算によるガンマ線検出数が更新前よりも小さくなるように画素値が更新され、計算によるガンマ線検出数が実測によるガンマ線検出数より小さい場合には、計算によるガンマ線検出数が更新前よりも大きくなるように画素値が更新される。
 sjは画素jに対する検出感度パラメータであり、sjによる除算は検出感度を考慮した画素値の補正を意味する。検出感度パラメータsjは、主に、画素jから見た検出器11の幾何学的効率、ガンマ線と検出器11及び12との相互作用確率、並びに、ガンマ線が放出されてから検出されるまでの減衰等の因子から構成される。これらの因子を用いた数値的演算を介して検出感度パラメータsjを設定できるほか、モンテカルロ・シミュレーションによって検出感度パラメータsjを得ても良い。
 画素jから見た検出器11の幾何学的効率は、画素jから見た検出器11の大きさ(立体角)を表す。画素jから放出されたガンマ線が検出器11に当たる確率は、画素jから見た検出器11の幾何学的効率に依存する。画素jと検出器11との距離が近いほど画素jから見た検出器11の大きさは大きくなり、画素jから放出されたガンマ線が検出器11に当たる確率が高まる。画素j及び検出器11間の距離が相応に離れている場合において、画素jから見た検出器11の幾何学的効率は、画素j及び検出器11間の距離の逆二乗にほぼ比例する。
 式(B1)中の総和(Σ)の意義を補足説明しておく。式(B1a)及び式(B1b)は、式(B1)の一部を抜粋したものである。式(B1a)において、i及びkは夫々イベント番号及び画素番号を示す。式(B1a)は、第iイベントにおける全画素についての“tikλk (l)”の総和を表す。式(B1b)の分子において、i及びjは夫々イベント番号及び画素番号を示す。式(B1b)は、式(B1a)を分母とし且つ“Yiij”を分子として持つ分数の全イベントについての総和を表す。全イベントとは、分布画像の各画素の画素値を求めるために用いられるイベントの全てを指す。
Figure JPOXMLDOC01-appb-M000005
[散乱角不確定性について]
 各イベントにおいて、前段検出器11に向かうガンマ線の到来方向はコンプトン円錐面の母線方向に限定される。しかし、コンプトン散乱角はコンプトン散乱時のエネルギによって求められるため、ガンマ線検出時の物理現象に由来する不確定性を持つ。この不確定性を示す、コンプトン散乱角の分布関数(コンプトン散乱角の誤差分布の分布関数)を散乱角分布関数と呼ぶ。コンプトン円錐面も散乱角分布関数に従った不確定性を有しており、ガンマ線源は散乱角分布関数に由来する広がりの中に確率的に存在すると推定される。
 本実施形態では、下記式(C1)で表される関数fO(θ)を散乱角分布関数として定義する。図9に、散乱角分布関数fO(θ)のグラフを示す。
Figure JPOXMLDOC01-appb-M000006
 式(C1)において、Zは、検出器11及び12を構成する半導体物質であってガンマ線と相互作用を起こす半導体物質(以下、検出器物質という)の原子番号を表し、nlは、その検出器物質中の第l番目の電子殻内の電子数を表す。式(C1)における総和演算記号“Σ”は、検出器物質中の全ての電子殻に対する総和を示す。σeは、検出器11がエネルギE1を検出する際のエネルギ検出分解能より求められる標準偏差を表す。σd,lは、検出器物質中の第l番目の電子殻の電子運動量に起因するドップラーブロードニングに起因した検出エネルギの標準偏差を表す。式(C1)の右辺の各パラメータの内、θC、θ、σe及びσd,l以外のパラメータは、画像化装置1にとって既知の所定値を持つ。
 図10を参照する。θCは上記式(A1)により求められるコンプトン散乱角を表す。θは、注目位置NPから放出されたガンマ線が検出器11内の相互作用位置R1にてコンプトン散乱したと仮定した場合におけるコンプトン散乱角を表す。即ち、検出器11及び12の相互作用位置R1及びR2を結ぶ直線(コンプトン円錐面の中心軸)と、注目位置NP及び位置R1間を結ぶ直線との成す角度が、角度θに相当する。角度(θ-θC)は、注目位置NPから放出されたガンマ線が検出器11内の相互作用位置R1にてコンプトン散乱した場合におけるコンプトン散乱角θと、上記式(A1)により求められるコンプトン散乱角θCとの誤差角度Δθを表している。
 散乱角分布関数fO(θ)の強度(即ち、関数fO(θ)の値)は、図9に示す如く、誤差角度Δθがゼロであるときに最大となり、誤差角度Δθの絶対値の増大に伴ってゼロに向かって減少してゆく。
 分布画像の再構成では、ガンマ線の測定データから推定される円錐状のガンマ線到来方向を画像空間ISに投影する逆投影演算を行う。但し、ガンマ線の測定データから推定される円錐状のガンマ線到来方向には上述の不確定性に基づく誤差が含まれるため、逆投影演算では、誤差による厚みを持った円錐を投影する必要がある。逆投影で考慮すべき円錐の厚みはコンプトン散乱角の不確定性によるものであるため、当該厚みを散乱角分布関数fO(θ)から定めれば良い。
 厚みを有するコンプトン円錐面をコンプトン円錐面領域と呼ぶ。主演算部30は、以下のように、イベントごとにコンプトン円錐面に厚みを持たせることでコンプトン円錐面領域を設定できる。
 “-180°<θ-θC≦180°”の範囲内で関数fO(θ)を角度θにて積分したときの積分値、“-θREF≦θ-θC≦θREF”の範囲内で関数fO(θ)を角度θにて積分したときの積分値を、夫々、INTA、INTBにて参照する。このとき、比“INTB/INTA”が1に近く且つ1未満の正の所定値(例えば95%)となる正の角度θREFを求める。
 図11に、或る1つのイベントに関して考慮される3つの円錐面311~313を示す。円錐面311~313は、全て、実測されたコンプトン散乱位置R1を頂点とし、且つ、コンプトン散乱後のガンマ線の飛行経路の直線上に(即ち、位置R1及びR2を通る直線上)中心軸を持った円錐面である。但し、円錐面312は角度θCを半頂角(頂角の半分;開角)とするコンプトン円錐面である一方、円錐面311、313は、夫々、角度(θC-θREF)を半頂角とする下限円錐面、角度(θC+θREF)を半頂角とする上限円錐面である。
 主演算部30は、下限円錐面311と上限円錐面313との間の領域をコンプトン円錐面領域として設定する。コンプトン円錐面312は、下限円錐面311と上限円錐面313との間の領域、即ちコンプトン円錐面領域に内包されることになる。(θREF×2)は上記厚みを角度で表現したものに相当する。位置R1からの距離がdの部分におけるコンプトン円錐面の厚みは、半径がdであって且つ中心角が(θREF×2)の弧の長さに相当する。
<<分布画像の生成に関する第1改良方法>>
 分布画像の生成に関する第1改良方法を説明する。第1改良方法並びに後述の各改良方法も、上述の基準方法と同じくLM-ML-EM法を元にしている。但し、第1改良方法では、下記式(D1)に基づき分布画像を生成する。
Figure JPOXMLDOC01-appb-M000007
 λj (l)、λj (l+1)及びtijの意義については上述した通りである。従って、tijλjは、第iイベントが画素値λjに相当する放射性核種の集積量(換言すれば存在量)がある画素jから放出されたガンマ線によって測定される期待値を表す。vij及びsijについては後述される。
 主演算部30は、第iイベントの測定データに基づき上記式(A1)に従ってコンプトン散乱角θCを算出し、コンプトン散乱角θCを半頂角とする第iイベントのコンプトン円錐面と散乱角分布関数とに基づき第iイベントのコンプトン円錐面領域を設定する。そして、主演算部30は、第iイベントのコンプトン円錐面領域と交差する画素を抽出し、抽出した各画素jに対してゼロより大きなtijを設定する。第iイベントのコンプトン円錐面領域と交差しない各画素jに対しては“tij=0”として良い。このようなシステム応答関数tijの設定処理はイベントごとに行われる。各イベントの測定データに基づくシステム応答関数tijの設定方法として公知の任意の方法を利用して良い。
 式(D1)中の総和(Σ)の意義を補足説明しておく。式(D1a)及び式(D1b)は、式(D1)の一部を抜粋したものである。式(D1a)において、i及びkは夫々イベント番号及び画素番号を示す。式(D1a)は、第iイベントにおける全画素についての“tikλk (l)”の総和を表す。式(D1b)の分子において、i及びjは夫々イベント番号及び画素番号を示す。式(D1b)は、式(D1a)とsijとの積を分母とし且つ“vijij”を分子として持つ分数の全イベントについての総和を表す。全イベントとは、分布画像の各画素の画素値を求めるために用いられるイベントの全てを指す。
Figure JPOXMLDOC01-appb-M000008
―――vijの導入について―――
 式(B1)による基準方法では、計算によって推定されるガンマ線検出数(Σkikλk)が極めて小さい場合などにおいて、画素値が非常に大きな値に更新されるオーバーフィッティングが生じる。これは、図12に示す如く、円錐状に推定されるガンマ線到来方向が画像空間ISに対してかすめるように交差した場合などに生じ、画像の辺縁部でアーチファクト(虚像)が生じる原因をとなる。
 図13等を参照し、この現象について説明を加える。図13において、円錐360は、第1イベントの測定データに基づき導出されたコンプトン円錐面を側面として持ち且つ母線の長さをdAとする円錐であり、円錐370は、第2イベントの測定データに基づき導出されたコンプトン円錐面を側面として持ち且つ母線の長さをdAとする円錐である。円錐360の底面付近において円錐360は画像空間IS内に収まっている一方で、円錐370の底面付近において円錐370の一部は画像空間ISから飛び出している。図13での円錐370の図示において、円錐370と画像空間ISとが交差する部分(即ち、円錐370の内、画像空間IS内に位置している部分)については実線で示し、円錐370と画像空間ISとが交差しない部分については破線(破線371に相当)で示している。
 今、第1イベントにおけるコンプトン散乱位置R11を中心に持ち且つ半径dAを持つ球面(以下、第1球面と呼ぶ)と、第2イベントにおけるコンプトン散乱位置R21を中心に持ち且つ半径dAを持つ球面(以下、第2球面と呼ぶ)とを想定する。図14(a)において、380は、第1球面と画像空間ISとの交差部分である曲面をXY平面上に投影して示したものである。図14(b)において、390は、第2球面と画像空間ISとの交差部分である曲面をXY平面上に投影して示したものである。第1球面と画像空間ISとの交差部分をXY平面上に投影したものの形状は実際には四角形にならないこともあるが、図14(a)では、簡単化のため、その形状を四角形と仮定している。図14(b)についても同様である。
 第1イベントの測定データに基づくコンプトン円錐面に対し散乱角分布関数に基づく厚みが付与されることで、第1イベントのコンプトン円錐面領域が設定される。図14(a)において、斜線領域381は、第1イベントのコンプトン円錐面領域と曲面380との交差領域(即ち、第1イベントのコンプトン円錐面領域と第1球面と画像空間ISとの交差領域)を表している。図14(b)において、斜線領域391は、第2イベントのコンプトン円錐面領域と曲面390との交差領域(即ち、第2イベントのコンプトン円錐面領域と第2球面と画像空間ISとの交差領域)を表している。図14(b)において、ドット領域392は、第2イベントのコンプトン円錐面領域と第2球面との交差領域の内、画像空間ISとは交差していない領域を表している。
 ここでは、第1イベントのコンプトン円錐面領域と第1球面との交差領域の内、画像空間ISに交差している領域の割合は100%であるとする。故に、領域381に属する画素が100個であって且つ第1イベントにおけるガンマ線源がコンプトン散乱位置R11から距離dAだけ離れた位置にあり且つ領域381内の散乱角分布関数の強度が一様であると仮定したならば、領域381に属する各画素に第1イベントのガンマ線源が存在する確率は均等に概ね“1/100”であると推定され、そのような推定に基づくシステム応答関数t1jが設定される。
 一方、第2イベントのコンプトン円錐面領域と第2球面との交差領域の内、画像空間ISに交差している領域の割合は50%であるとする。故に、領域391に属する画素が50個であって且つ第2イベントにおけるガンマ線源がコンプトン散乱位置R21から距離dAだけ離れた位置にあり且つ領域391内の散乱角分布関数の強度が一様であると仮定したならば、領域391に属する各画素に第2イベントのガンマ線源が存在する確率は均等に概ね“1/50”であると推定され、そのような推定に基づくシステム応答関数t2jが設定される。
 しかし、領域391及び392の全てが画像空間IS内に存在するとしたならば、各イベントにおいて概ね100個の画素の何れかからガンマ線が到来したと推定されるはずであり、従って、領域391に属する各画素に第2イベントのガンマ線源が存在する確率は均等に概ね“1/100”であると推定されるべきである。但し実際には、コンプトン円錐面領域と第2球面との交差領域の一部(領域392)が画像空間ISの外に位置しているため、領域391に属する各画素に第2イベントのガンマ線源が存在する確率は概ね“1/50”であると推定される。つまり、領域391に対しては、ガンマ線源がそこに存在すると過剰に強く推定されやすい。このような不正確な推定の実行は、コンプトン円錐面領域と第2球面との交差領域の一部(領域392)が画像空間ISの外に位置していることに由来する。
 このような傾向は、1つのイベントにおける複数の画素間でも発生しうる。図15において、円錐410は、第iイベントの測定データに基づき導出されたコンプトン円錐面を側面として持ち且つ母線の長さをdAとする円錐であり、円錐420は、第iイベントの測定データに基づき導出されたコンプトン円錐面を側面として持ち且つ母線の長さをdBとする円錐である(0<dA<dB)。円錐410及び420に対応するイベントは共通であり、また“0<dA<dB”であるため、円錐410の円錐面は円錐420の円錐面の一部である。円錐410の底面付近において円錐410は画像空間IS内に収まっている一方で、円錐420の底面付近において円錐420の一部は画像空間ISから飛び出している。図15での円錐420の図示において、円錐420と画像空間ISとが交差する部分(即ち、円錐420の内、画像空間IS内に位置している部分)については実線で示し、円錐420と画像空間ISとが交差しない部分については破線(破線421に相当)で示している。
 今、第iイベントにおけるコンプトン散乱位置Ri1を中心に持ち且つ半径dAを持つ球面(以下、球面Aと呼ぶ)と、第iイベントにおけるコンプトン散乱位置Ri1を中心に持ち且つ半径dBを持つ球面(以下、球面Bと呼ぶ)とを想定する。図16(a)において、440は、球面Aと画像空間ISとの交差部分である曲面をXY平面上に投影して示したものである。図16(b)において、450は、球面Bと画像空間ISとの交差部分である曲面をXY平面上に投影して示したものである。球面Aと画像空間ISとの交差部分をXY平面上に投影したものの形状は実際には四角形にならないこともあるが、図16(a)では、簡単化のため、その形状を四角形と仮定している。図16(b)についても同様である。
 第iイベントの測定データに基づくコンプトン円錐面に対し散乱角分布関数に基づく厚みが付与されることで、第iイベントのコンプトン円錐面領域が設定される。図16(a)において、斜線領域441は、第iイベントのコンプトン円錐面領域と曲面440との交差領域(即ち、第iイベントのコンプトン円錐面領域と球面Aと画像空間ISとの交差領域)を表している。図16(b)において、斜線領域451は、第iイベントのコンプトン円錐面領域と曲面450との交差領域(即ち、第iイベントのコンプトン円錐面領域と球面Bと画像空間ISとの交差領域)を表している。図16(b)において、ドット領域452は、第iイベントのコンプトン円錐面領域と球面Bとの交差領域の内、画像空間ISとは交差していない領域を表している。
 ここでは、第iイベントのコンプトン円錐面領域と球面Aとの交差領域の内、画像空間ISに交差している領域の割合は100%であるとする。故に、領域441に属する画素がNA個であって且つ第iイベントにおけるガンマ線源がコンプトン散乱位置Ri1から距離dAだけ離れた位置にあり且つ領域441内の散乱角分布関数の強度が一様であると仮定したならば、領域441に属する各画素に第iイベントのガンマ線源が存在する確率は均等に概ね“1/NA”であると推定され、そのような推定に基づくシステム応答関数tijが設定される。
 一方、第iイベントのコンプトン円錐面領域と球面Bとの交差領域の内、画像空間ISにも交差している領域の割合は50%であるとする。故に、領域451に属する画素がNB個であって且つ第iイベントにおけるガンマ線源がコンプトン散乱位置Ri1から距離dBだけ離れた位置にあり且つ領域451内の散乱角分布関数の強度が一様であると仮定したならば、領域451に属する各画素に第iイベントのガンマ線源が存在する確率は均等に概ね“1/NB”であると推定され、そのような推定に基づくシステム応答関数tijが設定される。
 しかし、領域451及び452の全てが画像空間IS内に存在するとしたならば概ね(2×NB)個の画素の何れかからガンマ線が到来したと推定されるはずであり、従って、領域451に属する各画素に第iイベントのガンマ線源が存在する確率は均等に概ね“1/(2×NB)”であると推定されるべきである。但し実際には、コンプトン円錐面領域と球面Bとの交差領域の一部(領域452)が画像空間ISの外に位置しているため、領域451に属する各画素に第iイベントのガンマ線源が存在する確率は概ね“1/NB”であると推定される。つまり、領域451に対しては、ガンマ線源がそこに存在すると過剰に強く推定されやすい。このような不正確な推定の実行は、コンプトン円錐面領域と球面Bとの交差領域の一部(領域452)が画像空間ISの外に位置していることに由来する。
 第1改良方法にて導入される確率パラメータvij(上記式(D1)参照)、このような現象(オーバーフィッティング)を抑制するための正則化因子として働く。vijは、下記式(D2)にて定義される。
Figure JPOXMLDOC01-appb-M000009
 上述の説明から明らかなように、Ri1は、第iイベントにて検出されたガンマ線と前段検出器11との相互作用位置(即ちコンプトン散乱位置)であり(図8参照)、第iイベントにおけるコンプトン円錐面の頂点位置を表す。rjは画素jの位置(中心位置)を示す。rは画像空間IS内の任意の位置(座標)を表す。Vkは画素kの体積を表す。“r∈Vk”は画素kに内包される位置r(換言すれば画素kの体積に包含される座標r)を表す。故に、vijは、第iイベントのコンプトン散乱位置Ri1を中心とし且つ画素jの位置(rj)及びコンプトン散乱位置Ri1間の距離を半径とする球面と交差する全画素kについてのpikの総和を表す。
 確率パラメータvijは、第iイベントでコンプトン散乱を起こしたガンマ線が画像空間IS内から到来した確率を示す。但し、確率パラメータは画素ごとに設定され、vijは画素jに対して設定される確率パラメータである。より詳細には、vijは、第iイベントでコンプトン散乱を起こしたガンマ線がコンプトン散乱位置Ri1を中心とし且つ画素jの位置(rj)及びコンプトン散乱位置Ri1間の距離を半径とする球面上の位置から到来したと仮定した場合に、当該ガンマ線が画像空間IS内から到来した確率を示す。
 コンプトン散乱角の不確定性により第iイベントについて推定されるガンマ線到来方向は幅を持った方向となるが、pikは、第iイベントについて推定されるガンマ線到来方向の内、画素kを通るガンマ線到来方向が占める割合を表す。つまり、pikは、第iイベントのガンマ線到来方向を確率変数として捉えたときの、画素kを通るガンマ線到来方向の確率密度を表す。式(D2)に示す如く、vijは、“|Ri1-rj|=|Ri1-r|”を満たす位置rを内包する全ての画素kについてのpikの総和であるため、pikはvijに対する画素kの寄与分を表す。
 pikは、コンプトン散乱角の不確定性に基づくものであり、式(D3)によって定義される。尚、式(D3)では、pikにおける変数kの代わりに変数jを用いてpijを定義している。式(D3)中の角度θi(r)は式(D3a)により定義される。
Figure JPOXMLDOC01-appb-M000010
 fO,iは、第iイベントにおける散乱角分布関数(コンプトン散乱角の誤差分布の分布関数)を表す。式(D3)では、fO,iを、位置rに依存する角度θi(r)の関数として捉えている。位置rは上述したように画像空間IS内の任意の位置を表し、図10の注目位置NPに相当すると考えて良い。θi(r)は、第iイベントが位置rから放出されたガンマ線によって生じたと仮定した場合における、ガンマ線のコンプトン散乱角を表す。Rii及びRi2は、夫々、第iイベントにて検出されたガンマ線と検出器11の相互作用位置(即ち検出されたコンプトン散乱位置)及び第iイベントにて検出されたガンマ線と検出器12の相互作用位置(即ち検出された光電吸収位置)である。Vjは画素jの体積を表す。式(D3)では、位置rを積分変数として用い、画素jに属する位置rに対して“fO,ii(r)]/{|r-Ri122πsin[θi(r)]}”を積分することにより、pijを得る。
 fO,ii(r)]としての散乱角分布関数fO(θ)を考える場合、上記式(C1)中のθCは、第iイベントに関し式(A1)を用いて算出されたコンプトン散乱角である。fO,ii(r)]としての散乱角分布関数fO(θ)を考える場合、上記式(C1)中のθは、θi(r)を表し、位置Ri1及びRi2を結ぶ直線(コンプトン円錐面の中心軸)と位置Ri1及びrを結ぶ直線との成す角度に相当する。
 fO,ii(r)]の値は、位置rにおけるコンプトン散乱角の誤差分布の確率密度を表しており、これを画素jの体積Vjにわたって積分することでpijが求められる。しかし、単純な積分では、画素jから見た検出器11の幾何学的効率が考慮されない。また、fO,iは角度の関数であるため、角度の違いによる立体角を考慮する必要がある。そこで、式(D3)では、fO,ii(r)]を、幾何学的効率の近似“|r-Ri1-2”の逆数及び角度θi(r)の微小立体角“2πsin[θi(r)]”で除したものを積分している(図17参照)。
 尚、各イベントにおいて、コンプトン円錐面領域に属さない画素に対する関数fO,iの値は十分に小さいと考えられる。従って、式(D2)の右辺は、“|Ri1-rj|=|Ri1-r|”を満たす位置rを内包し且つ第iイベントのコンプトン円錐面領域に属する全ての画素kについてのpikの総和である、と考えて良い。位置Ri1を中心とし且つ距離|Ri1-rj|を半径とする球面が、画素k内の任意の位置を通過して画素kと交差するとき、当該画素kについて“|Ri1-rj|=|Ri1-r|”が満たされる。
 このように、第1改良方法に係る主演算部30は、各イベントでコンプトン散乱したガンマ線が画像空間IS内から到来した確率を示す確率パラメータvijを、各イベントでの測定データ(Ri1、Ri2、Ei1、Ei2)を元にした演算により分布画像の画素ごとに設定し、画素ごとに設定した確率パラメータvijを用いて分布画像を生成する(即ち各画素の画素値を求める)。
 各イベントでの測定データを元にして確率パラメータvijが設定されるべき画素(対象画素)は分布画像の全画素である必要はない(但し、分布画像の全画素であっても構わない)。各イベントにおいて、コンプトン円錐面領域内の各画素を確率パラメータvijが設定されるべき対象画素に含めておけば足る。即ち、各イベントにおいて、コンプトン円錐面領域内の各画素に対してのみ、式(D2)に基づきvijを算出すれば足る。各イベントにおいて、コンプトン円錐面領域外の各画素に対しては“vij=0”としておいても良い。
 これらを考慮すれば、第1改良方法に係る主演算部30は、確率パラメータvijを、各イベントでの測定データを元にした演算により分布画像内の複数の対象画素の夫々に対して個別に設定し、設定した確率パラメータvijを用いて分布画像を生成する(即ち各画素の画素値を求める)、とも言える。主演算部30は、イベントごとに且つ対象画素ごとに確率パラメータvijを設定することになる。つまり、主演算部30は、複数のイベントと複数の対象画素の組み合わせの夫々に注目して、組み合わせごとに確率パラメータvijを設定することになる。
 第iイベント及び1つの画素jについて注目した場合(第iイベント及び画素jの組み合わせに注目した場合)、主演算部30は、検出されたコンプトン散乱位置Ri1を中心とし且つ注目画素jの位置(rj)及びコンプトン散乱位置Ri1間の距離を半径(例えば図15のdA又はdB)とする球面と、注目された第iイベントのコンプトン円錐面領域と、画像空間ISとの交差状態に基づき、注目された第iイベントの注目画素jに対する確率パラメータvijを設定する。例えば、図16(a)の斜線領域441又は図16(b)の斜線領域451が、上記球面とコンプトン円錐面領域と画像空間ISの交差領域に相当し、図16(a)及び(b)の状況間では、それらの交差状態が互いに異なる。
 より具体的には、第iイベント及び1つの画素jについて注目した場合、主演算部30は、検出されたコンプトン散乱位置Ri1を中心とし且つ注目画素jの位置(rj)及びコンプトン散乱位置Ri1間の距離を半径(例えば図15のdA又はdB)とする球面と、注目された第iイベントのコンプトン円錐面領域と、画像空間ISとの交差領域(例えば図16(a)の斜線領域441又は図16(b)の斜線領域451)を求める。次に、その交差領域に配置された各画素を抽出し、抽出した各画素に対して散乱角分布関数(fO,i)の強度に応じた指標(pij)を導出する(式(D3)参照)。そして、抽出した各画素について導出された指標の総和に基づき(式(D2)参照)、注目された第iイベントの注目画素jに対する確率パラメータvijを設定する。
 vijは積分範囲を全散乱角の範囲とした散乱角不確定性の積分と考えられるため、vijの理想的な値は1となる。しかし、コンプトン円錐面領域の一部が画像空間ISの外にある場合には、vijは1より小さくなる。従って、コンプトン円錐面領域が画像空間ISと僅かしか交差しなかったイベントは、vijの導入によって再構成画像(分布画像)に対する寄与が低減することになり、アーチファクトが低減される。更に、図15並びに図16(a)及び(b)を参照して説明したように、1つのイベントにおいても注目画素が異なれば上記交差状態が異なってくる。第1改良方法では、これを考慮し、各イベントにおいて画素ごとに最適なvijを設定するためアーチファクトの低減効果が高まる。
―――sijの導入について―――
 検出感度を構成する因子で最も大きな影響を持つのは、個々の画素jから見た検出器11の幾何学的効率であり、その幾何学的効率は、近似的には画素j及び検出器11間の距離の逆二乗に比例する。しかし、式(B1)による基準方法では、検出感度(sj)が画素jの位置に対してのみ依存すると仮定しているため、コンプトン散乱位置(検出器11での相互作用位置)の違いによる幾何学的効率の違いが平均化される。
 即ち、図18を参照し、コンプトン散乱位置が位置510であるときと位置520であるときとで、コンプトン散乱位置及び画素j間の距離(511、521)は相違し、結果、幾何学的効率が異なるはずであるのに、検出感度(sj)が画素jの位置に対してのみ依存すると仮定したならば、コンプトン散乱位置の違いによる幾何学的効率の違いが平均化されることになる。検出器11に近い画素ほど、この平均化による影響が大きく(距離511及び521の違いが出やすいため)、検出感度パラメータを用いた感度補正の正確性が低下する。
 そこで、第1改良方法では、イベントごとのコンプトン散乱位置Ri1をも考慮して画素ごとの検出感度パラメータを設定する。上記式(D1)におけるsijは、第iイベントのコンプトン散乱に対する検出器11の検出感度を示す検出感度パラメータである。但し、当該検出感度パラメータは画素ごとに設定され、sijは画素jに対して設定される検出感度パラメータである。上記式(D1)において、sjjによる除算は検出感度を考慮した画素値の補正を意味する。図19(a)及び(b)に、基準方法と第1改良方法との比較図を示す。検出感度パラメータsijは下記式(D4)によって定義される。
Figure JPOXMLDOC01-appb-M000011
 ここで、σt(E0)は、E0のエネルギを持ったガンマ線と検出器物質との全相互作用断面積を表し、di,j1は、第iイベントにおける画素jからのガンマ線の検出器11内の飛程を表す(図20参照)。即ち、di,j1は、第iイベントにおける画素jからのガンマ線がコンプトン散乱を起こすまでに検出器11内を通過した距離を表す。そのため、“exp[-σt(E0)di,j1]”は、第iイベントにおいて画素jからのガンマ線がコンプトン散乱位置Ri1まで検出器11内の検出器物質と相互作用せずに到達する確率(換言すれば、第iイベントにおける画素jからのガンマ線がコンプトン散乱するまでに検出器11内で生じる減衰)を表す。σt(E0)の値は、ガンマ線源201が放出するガンマ線の初期エネルギE0と検出器物質の物性とに基づき一意に定まり、画像化装置1にとって既知である。飛程di,j1は、画素jの位置rj及び第iイベントのコンプトン散乱位置Ri1と、検出器11の形状とから定まる。
 ここでは、図21(a)及び(b)に示す如く、検出器11は、相互作用位置をX軸、Y軸、Z軸方向において、夫々、m、n及びo段階で検出できるものとする(m、n及びoは2以上の任意の整数)。そうすると、検出器11は(m×n×o)個の独立した検出器要素の集合体であると考えることができる。このように考えると、各々の検出器要素が画像空間IS内の全画素に対して検出感度を持つことになる。Ri1は、第iイベントでのコンプトン散乱位置を内包する検出器要素の中心位置を表すと考えて良い。そうすると、“|Ri1-rj-2”は、画素jから見た、位置Ri1を中心に持つ検出器要素の幾何学的効率の近似を表す。
 このように、第1改良方法に係る主演算部30は、各イベントでのコンプトン散乱の検出感度を示す検出感度パラメータsijを各イベントでの測定データ(Ri1を含む)に元にした演算により分布画像の画素ごとに設定し、画素ごとに設定した検出感度パラメータsijを用いて分布画像を生成する(即ち各画素の画素値を求める)。
 各イベントでの測定データを元にして検出感度パラメータsijが設定されるべき画素(対象画素)は分布画像の全画素である必要はない(但し、分布画像の全画素であっても構わない)。各イベントにおいて、コンプトン円錐面領域内の各画素を検出感度パラメータsijが設定されるべき対象画素に含めておけば足る。即ち、各イベントにおいて、コンプトン円錐面領域内の各画素に対してのみ、式(D4)に基づきsijを算出すれば足る。各イベントにおいて、コンプトン円錐面領域外の各画素に対するsijには一定値を代入しておいて良い。
 これらを考慮すれば、第1改良方法に係る主演算部30は、検出感度パラメータsijを、各イベントでの測定データを元にした演算により分布画像内の複数の対象画素の夫々に対して個別に設定し、設定した検出感度パラメータsijを用いて分布画像を生成する(即ち各画素の画素値を求める)、とも言える。主演算部30は、イベントごとに且つ対象画素ごとに検出感度パラメータsijを設定することになる。つまり、主演算部30は、複数のイベントと複数の対象画素の組み合わせの夫々に注目して、組み合わせごとに検出感度パラメータsijを設定することになる。
 第iイベント及び1つの画素jについて注目した場合(第iイベント及び画素jの組み合わせに注目した場合)、主演算部30は、検出されたコンプトン散乱位置Ri1と注目画素jの位置rjに応じて検出感度パラメータ(sij)を設定することになる。より具体的には、第iイベント及び1つの画素jについて注目した場合、主演算部30は、コンプトン散乱位置Ri1と注目画素jの位置rjに基づく検出器11内におけるコンプトン散乱前のガンマ線の飛程di,j1と、コンプトン散乱位置Ri1及び注目画素jの位置rj間の距離|Ri1-rj|とに基づき、注目された第iイベントの注目画素jに対する検出感度パラメータsijを設定する。
 例えば、ガンマ線源と検出器11との距離が比較的小さいとき、当該ガンマ線源からのガンマ線は検出器11に当たりやすいので検出器11で検出されやすく、ガンマ線源と検出器11との距離が比較的大きいとき、当該ガンマ線源からのガンマ線は検出器11に当たりにくいので検出器11で検出されにくい。仮に、このような検出されやすさ/検出されにくさを考慮せずに分布画像を生成すると、検出器11からの距離が遠いガンマ線源の存在が分布画像に反映されにくくなる。検出感度パラメータの導入によって、検出されやすさ/検出されにくさが考慮された感度補正が成される。この際、各画素の位置だけでなく検出器11内での相互作用位置Ri1をも考慮して検出感度パラメータ(sij)を設定することにより、相互作用位置Ri1を考慮せずに検出感度パラメータを設定する方法(図19(a))と比べて、より妥当な感度補正を行うことが可能となる。結果、分布画像がガンマ線源の分布をより正確に表すようになる。また、第1改良方法では、相互作用位置Ri1を用いて検出器11内でのガンマ線減衰(exp[-σt(E0)di,j1])をも考慮するため、分布画像の妥当性が増す(ガンマ線源分布の推定精度が増す)。
―――tijについて―――
 システム応答関数tijは下記式(D5)に従って決定される。
Figure JPOXMLDOC01-appb-M000012
 式(D5)における“exp[-σt(E0)di,j1]”は、上記式(D4)に示したものと同じものであり、第iイベントにおいて画素jからのガンマ線がコンプトン散乱位置Ri1まで検出器11内の検出器物質と相互作用せずに到達する確率を表す。式(D5)における“dσc,ij/dΩ”は、画素jから放出されたガンマ線が第iイベントを形成するガンマ線として検出された場合におけるコンプトン散乱微分断面積を表す。即ち、“dσc,ij/dΩ”は、第iイベントにおいて画素jからのガンマ線がコンプトン散乱を生じさせる確率を表す。この確率はコンプトン散乱角によって変化するため、“dσc,ij/dΩ”は画素jに依存する。
 E0’はコンプトン散乱後のガンマ線のエネルギを表す。故に、式(D5)におけるE0’は、ガンマ線の初期エネルギE0から第iイベントにおける検出エネルギE1を差し引いたものに相当する。σt(E0’)は、E0’のエネルギを持ったガンマ線と検出器物質との全相互作用断面積を表し、di,12は、第iイベントにおけるコンプトン散乱後のガンマ線が検出器12にて光電吸収されるまでの、ガンマ線の検出器11及び12内の飛程を表す。即ち、di,12は、第iイベントにおけるコンプトン散乱後のガンマ線が検出器12にて光電吸収されるまでに検出器11及び12内を通過した距離を表す。そのため、“exp[-σt(E0’)di,12]”は、第iイベントにおけるコンプトン散乱後のガンマ線が光電吸収位置Ri2まで検出器11及び12内の検出器物質と相互作用せずに到達する確率(換言すれば、第iイベントにおけるコンプトン散乱後のガンマ線が位置Ri2にて検出器12に光電吸収されるまでに検出器11及び12内で生じる減衰)を表す。
 |Ri2-Ri1-2は、コンプトン散乱位置Ri1から見た光電吸収位置Ri2の幾何学的効率(立体角)の近似を表している。σp(E0’)は、E0’のエネルギを持ったガンマ線と検出器物質との光電吸収断面積を表す、即ちコンプトン散乱後のガンマ線が光電吸収を生じさせる確率を表す。
 主演算部30は、イベントごとの測定データ(Ri1、Ri2、Ei1、Ei2)に基づきイベントごとに且つ画素ごとにシステム応答関数tijを導出及び設定することができる。式(D5)の右辺中の物理量の内、第iイベントの測定データに依存しない物理量の値は画像化装置1にとって既知である。
[分布画像の生成動作フローチャート]
 図22は、各イベントの測定データに基づく分布画像の生成動作フローチャートである。図22を参照して、各イベントの測定データに基づく分布画像の生成動作手順を説明する。図4による各イベントの測定データの取得後、ステップS21に至る。ステップS21において、主演算部30は、分布画像の全画素の画素値λk (0)に対し所定の共通初期値(例えば1)を代入すると共に、反復計算の回数を示す変数lにゼロを代入する。主演算部30は、続くステップS22にて変数iに1を代入した後、ステップS23~S28の処理を順次実行する。
 具体的には、主演算部30は、ステップS23にて第iイベントの測定データ(Ri1、Ri2、Ei1、Ei2)を記憶部44から読み出し、ステップS24において、第iイベントの測定データに基づき上記式(A1)に従ってコンプトン散乱角θCを算出すると共にコンプトン円錐面を設定する。続くステップS25において、主演算部30は、散乱角分布関数fO(θ)に基づき第iイベントのコンプトン円錐面に厚みを持たせることで第iイベントのコンプトン円錐面領域を設定する。その後、ステップS26において、主演算部30は、第iイベントのコンプトン円錐面領域内の各画素を対象画素として捉え、対象画素ごとに、上述の各式等に従いつつ測定データに基づいてvij、sij及びtijを算出及び設定する。コンプトン円錐面領域外の各画素(非対象画素)に対してはシステム応答関数tijをゼロに設定すれば良い。
 ステップS26に続くステップS27において、主演算部30は、第(i-1)イベントについてまで求めた式(D1b)の総和計算結果に、第iイベントについて求めたvij、sij及びtijに基づく計算結果を加算することで、第1~第iイベントについての式(D1b)の総和計算結果(即ち、式(D1b)の値)を求める。その後、ステップS28において、全イベントについての測定データが記憶部44から読み出されたか否かが確認される。未だ読み出されていない測定データが存在する場合には(ステップS28のN)、ステップS29にて変数iに1を加算してからステップS23に戻り、ステップS23~S28の処理が繰り返される。全イベントについての測定データが記憶部44から読み出された場合には(ステップS28のY)ステップS30に進む。
 ステップS30において、主演算部30は画素値の更新処理を行う。更新処理では、式(D1)に従い、更新された画素値である画素値λk (l+1)を求める。即ち、ステップS27にて求められた最新の式(D1b)の総和計算結果にλk (l)を乗じることで、分布画像の各画素の画素値λk (l+1)を求める。その後のステップS31において、主演算部30は、所定の更新終了条件の成否を確認する。更新終了条件が成立する場合には、ステップS33に進んで、最新の各画素値λkを持つ分布画像が記憶部44に記憶された後、図22の動作を終える。更新終了条件が成立しない場合にはステップS32にて変数lに1を加算してからステップS22に戻り、ステップS22以降の各処理が繰り返される。
 更新終了条件は、例えば、ステップS30の更新処理の実行回数が所定回数に達した時点で成立する。或いは例えば、ステップS22~S32の処理の繰り返し実行の中で、所定の終了指示操作が操作部43に入力されたとき、更新終了条件が成立しても良い。更に或いは例えば、ステップS30の更新処理による各画素値λkの変化量が充分に小さくなったとき(例えば、分布画像の全画素における画素値λk (l+1)の総和と画素値λk (l)の総和との差が所定値以下となったとき)、更新終了条件が成立しても良い。
<<分布画像の生成に関する第2改良方法>>
 分布画像の生成に関する第2改良方法を説明する。第2改良方法では、下記式(E1)に基づき分布画像が生成される。式(E1)では、基準方法を元にしてvijが導入されてはいるが、sijの導入は見送られている。sijの代わりにsjが画素値λjの算出式に組み込まれている点を除き、第2改良方法は第1改良方法と同様である。検出感度パラメータsjは、例えば下記式(E2)に従って画素ごとに求められる。
Figure JPOXMLDOC01-appb-M000013
 上述したように、σt(E0)はE0のエネルギを持ったガンマ線と検出器物質との全相互作用断面積を表し、rjは画素jの位置(中心位置)を表す。式(E2)において、pは、検出器11を形成する第1~第(m×n×o)番目の検出器要素(図21(a)及び(b)参照)中の或る検出器要素の番号を表しており、Rpは第p番目の検出器要素の中心位置を表している。dp,j1は、位置rjから位置Rpまでの直線経路の内、検出器11内に位置する経路の長さである。つまり、dp,j1は、図20のdi,j1に類似し、画素jからのガンマ線が第p番目の検出器要素に至るまでに検出器11内を通過した距離を表す。そのため、式(E2)の“exp[-σt(E0)dp,j1]”は、画素jからのガンマ線が第p番目の検出器要素まで検出器11内の検出器物質と相互作用せずに到達する確率を表す。“|Rp-rj-2”は、画素jから見た、第p番目の検出器要素の幾何学的効率の近似を表す。全検出器要素に対するガンマ線の検出感度“exp[-σt(E0)dp,j1]・|Rp-rj-2”の総和が、式(E2)のsjとなる。
 sijの定義式(D4)との比較からも理解されるように、式(E2)のsjは、全ての検出器要素に対して設定される検出感度(ガンマ線検出感度)の合計として計算される。検出器要素の個数は固定値であるので、この検出感度の合計は実質的に全検出器要素の検出感度の平均値とみなせる。各検出器要素の大きさは、検出器11の位置検出精度のX、Y及びZ軸成分の積(例えば、3×3×1[mm3])とされる。
 尚、検出器11の大きさに比べて画像空間ISが充分に小さい場合などにおいては、検出感度が一様であると考えて、全画素に対し“sj=1”と設定することも可能である。
 第2改良方法では、基準方法との比較において、vijの導入による有利な効果が奏される。但し、sijをも導入した第1改良方法の方が、第2改良方法よりも正確な分布画像の生成に適している。
<<分布画像の生成に関する第3改良方法>>
 分布画像の生成に関する第3改良方法を説明する。第3改良方法では、下記式(F1)又は式(F2)に基づき分布画像が生成される。式(F1)又は式(F2)では、基準方法を元にしてsijが導入されてはいるが、vijの導入は見送られている。vijの代わりにYi又はviが画素値λjの算出式に組み込まれている点を除き、第3改良方法は第1改良方法と同様である。
Figure JPOXMLDOC01-appb-M000014
 上述したように、Yiは全イベントにおいて“1”であって良い。確率パラメータviはイベントごとに設定され、1つの確率パラメータviは各画素に対して共通に適用される。例えば、確率パラメータviは下記式(F3)に従って算出される。式(F3)の右辺は、第iイベントのコンプトン円錐面領域内の全画素についてのpijの総和を表している。
Figure JPOXMLDOC01-appb-M000015
 第3改良方法では、基準方法との比較において、sijの導入による有利な効果が奏される。但し、vijをも導入した第1改良方法の方が、第3改良方法よりも正確な分布画像の生成に適している。vijの代わりにviを用いることでも基準方法との比較においてアーチファクトの低減が図られるが、図15並びに図16(a)及び(b)を参照した上述の説明から理解されるように、各イベントにおいて画素ごとに最適なvijを設定する第1改良方法の方がアーチファクトの低減効果が高い。
[実験結果]
 基準方法、第1改良方法及び第2改良方法を利用した実験の内容及び結果を説明する。図23は、当該実験で利用された担癌マウスの概略平面図である。3つの腫瘍組織A431、4T1及びC6が移植された担癌マウスに64Cu標識抗HER2抗体を投与したものを撮像対象200とする。64Cu標識抗HER2抗体は、64Cuを内包する抗体であって、HER2を含む腫瘍に集積する抗体である。3つの腫瘍組織の内、腫瘍組織A431においてHER2が最も高発現であり、故に腫瘍組織A431に対して高い抗体の集積が期待される。本実験では、担癌マウス中の64Cuがガンマ線源201として機能し、64Cuより放出される陽電子の対消滅に伴う511keV(キロエレクトロンボルト)のガンマ線が、ガンマ線源201からのガンマ線として機能する。当該実験において、コンプトンカメラ10を用いた撮像時間は9時間であり、測定されたイベント数は1.9×106である。
 基準方法、第2改良方法、第1改良方法への実験の結果として得られる三次元分布画像を、夫々、三次元分布画像610、620、630と呼ぶ(図24参照)。当該実験における反復計算の反復回数は一律に60回である。つまり、当該実験の条件下で、三次元分布画像610、620、630は、夫々、式(B1)、式(E1)、式(D1)におけるλj (60)を画素jの画素値として有する三次元分布画像である。
 図25に、三次元分布画像610に基づく二次元のMIP画像611~613、三次元分布画像620に基づく二次元のMIP画像621~623、及び、三次元分布画像630に基づく二次元のMIP画像631~633を示す。MIP画像611、612、613は、最大値投影法(maximum intensity projection)を用いて三次元分布画像610を、夫々、XY面、ZX面、YZ面に投影して得られる平面画像である。同様に、MIP画像621、622、623は、最大値投影法を用いて三次元分布画像620を、夫々、XY面、ZX面、YZ面に投影して得られる平面画像である。同様に、MIP画像631、632、633は、最大値投影法を用いて三次元分布画像630を、夫々、XY面、ZX面、YZ面に投影して得られる平面画像である。つまり例えば、三次元分布画像610における、X及びY軸方向の位置が共通の画素の群の内、最大の画素値を有する画素をXY面に投影する単位操作を、三次元分布画像610の全体に亘って行うことにより、投影面であるXY面にMIP画像611が形成される。MIP画像612及び613についても同様であり、また、MIP画像621~623及び631~633についても同様である。
 基準方法による分布画像では、画像空間ISの辺縁部においてオーバーフィッティングに起因した強いアーチファクトが出現しているのに対し(特に画像611を参照)、そのようなアーチファクトが第2改良方法及び第1改良方法では低減している(特に画像621及び631参照)。また、基準方法による分布画像では腫瘍組織の輪郭のぼけが顕著であるのに対し、そのようなぼけが第2改良方法及び第1改良方法では低減している。また、第2改良方法の分布画像で観測される背景のノイズ(例えば画像621中の領域NZ1及び画像622中のNZ2内のノイズ)が第1改良方法では低減しており、図25からは認識しづらいかもしれないが、第1改良方法の分布画像における各腫瘍組織(特にA431)のコントラストは第2改良方法のそれよりも高くなっていることが分かった。
 次に、各分布画像の画素の強度(画素値)が実際の組織に集積した放射能を反映したものであるかを評価するために、VOI(volume of interest)解析を行った。VOI解析において、図26に示す如く、三次元分布画像610に腫瘍組織A431、4T1及びC6が存在する領域610[A431]、610[4T1]及び610[C6]を設定し、且つ、三次元分布画像620に腫瘍組織A431、4T1及びC6が存在する領域620[A431]、620[4T1]及び620[C6]を設定し、且つ、三次元分布画像630に腫瘍組織A431、4T1及びC6が存在する領域630[A431]、630[4T1]及び630[C6]を設定する。腫瘍組織A431に対して設定される領域(即ち、領域610[A431]、620[A431]及び630[A431])の大きさは、画像610、620及び630間で同じである。腫瘍組織4T1及びC6に対して設定される領域についても同様である。
 図27は、撮像実験後における担癌マウスの各組織からの64Cuの放射能の計測結果を示している。図28において、バー651~653から成る棒グラフでは、バー651にて示される腫瘍組織A431からの計測放射能を1とする正規化を行った上で、腫瘍組織4T1、C6からの計測放射能を夫々バー652、653にて示す。
 図28のバー671~673が示す相対強度は、夫々、第1改良方法に基づく領域630[A431]、630[4T1]、630[C6]内の画素値の合計値を表す。但し、バー671~673が示す相対強度に関し、領域630[A431]に対応する相対強度を1.0とする正規化を行っている。
 図28のバー681~683が示す相対強度は、夫々、第2改良方法に基づく領域620[A431]、620[4T1]、620[C6]内の画素値の合計値を表す。但し、バー681~683が示す相対強度に関し、領域620[A431]に対応する相対強度を1.0とする正規化を行っている。
 図28のバー691~693が示す相対強度は、夫々、基準方法に基づく領域610[A431]、610[4T1]、610[C6]内の画素値の合計値を表す。但し、バー691~693が示す相対強度に関し、領域610[A431]に対応する相対強度を1.0とする正規化を行っている。
 図28のグラフにおいて、腫瘍組織A431の相対強度(1.0)から見た腫瘍組織4T1及びC6の相対強度が計測放射能についてのそれら(バー651の強度から見たバー652及び653の強度)に近いほど、画像の定量性が優れていると言える。
 バー691の強度から見たバー692及び693の強度との比較において、バー671の強度から見たバー672及び673の強度、並びに、バー681の強度から見たバー682及び683の強度は、共に、バー651の強度から見たバー652及び653の強度に近い。つまり、基準方法との比較において、第1又は第2改良方法では画像の定量性の改善が図られている。
 <<変形等>>
 本発明の実施形態は、特許請求の範囲に示された技術的思想の範囲内において、適宜、種々の変更が可能である。以上の実施形態は、あくまでも、本発明の実施形態の例であって、本発明ないし各構成要件の用語の意義は、以上の実施形態に記載されたものに制限されるものではない。上述の説明文中に示した具体的な数値は、単なる例示であって、当然の如く、それらを様々な数値に変更することができる。
 画像化装置1そのものである対象装置又は画像化装置1の一部(例えば主演算部30)である対象装置を、集積回路等のハードウェア、或いは、ハードウェアとソフトウェアの組み合わせによって構成することができる。対象装置にて実現される機能の全部又は一部である任意の特定の機能をプログラムとして記述して、該プログラムを対象装置に搭載可能なメモリ(フラッシュメモリ等)に保存しておいても良い。そして、該プログラムをプログラム実行装置(例えば、対象装置に搭載可能なマイクロコンピュータ)上で実行することによって、その特定の機能を実現するようにしてもよい。上記プログラムは任意の記録媒体に記憶及び固定されうる。上記プログラムを記憶及び固定する記録媒体は対象装置と異なる機器(サーバ機器等)に搭載又は接続されても良い。
 例えば、以下のように考えることができる。画像化装置1は、イベントを検出するイベント検出手段及び分布画像を生成する演算手段を備えている。イベント検出手段は主として例えばイベント測定データ取得部20により実現され、演算手段は主として例えば主演算部30により実現される(図3参照)。演算手段(30)は、イベントごとにコンプトン円錐面を設定する円錐面設定手段(S21;図22)、イベントごとにコンプトン円錐面領域を設定する円錐面領域設定手段(S22;図22)、イベントごとに且つ画素(対象画素)ごとに確率パラメータvijを設定する確率パラメータ設定手段(S23;図22)、及び、イベントごとに且つ画素(対象画素)ごとに検出感度パラメータsijを設定する検出感度パラメータ設定手段(S23;図22)を内包していると言える。
  1 ガンマ線源分布画像化装置
 10 コンプトンカメラ
 11 前段検出器
 12 後段検出器
 20 イベント測定データ取得部
 30 主演算部
200 撮像対象
201 ガンマ線源
1、Ri1 相互作用位置(コンプトン散乱位置)
2、Ri2 相互作用位置(全エネルギ吸収位置)
θC コンプトン散乱角
IS 画像空間

Claims (16)

  1.  後段検出器と、
     ガンマ線源と前記後段検出器との間に配置された前段検出器と、
     前記前段検出器でコンプトン散乱されたガンマ線が前記後段検出器にて光電吸収されるイベントを検出するイベント検出手段と、
     複数のイベントにおける各検出器とガンマ線との相互作用の測定データに基づき、前記ガンマ線源の空間分布を示す分布画像を、前記ガンマ線源を内包する画像空間内の画像として生成する演算手段と、を備えた画像化装置において、
     前記演算手段は、各イベントでコンプトン散乱したガンマ線が前記画像空間内から到来した確率を示す確率パラメータを各イベントの前記測定データに基づき前記分布画像内の複数の画素の夫々に対して個別に設定し、設定した前記確率パラメータを用いて前記分布画像を生成する
    ことを特徴とする画像化装置。
  2.  前記演算手段は、
     前記測定データに基づいて、前記イベントごとに、コンプトン散乱位置を頂点とし且つコンプトン散乱角を半頂角とし且つコンプトン散乱後のガンマ線の飛行経路の直線上に中心軸を持った円錐面を設定する円錐面設定手段と、
     前記イベントごとに、前記コンプトン散乱角の不確定性を示す散乱角分布関数に基づき前記円錐面に厚みを持たせることで円錐面領域を設定する円錐面領域設定手段と、
     前記イベントごとに且つ前記画素ごとに、前記コンプトン散乱位置を中心とし当該画素の位置及び前記コンプトン散乱位置間の距離を半径とする球面と前記円錐面領域と前記画像空間との交差状態に基づいて、前記確率パラメータを設定する確率パラメータ設定手段と、を有する
    ことを特徴とする請求項1に記載の画像化装置。
  3.  前記確率パラメータ設定手段は、前記イベントごとに且つ前記画素ごとに、前記球面と前記円錐面領域と前記画像空間との交差領域に配置された各画素を抽出して、抽出した各画素に対して前記散乱角分布関数の強度に応じた指標を導出し、抽出した各画素について導出された指標の総和に基づき前記確率パラメータを設定する
    ことを特徴とする請求項2に記載の画像化装置。
  4.  前記演算手段は、前記イベントごとに前記円錐面領域内の各画素に対して前記確率パラメータを設定する
    ことを特徴とする請求項2又は3に記載の画像化装置。
  5.  前記演算手段は、各イベントでのコンプトン散乱の検出感度を示す検出感度パラメータを各イベントの前記測定データに基づき前記複数の画素の夫々に対して個別に設定し、設定した前記検出感度パラメータ及び前記確率パラメータを用いて前記分布画像を生成する
    ことを特徴とする請求項1~4の何れかに記載の画像化装置。
  6.  前記演算手段は、前記イベントごとに且つ前記画素ごとに、コンプトン散乱位置と当該画素の位置とに応じて前記検出感度パラメータを設定する検出感度パラメータ設定手段を有する
    ことを特徴とする請求項5に記載の画像化装置。
  7.  前記検出感度パラメータ設定手段は、前記イベントごとに且つ前記画素ごとに、前記コンプトン散乱位置及び当該画素の位置に基づく前記前段検出器内におけるコンプトン散乱前のガンマ線の飛程と、前記コンプトン散乱位置及び当該画素の位置間の距離とに基づき、前記検出感度パラメータを設定する
    ことを特徴とする請求項6に記載の画像化装置。
  8.  前記演算手段は、各イベントでのコンプトン散乱の検出感度を示す検出感度パラメータを各イベントの前記測定データに基づき前記複数の画素の夫々に対して個別に設定する検出感度パラメータ設定手段を有して、前記検出感度パラメータ及び前記確率パラメータを用いて前記分布画像を生成し、
     前記検出感度パラメータ設定手段は、前記イベントごとに且つ前記画素ごとに、前記コンプトン散乱位置及び当該画素の位置に基づく前記前段検出器内におけるコンプトン散乱前のガンマ線の飛程と、前記コンプトン散乱位置及び当該画素の位置間の距離とに基づき、前記検出感度パラメータを設定し、
     前記演算手段は、
    Figure JPOXMLDOC01-appb-M000001
    に従う反復計算を介して、前記分布画像における各画素の画素値を導出し、
     λj (l)、λj (l+1)は、夫々、第l回目、第(l+1)回目の反復計算によって得られる、前記分布画像の第j番目の画素の画素値を表し、
     tij、vij、sijは、夫々、第iイベントの前記測定データに基づき、第iイベントにおける第j番目の画素に対して設定されるシステム応答関数、前記確率パラメータ、前記検出感度パラメータを表す
    ことを特徴とする請求項2~4の何れかに記載の画像化装置。
  9.  各イベントにおいて、前記測定データにより、前記前段検出器内におけるガンマ線のコンプトン散乱位置及びコンプトン散乱により当該ガンマ線が失ったエネルギ、並びに、前記後段検出器内におけるガンマ線の光電吸収位置及び光電吸収エネルギが示される
    ことを特徴とする請求項1~8の何れかに記載の画像化装置。
  10.  後段検出器と、
     ガンマ線源と前記後段検出器との間に配置された前段検出器と、
     前記前段検出器でコンプトン散乱されたガンマ線が前記後段検出器にて光電吸収されるイベントを検出するイベント検出手段と、
     複数のイベントにおける各検出器とガンマ線との相互作用の測定データに基づき、前記ガンマ線源の空間分布を示す分布画像を生成する演算手段と、を備えた画像化装置において、
     前記演算手段は、各イベントでのコンプトン散乱の検出感度を示す検出感度パラメータを各イベントの前記測定データに基づき前記分布画像内の複数の画素の夫々に対して個別に設定し、設定した前記検出感度パラメータを用いて前記分布画像を生成する
    ことを特徴とする画像化装置。
  11.  前記演算手段は、前記イベントごとに且つ前記画素ごとに、コンプトン散乱位置と当該画素の位置とに応じて前記検出感度パラメータを設定する検出感度パラメータ設定手段を有する
    ことを特徴とする請求項10に記載の画像化装置。
  12.  前記検出感度パラメータ設定手段は、前記イベントごとに且つ前記画素ごとに、前記コンプトン散乱位置及び当該画素の位置に基づく前記前段検出器内におけるコンプトン散乱前のガンマ線の飛程と、前記コンプトン散乱位置及び当該画素の位置間の距離とに基づき、前記検出感度パラメータを設定する
    ことを特徴とする請求項11に記載の画像化装置。
  13.  前記演算手段は、
     前記測定データに基づいて、前記イベントごとに、コンプトン散乱位置を頂点とし且つコンプトン散乱角を半頂角とし且つコンプトン散乱後のガンマ線の飛行経路の直線上に中心軸を持った円錐面を設定する円錐面設定手段と、
     前記イベントごとに、前記コンプトン散乱角の不確定性を示す散乱角分布関数に基づき前記円錐面に厚みを持たせることで円錐面領域を設定する円錐面領域設定手段と、を有し、前記イベントごとに前記円錐面領域内の各画素に対して前記検出感度パラメータを設定する
    ことを特徴とする請求項10~12の何れかに記載の画像化装置。
  14.  各イベントにおいて、前記測定データにより、前記前段検出器内におけるガンマ線のコンプトン散乱位置及びコンプトン散乱により当該ガンマ線が失ったエネルギ、並びに、前記後段検出器内におけるガンマ線の光電吸収位置及び光電吸収エネルギが示される
    ことを特徴とする請求項10~13の何れかに記載の画像化装置。
  15.  後段検出器及びガンマ線源と前記後段検出器との間に配置された前段検出器を備えたコンプトンカメラに利用される方法であって、前記前段検出器でコンプトン散乱されたガンマ線が前記後段検出器にて光電吸収される各イベントでの各検出器とガンマ線との相互作用の測定データに基づき、前記ガンマ線源の空間分布を示す分布画像を、前記ガンマ線源を内包する画像空間内の画像として生成する画像化方法において、
     各イベントでコンプトン散乱したガンマ線が前記画像空間内から到来した確率を示す確率パラメータを各イベントの前記測定データに基づき前記分布画像内の複数の画素の夫々に対して個別に設定し、設定した前記確率パラメータを用いて前記分布画像を生成する
    ことを特徴とする画像化方法。
  16.  後段検出器及びガンマ線源と前記後段検出器との間に配置された前段検出器を備えたコンプトンカメラに利用される方法であって、前記前段検出器でコンプトン散乱されたガンマ線が前記後段検出器にて光電吸収される各イベントでの各検出器とガンマ線との相互作用の測定データに基づき、前記ガンマ線源の空間分布を示す分布画像を生成する画像化方法において、
     各イベントでのコンプトン散乱の検出感度を示す検出感度パラメータを各イベントでの前記測定データに基づき前記分布画像内の複数の画素の夫々に対して個別に設定し、設定した前記検出感度パラメータを用いて前記分布画像を生成する
    ことを特徴とする画像化方法。
PCT/JP2015/074427 2014-09-05 2015-08-28 画像化装置及び方法 WO2016035706A1 (ja)

Priority Applications (3)

Application Number Priority Date Filing Date Title
JP2016546612A JP6607576B2 (ja) 2014-09-05 2015-08-28 画像化装置及び方法
EP15837605.3A EP3190431B1 (en) 2014-09-05 2015-08-28 Imaging device and method
US15/509,112 US9784857B2 (en) 2014-09-05 2015-08-28 Imaging device and method

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2014181498 2014-09-05
JP2014-181498 2014-09-05

Publications (1)

Publication Number Publication Date
WO2016035706A1 true WO2016035706A1 (ja) 2016-03-10

Family

ID=55439769

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2015/074427 WO2016035706A1 (ja) 2014-09-05 2015-08-28 画像化装置及び方法

Country Status (4)

Country Link
US (1) US9784857B2 (ja)
EP (1) EP3190431B1 (ja)
JP (1) JP6607576B2 (ja)
WO (1) WO2016035706A1 (ja)

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2018036234A (ja) * 2016-09-02 2018-03-08 国立大学法人群馬大学 情報処理装置、イメージング方法、及び、イメージングプログラム
JP2020016508A (ja) * 2018-07-24 2020-01-30 キヤノン株式会社 放射線検出器
JP2020529903A (ja) * 2017-06-15 2020-10-15 コンセッホ スペリオル デ インベスティガシオンス サイエンティフィカス(シーエスアイシー) 腫瘍学的診断およびリアルタイムガイド下生検に適したデュアル画像システム
CN112512427A (zh) * 2018-08-07 2021-03-16 美国西门子医疗系统股份有限公司 用于医学成像的自适应康普顿摄像机
JP7578450B2 (ja) 2020-10-07 2024-11-06 キヤノンメディカルシステムズ株式会社 核医学診断装置

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP6551003B2 (ja) * 2015-07-24 2019-07-31 三菱重工業株式会社 放射線測定装置及び放射線測定方法
CN111880211B (zh) * 2020-06-08 2022-11-15 中国科学院高能物理研究所 一种利用康普顿散射事例统计进行放射性核素识别的方法

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2010048699A (ja) * 2008-08-22 2010-03-04 Gunma Univ コンプトンカメラ

Family Cites Families (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5861627A (en) * 1997-06-24 1999-01-19 The University Of Utah Image reconstruction for compton camera including spherical harmonics
US20050139775A1 (en) * 2003-12-26 2005-06-30 Riken Gamma-ray detector and gamma-ray image pickup apparatus
WO2006058432A1 (en) * 2004-12-01 2006-06-08 Triumf, Operating As A Joint Venture By The Governors Of The University Of Alberta, The University Of British Columbia, Carleton System for selecting true coincidence events in positron emission tomography
JP4486623B2 (ja) 2006-08-11 2010-06-23 独立行政法人理化学研究所 コンプトン撮像カメラ
JPWO2008035708A1 (ja) * 2006-09-21 2010-01-28 株式会社日立メディコ 核医学診断装置
US8515011B2 (en) * 2009-06-02 2013-08-20 Mayo Foundation For Medical Education And Research System and method for dose verification radiotherapy
US10088580B2 (en) * 2012-05-31 2018-10-02 Minnesota Imaging And Engineering Llc Detector systems for radiation imaging

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2010048699A (ja) * 2008-08-22 2010-03-04 Gunma Univ コンプトンカメラ

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
TAKAHIRO IDA ET AL.: "Development of a list- mode iterative reconstruction for Ge semiconductor Compton camera with an accurate model of the uncertainty of the Compton scattering angle", JSAP AUTUMN MEETING KOEN YOKOSHU, vol. 75 th, no. 18a-B2-9, 1 September 2014 (2014-09-01), XP009500813 *
YUKIO OKABE ET AL.: "Study of Implementation of Maximum-likelihood Reconstruction for Compton Camera", IEICE TECHNICAL REPORT, vol. 109, no. 407, 21 January 2010 (2010-01-21), pages 47 - 51, XP055423561 *

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2018036234A (ja) * 2016-09-02 2018-03-08 国立大学法人群馬大学 情報処理装置、イメージング方法、及び、イメージングプログラム
JP2020529903A (ja) * 2017-06-15 2020-10-15 コンセッホ スペリオル デ インベスティガシオンス サイエンティフィカス(シーエスアイシー) 腫瘍学的診断およびリアルタイムガイド下生検に適したデュアル画像システム
JP7226827B2 (ja) 2017-06-15 2023-02-21 コンセッホ スペリオル デ インベスティガシオンス サイエンティフィカス(シーエスアイシー) 被験者の第1の画像および第2の画像を生成するシステム及びシステムの作動方法
JP2020016508A (ja) * 2018-07-24 2020-01-30 キヤノン株式会社 放射線検出器
JP7182930B2 (ja) 2018-07-24 2022-12-05 キヤノン株式会社 放射線検出器
CN112512427A (zh) * 2018-08-07 2021-03-16 美国西门子医疗系统股份有限公司 用于医学成像的自适应康普顿摄像机
CN112512427B (zh) * 2018-08-07 2024-05-24 美国西门子医疗系统股份有限公司 用于医学成像的自适应康普顿摄像机
JP7578450B2 (ja) 2020-10-07 2024-11-06 キヤノンメディカルシステムズ株式会社 核医学診断装置

Also Published As

Publication number Publication date
US20170261624A1 (en) 2017-09-14
EP3190431A4 (en) 2018-04-18
EP3190431A1 (en) 2017-07-12
JP6607576B2 (ja) 2019-11-20
JPWO2016035706A1 (ja) 2017-07-13
EP3190431B1 (en) 2019-03-20
US9784857B2 (en) 2017-10-10

Similar Documents

Publication Publication Date Title
JP6607576B2 (ja) 画像化装置及び方法
JP5582370B2 (ja) ガンマ線を利用する画像化装置、画像信号処理装置およびガンマ線測定データの画像処理方法
US7888651B2 (en) Method and system for using tissue-scattered coincidence photons for imaging
Tashima et al. 3D Compton image reconstruction method for whole gamma imaging
Etxebeste et al. CCMod: a GATE module for Compton camera imaging simulation
Pedemonte et al. A machine learning method for fast and accurate characterization of depth-of-interaction gamma cameras
Feng et al. Influence of Doppler broadening model accuracy in Compton camera list-mode MLEM reconstruction
JPWO2017209059A1 (ja) ガンマ線画像取得装置およびガンマ線画像取得方法
Ida et al. Accurate modeling of event-by-event backprojection for a germanium semiconductor Compton camera for system response evaluation in the LM-ML-EM image reconstruction method
Yao et al. Study of 3D fast Compton camera image reconstruction method by algebraic spatial sampling
Ding et al. Charged‐particle emission tomography
Wu et al. An accurate probabilistic model with detector resolution and Doppler broadening correction in list-mode MLEM reconstruction for Compton camera
Borja-Lloret et al. Influence of the background in Compton camera images for proton therapy treatment monitoring
Wang et al. Enhancing spatial resolution of 18F positron imaging with the Timepix detector by classification of primary fired pixels using support vector machine
Hirano et al. Monte Carlo simulation of scintillation photons for the design of a high-resolution SPECT detector dedicated to human brain
El Bitar et al. A detector response function design in pinhole SPECT including geometrical calibration
Lin et al. Recovering the triple coincidence of non-pure positron emitters in preclinical PET
Saaidi et al. Monte Carlo simulation of two siemens biograph PET/CT system using gate: image quality performance
Calderón et al. Design, development, and modeling of a compton camera tomographer based on room temperature solid statepixel detector
Giovagnoli Image reconstruction for three-gamma PET imaging
Liang et al. Spatial resolution recovery utilizing multi-ray tracing and graphic processing unit in PET image reconstruction
Chiang et al. Time of flight dual photon emission computed tomography
Feng Modeling and regularization in tomographic reconstruction for Compton camera imaging
Kagaya et al. Evaluating the capability of detecting recoil-electron tracks using an electron-tracking Compton camera with a silicon-on-insulator pixel sensor
Kalaitzidis From Monte Carlo PET Simulations to Reconstructed Images: Modelling and Optimisation for 68Ga Theragnostics

Legal Events

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

Ref document number: 15837605

Country of ref document: EP

Kind code of ref document: A1

ENP Entry into the national phase

Ref document number: 2016546612

Country of ref document: JP

Kind code of ref document: A

REEP Request for entry into the european phase

Ref document number: 2015837605

Country of ref document: EP

WWE Wipo information: entry into national phase

Ref document number: 2015837605

Country of ref document: EP

NENP Non-entry into the national phase

Ref country code: DE

WWE Wipo information: entry into national phase

Ref document number: 15509112

Country of ref document: US