WO2022158575A1 - 断層撮像のためのシステム、方法、プログラム、及びプログラムを記録した記録媒体 - Google Patents

断層撮像のためのシステム、方法、プログラム、及びプログラムを記録した記録媒体 Download PDF

Info

Publication number
WO2022158575A1
WO2022158575A1 PCT/JP2022/002247 JP2022002247W WO2022158575A1 WO 2022158575 A1 WO2022158575 A1 WO 2022158575A1 JP 2022002247 W JP2022002247 W JP 2022002247W WO 2022158575 A1 WO2022158575 A1 WO 2022158575A1
Authority
WO
WIPO (PCT)
Prior art keywords
matrix
system matrix
inverse
square
elements
Prior art date
Application number
PCT/JP2022/002247
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 JP2022576763A priority Critical patent/JPWO2022158575A1/ja
Priority to EP22742694.7A priority patent/EP4283343A1/en
Priority to US18/262,230 priority patent/US20240153162A1/en
Publication of WO2022158575A1 publication Critical patent/WO2022158575A1/ja

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/00—2D [Two Dimensional] image generation
    • G06T11/003—Reconstruction from projections, e.g. tomography
    • G06T11/006—Inverse problem, transformation from projection-space into object-space, e.g. transform methods, back-projection, algebraic methods
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/52—Devices using data or image processing specially adapted for radiation diagnosis
    • A61B6/5205—Devices using data or image processing specially adapted for radiation diagnosis involving processing of raw data to produce diagnostic data
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/02—Arrangements for diagnosis sequentially in different planes; Stereoscopic radiation diagnosis
    • A61B6/03—Computed tomography [CT]
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/58—Testing, adjusting or calibrating thereof
    • A61B6/582—Calibration
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/00—2D [Two Dimensional] image generation
    • G06T11/003—Reconstruction from projections, e.g. tomography
    • G06T11/005—Specific pre-processing for tomographic reconstruction, e.g. calibration, source positioning, rebinning, scatter correction, retrospective gating
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/58—Testing, adjusting or calibrating thereof
    • A61B6/582—Calibration
    • A61B6/585—Calibration of detector units
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01T—MEASUREMENT OF NUCLEAR OR X-RADIATION
    • G01T1/00—Measuring X-radiation, gamma radiation, corpuscular radiation, or cosmic radiation
    • G01T1/29—Measurement performed on radiation beams, e.g. position or section of the beam; Measurement of spatial distribution of radiation
    • G01T1/2914—Measurement of spatial distribution of radiation
    • G01T1/2985—In depth localisation, e.g. using positron emitters; Tomographic imaging (longitudinal and transverse section imaging; apparatus for radiation diagnosis sequentially in different planes, steroscopic radiation diagnosis)
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2211/00—Image generation
    • G06T2211/40—Computed tomography
    • G06T2211/424—Iterative
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2211/00—Image generation
    • G06T2211/40—Computed tomography
    • G06T2211/436—Limited angle

Definitions

  • the present invention relates to a system for tomographic imaging, a program, and a recording medium on which the program is recorded.
  • Tomography which utilizes penetrating electromagnetic waves and particle beams such as X-rays, gamma rays, and light waves, or seismic waves and any other penetrating waves or particles, has been put to practical use.
  • Typical examples thereof are X-ray CT scanners, SPECT (Single Photon Emission Computed Tomography), Optical Tomography, and Electron Tomography (TEMT) used in medical and industrial applications.
  • Radon transform is the mathematical principle by which tomographic images of objects can be obtained by these imaging techniques.
  • the process of extracting the attenuation rate of each part inside the object as imaging information to the outside is described by the Radon transform
  • the process of reconstructing the tomographic information of the object from the imaging information extracted to the outside is described by the Radon inverse transform.
  • the Radon transform/inverse Radon transform in continuous coordinates and continuous orientations is performed in sampled coordinates and orientations to deal with a finite number of pixels in real images.
  • the Radon inversion (discrete Radon inversion) for image reconstruction under sampled coordinates and orientations is reduced to solving a multi-dimensional system of linear equations, ie matrix operations.
  • ART Algebraic Reconstruction Technique
  • ML-EM Maximum Likelihood - Expectation Maximization
  • ASIRT Additive - Simultaneous Iterative Reconstruction Techniques
  • MSIRT Multiplicative - Simultaneous Iterative Reconstruction Techniques
  • ART is superior in terms of the quality of reconstructed images, but even if the conventional ART is adopted in a situation where the number of imaging pixels is increased, the number of elements (variables to be determined) of the simultaneous linear equations to be solved is limited. As the number increases, enormous computational resources are required, and the processing takes a long time.
  • Non-Patent Document 5 has a scanning angle range of only 100°
  • the TEMT disclosed in Non-Patent Document 6 has a scanning angle range of only 160°.
  • the present invention provides a system, method, program, and recording medium storing the program that can generate a reconstructed image supported by an algebraically exact solution of the discrete Radon inverse transform with smaller computational resources.
  • One of the purposes is to
  • the present invention can generate reconstructed images backed by an algebraically exact solution of the discrete Radon inverse transform with less computational resources in a scan angle range of less than (180°-180°/number of projections). It is an object to provide a system, method, program, and recording medium recording the program.
  • One aspect of the present invention is a computer-implemented method for determining operating parameters of a tomographic imaging system that reconstructs tomographic images by an inverse discrete Radon transform in a scan angle range of less than (180°-180°/number of projections). based on the values of specified one or two operating parameters of the number of detector elements, the number of projections, and a scan angle range of less than (180°-180°/number of projections), the system matrix generating a square system matrix from the determined system matrix; determining whether or not an inverse matrix exists for the square system matrix; and repeatedly determining values of operating parameters corresponding to operating parameters other than the specified one or two for which there is an inverse of the square system matrix.
  • the step of generating the square system matrix can be a step of first removing duplicate row vector or column vector elements from the determined system matrix to generate a square system matrix.
  • the step of generating the square system matrix first removes duplicate row vector or column vector elements from the determined system matrix, and then randomly selects row vector or column vector elements for the remaining row vector or column vector elements. Removing the vector may be the step of generating a square system matrix.
  • One aspect of the present invention provides a program for causing a computer to execute the method.
  • One aspect of the present invention provides a computer-readable recording medium on which the program is recorded.
  • One aspect of the present invention is an apparatus for determining operating parameters of a tomographic imaging system that reconstructs tomographic images by inverse discrete Radon transform in a scan angle range of less than (180°-180°/number of projections), , the number of detector elements, the number of projections, and a scan angle range of less than (180°-180°/number of projections).
  • a square system matrix generation unit that generates a square system matrix from the system matrix; and an inverse matrix determination unit that determines whether or not the square system matrix has an inverse matrix. Generating a system matrix and determining the presence or absence of the inverse matrix are repeated as necessary for operating parameters other than the specified one or two for which an inverse matrix of the square system matrix exists.
  • An apparatus is provided for determining the value of an operating parameter that
  • the square system matrix generation unit may first remove duplicate row vector or column vector elements from the determined system matrix to generate a square system matrix.
  • the square system matrix generation unit first removes duplicate row vector or column vector elements from the determined system matrix, and then randomly generates a row vector or column vector for the remaining row vector or column vector elements. can be removed to produce a square system matrix.
  • One aspect of the invention is a computer-implemented method for determining operating parameters of a tomographic imaging system that reconstructs a tomographic image by an inverse discrete Radon transform over a scan angle range of predetermined values less than 180°. determining a system matrix based on the values of specified one or two operating parameters of the number of detector elements, the number of projections, and a scan angular range of predetermined values less than 180°; generating a square system matrix from a matrix; determining whether an inverse matrix exists for the square system matrix; and repeating each of the above steps as necessary to obtain an inverse matrix of the square system matrix determining a value of an operating parameter corresponding to operating parameters other than the specified one or two for which there is a .
  • One aspect of the present invention is an apparatus for determining operating parameters of a tomographic imaging system for reconstructing a tomographic image by means of an inverse discrete Radon transform in a scan angle range of predetermined values of less than 180°, wherein the number of detector elements is , the number of projections, and a scan angular range of predetermined values less than 180°; a square system matrix generation unit that generates a system matrix; and an inverse matrix determination unit that determines whether or not the square system matrix has an inverse matrix, determines the system matrix, generates the square system matrix, and determines the inverse matrix.
  • Waves include any waves that can be considered to be attenuated by propagation, such as electromagnetic waves and sound waves, and include light waves (infrared, visible, ultraviolet) and vibrations (earthquakes, etc.).
  • Particles include any particle that can be considered to be attenuated by propagation in a stream of particle beams, including neutrons and other radiation.
  • One or both of these waves and particles also includes those that have both wave and particle properties in terms of quantum mechanics, such as electrons and photons.
  • quantum mechanics such as electrons and photons.
  • waves or particles may be generated by man-made sources of their emission (irradiators), or may be naturally occurring, e.g. Even those that occur as properties of the object such as radiation emitted from each part of the object (referred to in this application as "waves or particles generated without irradiation equipment") can be used to implement the present disclosure .
  • the captured image at a scan angle of 0° and the captured image at a scan angle of 180° are the same. It is used to mean the scan angle range from the angle 0° to the scan end angle (180°-180°/number of projections). Therefore, in this specification and claims, a scan angle range smaller than the "scan angle range of 180°" in common expression is defined as "a scan angle of less than (180°-180°/number of projections) It is expressed as "range”.
  • a system, a receiver, a transmitter, a method, a program capable of generating a reconstructed image backed by an algebraically rigorous solution of the inverse discrete Radon transform with smaller computational resources.
  • a recording medium recording the program can be provided.
  • FIG. 1 is a perspective view showing the relationship between an object from which an image is acquired and a detection device in a tomographic imaging system as an example of an embodiment of the present invention
  • FIG. FIG. 2 is a schematic diagram for explaining a schematic configuration including a planar arrangement of a cut plane of an object from which an image is acquired in a tomographic imaging system as an example of an embodiment of the present invention, and a tomographic image according to the embodiment of the present invention
  • 1 is a diagram showing the overall configuration of an imaging system
  • FIG. It is an explanatory view explaining the principle of the algebraic technique of the conventional discrete Radon inverse transform.
  • FIG. 10 is a diagram showing an example of the relationship between pixels of a tomographic image, wave or particle beam irradiation, and the positions of the detection elements when the number of detection elements is greater than the resolution.
  • FIG. 10 is a diagram showing the relationship between pixels of an example tomographic image, wave or particle beam irradiation, and detection element positions when the number of projections is greater than the resolution.
  • FIG. 10 is a diagram showing the relationship between pixels of an example tomographic image, wave or particle beam irradiation, and detection element positions when the number of detection elements and the number of projections are larger than the resolution.
  • FIG. 10 illustrates the resolution and scan angle range for which exact solutions exist for additional numbers of detector elements and projections for resolutions up to 16; FIG.
  • FIG. 10 illustrates the resolution and scan angle range for which exact solutions exist for the additional number of detector elements and projections for resolutions up to 64;
  • FIG. 4 is a diagram showing an example of an efficient row vector selection method in matrix decimation processing according to the embodiment of the present invention;
  • 1 is a block diagram showing the functional configuration of an operating parameter calculation device according to an embodiment of the present invention;
  • FIG. 1 is a block diagram showing functional configurations of a discrete Radon inversion matrix generation device and a reconstructed image generation device according to an embodiment of the present invention;
  • FIG. 1 is a diagram showing a hardware configuration of an operating parameter calculation device 30 according to an embodiment of the present invention;
  • 4 is a flowchart of an example of operation parameter calculation processing according to the embodiment of the present invention
  • 4 is a flowchart of an example of operation parameter calculation processing according to an embodiment of the present invention
  • 4 is a flowchart of an example of a discrete Radon inversion matrix generation process according to an embodiment of the present invention
  • 4 is a flowchart of an example of a discrete Radon inversion matrix generation process according to an embodiment of the present invention
  • 6 is a flowchart of an example of reconstructed image generation processing according to the embodiment of the present invention
  • 6 is a flowchart of an example of reconstructed image generation processing according to the embodiment of the present invention
  • It is a figure which shows the verification data (8-bit numerical value phantom) used for verification.
  • FIG. 10 is a diagram showing an example of a sinogram when the number of detection elements is made larger than the resolution; It is a reconstructed image obtained in an example in which the number of detection elements is made larger than the resolution.
  • 15B is a reconstructed image obtained using a Shepp-Logan filter in the conventional FBP method based on the sinogram of FIG. 15A;
  • FIG. 15A is a reconstructed image obtained by performing 10 iterative approximation calculations (ML-EM method) as a conventional algebraic reconstruction method based on the sinogram of FIG. 15A.
  • FIG. 15B is a reconstructed image obtained without a filter in the FBP method as a reference based on the sinogram of FIG. 15A;
  • FIG. 10 is a diagram showing an example of a sinogram when the number of projections is made larger than the resolution; It is a reconstructed image obtained in an embodiment in which the number of projections is made larger than the resolution.
  • 16B is a reconstructed image obtained using a Shepp-Logan filter in the conventional FBP method based on the sinogram of FIG. 16A;
  • FIG. 16A is a reconstructed image obtained by performing 10 iterative approximation calculations (ML-EM method) as a conventional algebraic reconstruction method based on the sinogram of FIG. 16A.
  • FIG. 16B is a reconstructed image obtained without a filter in the FBP method as a reference, based on the sinogram of FIG. 16A;
  • FIG. 10 is a diagram showing an example of a sinogram when the number of detection elements and the number of projections are made larger than the resolution; It is a reconstructed image obtained in an example in which the number of detection elements and the number of projections are made larger than the resolution.
  • 17B is a reconstructed image obtained using a Shepp-Logan filter in the conventional FBP method based on the sinogram of FIG. 17A;
  • FIG. 17A is a reconstructed image obtained by performing 10 iterative approximation calculations (ML-EM method) as a conventional algebraic reconstruction method based on the sinogram of FIG. 17A.
  • FIG. 17B is a reconstructed image obtained without a filter in the FBP method as a reference, based on the sinogram of FIG. 17A;
  • FIG. 1A is a perspective view showing the relationship between an object 2 from which an image is acquired and a detection device 20 in a tomographic imaging system 100 as an example of this embodiment.
  • FIG. 1B is a schematic diagram for explaining a schematic configuration including a planar arrangement on a cross section of the object 2 from which an image is acquired in the tomographic imaging system 100 as an example of the present embodiment.
  • a group of pixels having N ⁇ N pixels forming a tomographic image is defined on one plane of this space, and this integer N is the resolution.
  • this integer N is matched to the number of detector elements 22 in detector 20 .
  • a beam 4 of waves or particles is contained in that plane and is emitted from an illumination device 10 .
  • the detection device 20 and the illumination device 10 can be rotated with respect to the object 2 while their relative positions are fixed.
  • the tomographic imaging system 100 includes a control unit 501, which will be described later, for the operation of the apparatus including such control. This rotation is relative to the group of pixels or object 2 . For this reason, it is not necessarily the case that only the detection device 20 and the irradiation device 10 are rotated. It may be rotated.
  • the operation of irradiating the wave or particle beam 4 with the irradiation device 10 and detecting the beam 4 on the opposite side of the object 2 with the detection device 20 is performed by dividing the scanning range in the detection direction ⁇ by the number of pixels N in the conventional method. It is sampled for the detection direction in increments. The sampling number for this detection direction is hereinafter also referred to as the "projection number". In the case where the beam 4 is parallel as shown in FIG. The detection signal is sampled at (180/N)° steps.
  • a typical scanning operation in the detection direction ⁇ of the tomographic imaging system 100, including the present embodiment, is a constant speed rotation that increases or decreases at a constant speed, and the timing of acquiring the value of the detection device 20 is electrically controlled. Sampling is performed by controlling However, this embodiment is not limited to such a thing.
  • FIG. 2 is an explanatory diagram for explaining the principle of a conventional algebraic technique for inverse discrete Radon transform.
  • This figure is an elementary one showing the essence of the process of irradiating waves or particles and reconstructing a tomographic image using intensity data after transmission.
  • the numerical value indicated by the arrow indicates the value (for example, attenuation) obtained when the detector element of the detection device is positioned at the tip of the arrow. It is indicated numerically. That is, the arrow represents a beam of waves or particles, and the pixels through which the arrow passes and the direction of the arrow indicate irradiation and detection (hereinafter simply referred to as "irradiation").
  • the total value at the tip of the arrow corresponds to the attenuation component detected by each detection element when the detection element of the detection device is at the tip of the arrow. Determining the four values a 1 to a 4 only by the values indicated by the arrows, ie, the values obtained from the detection device, is performed by executing the inverse discrete Radon transform based on the values indicated by the detection elements of the detection device. corresponds to the process of reconstructing the
  • Equation (1) remains unchanged even if new a' 1 to a' 4 shifted by an arbitrary number A are substituted.
  • the equation (1) has four elements (variables to be determined), while the actual number of simultaneous equations that can be derived by the detection indicated by the arrow is three.
  • the four values a1 to a4 are all indeterminate and cannot be specified.
  • w nm are the elements of matrix W (referred to as the “system matrix”), which is a square matrix with N 2 rows and N 2 columns.
  • the system matrix is a square matrix with N 2 rows and N 2 columns.
  • the values of a 1 to a 4 can be obtained from the values of p 1 to p 4 . That is, it is possible to uniquely determine a 1 to a 4 only by an algebraically rigorous method of applying this inverse matrix to the measured values p 1 to p 4 .
  • This difference suggests that, if appropriate measures are taken at the irradiation stage, there is still room for algebraic methods to solve even multidimensional simultaneous linear equations of a more general form for reconstruction processing of real tomographic images. That is, there is room for obtaining non-indefinite exact solutions by algebraic methods.
  • the inventors have found that even if the scan angle range is less than (180°-180°/number of projections), it may be possible to obtain an exact solution.
  • the three operation parameters ie, the number of detection elements, the number of projections, which is the number of samples in the detection direction, and the scanning angle range.
  • the number of samples such as the number of detection elements and the number of projections
  • the scanning angle range should be adjusted.
  • the probability of existence of the exact solution increases. It was found that as the probability increases, the range in which the exact solution exists increases, and the range in which the exact solution exists expands in the direction where the scan angle range is small.
  • the beam 4 has an extent that covers the entire pixel range at the above coordinates and has a certain width ⁇ in the portion directed to each detector element 22 .
  • the contribution that the pixel has to the beam 4 on the detector element 22 in question is typically the width ⁇ of the beam 4 towards one of the detector elements 22 and the overlap of each pixel (vertical and horizontal ⁇ ). . In FIG. 1B, this overlap is indicated by hatching for one pixel.
  • Contributions of all pixels of N ⁇ N pixels are determined for each detection direction ⁇ and for each detection element 22 .
  • the system matrix W holds this contribution as a numerical value.
  • the system matrix W the array of sampling points for the space corresponding to the total number of pixels is arranged in the row direction, and the array of sampling points for irradiation and detection corresponding to the step of the detection direction and the detection element is arranged in the column direction. It is normal to take When the conventional method is applied to the setting of the number of pixels in FIG. 1B, the system matrix W becomes a square matrix having N ⁇ N elements in both row and column directions, and the total number of elements is N 4 . Become.
  • the number of pixels is N pixels ⁇ N pixels
  • the scanning range in the detection direction ⁇ is N samples
  • the number of detection elements 22 is N elements.
  • the value of N has increased with the times in order to improve image quality. is not limited to a specific value.
  • N is about 10 3
  • the system matrix has a size of about 10 6 in both row and column directions, and the number of elements is about 10 12 .
  • X in equation (4) represents an attenuation image, which is an image representing attenuation (including absorption and scattering) at each pixel position, as a column vector.
  • the attenuation image in the actual space is the two-dimensional image of N pixels ⁇ N pixels shown in FIG. 1B. is X.
  • the column vector X' on the left side obtained by applying the system matrix W to the attenuation image X in Equation (4) corresponds to an image called a sinogram.
  • the column vector X' also indicates the detection intensity obtained for each of the sampling points and the detection elements in the detection direction.
  • the sinogram is obtained by rearranging the elements of the column vector X' two-dimensionally with the axis of the detection element and the axis of the increments in the detection direction.
  • Equation (4) expresses illumination and detection.
  • the column vector X corresponds to each value of a 1 to a 4
  • the column vector X' corresponds to the value attached to each arrow
  • the system matrix W corresponds to each arrow and each It can be said that it reflects the correspondence between pixels.
  • Reconstruction of a tomographic image is a process of calculating a column vector X from the system matrix W and the sinogram column vector X', and obtaining an attenuation image from the column vector X, on the premise that equation (4) is realized.
  • X W -1X ' (5)
  • W ⁇ 1 is the inverse matrix of W
  • Equation (5) is obtained by applying W ⁇ 1 to both sides of Equation (4) and transposing the left and right sides.
  • the square matrix W does not always have the inverse matrix W ⁇ 1 . The fact that the values of a 1 to a 4 are not fixed in FIG.
  • FIG. 2 corresponds to the fact that W -1 does not exist in the way of irradiation and detection.
  • this embodiment can be used when the irradiation device 10 and the detection device 20 that utilize non-parallel beams such as fan beams and cone beams are used. can also be applied with slight modifications obvious to those skilled in the art. Also, the method of the present embodiment can be applied to settings other than the settings such as the number of pixels in FIG. 1B.
  • equations (4) and (5) can also be expressed in other mathematically equivalent forms, e.g., for the system matrix W, space
  • the array of sampling points for is taken in the column direction
  • the array of sampling points for irradiation and detection corresponding to the increments in the detection direction and the detection elements is taken in the row direction.
  • X in equation (4) can also be expressed as a row vector according to the position of each pixel.
  • the equations corresponding to equations (4) and (5) are obtained by transposing each side. This transpose operation has the expressive effect of swapping the rows and columns of a matrix or vector, changing the order in which the matrix operated on the vector operates from the right rather than from the left.
  • the inventors have found that an exact solution can sometimes be obtained even if the scan angle range is less than (180°-180°/number of projections).
  • the inventors found that it is sufficient to adjust the three operating parameters: the number of detection elements, the number of projections, which is the number of samples in the detection direction, and the scanning angle range.
  • the number of samples such as the number of detection elements and the number of projections
  • the scanning angle range is made larger than the resolution of the tomographic image.
  • the probability of existence of the exact solution increases. It was found that as the probability increases, the range in which the exact solution exists increases, and the range in which the exact solution exists expands in the direction where the scan angle range is small.
  • the number of detection elements is larger than the resolution
  • the number of projections is larger than the resolution
  • a case where the number of detection elements and the number of projections are larger than the resolution will be described below.
  • the number of detection elements is larger than the resolution
  • the resolution is 2, that is, the tomographic image is 2 pixels ⁇ 2 pixels
  • the number of detection elements is 4
  • the number of projections is 2
  • the scanning angle range is 22°. A case will be described.
  • Equation (6) the geometry will be as shown in FIG. 3, so that the values p 1 -p 8 detected by each detector element are the attenuation of the beam by the object 2 at each pixel location.
  • equation (8) the system matrix W is given by equation (8).
  • the system matrix W is not a square matrix, a 1 to a 4 cannot be obtained by multiplying the inverse matrix from the left in Equation (7). Therefore, the system matrix W is squared by thinning out four row vector elements among the eight row vector elements forming the system matrix W.
  • FIG. This squared matrix is called the square system matrix W SQ .
  • the square system matrix W SQ has nine inverse combinations. One of them is the case of choosing the dashed row vector, in which case equation (9) holds for the resulting square system matrix W SQ . Then, by multiplying the inverse matrix W SQ -1 from the left, a 1 to a 4 can be obtained as in equation (10), and a tomographic image can be obtained.
  • the resolution is 2, that is, the tomographic image is 2 pixels x 2 pixels, the number of detection elements is 2, the number of projections is 6, and the scanning angle range is 65.2°. A case will be described.
  • Equation (11) the geometry is as shown in FIG. 4, so that the values p 1 -p 8 detected by each detector element are the attenuation of the beam by object 2 at each pixel location.
  • equation (12) the system matrix W is given by equation (13).
  • the square system matrix W SQ has 160 inverse combinations. One of them is the case of choosing the dashed row vector, in which case equation (14) holds for the resulting square system matrix W SQ .
  • equation (14) holds for the resulting square system matrix W SQ .
  • the resolution is 2, that is, the tomographic image is 2 pixels ⁇ 2 pixels, the number of detection elements is 4, the number of projections is 4, and the scanning angle range is A case of 54° will be described.
  • Equation (16) the geometry will be as shown in FIG. 5, so the values p 1 -p 16 detected by each detector element will be the attenuation of the beam by the object 2 at each pixel location.
  • equation (17) the system matrix W is given by equation (18).
  • W SQ the system matrix W SQ has an inverse matrix.
  • equation (19) holds for the resulting square system matrix W SQ .
  • 6A is a diagram showing the resolution and scan angle range for which there is an exact solution for the number of detection elements and the number of projections added to the resolution for resolutions up to 16, and FIG. 6B for detection for resolutions up to 64.
  • FIG. 10 illustrates the resolution and scan angle range for which exact solutions exist for the number of elements and number of projections added to the resolution;
  • FIG. 6B the illustration when both the number of detection elements and the number of projections are added is omitted.
  • 6A and 6B are obtained by changing the number of additional detection elements, the number of additional projections, the scanning angle range, and the resolution under the same computational experiment conditions as in (Computer Simulation Verification (Example)) described later. . As can be seen from FIGS.
  • the scan angle range is less than (180°-180°/number of projections), and an exact solution can be obtained by appropriately selecting the number of detection elements, the number of projections, and the scan angle range. can be done. Also, if the number of samples, such as the number of detection elements and the number of projections, is larger than the resolution of the tomographic image, the probability of existence of the exact solution increases. Along with this, the range in which the exact solution exists increases, and the range in which the exact solution exists expands in the direction in which the scan angle range is small. Thus, for example, an exact solution can be obtained in a scan angle range of less than the desired (180°-180°/number of projections) by appropriately choosing the number of detector elements and the number of projections.
  • the scan angle range must be 160° or less.
  • the reconstructed image obtained by this embodiment is a high-quality image based on an exact solution, it is useful as learning data for a machine learning system that predicts a reconstructed image from an input tomographic image. be.
  • the above is generalized and explained as follows.
  • N is a positive integer
  • the number of detecting elements 22 of the detecting device 20 is N, but in this embodiment is (N+j) (j is an integer equal to or greater than 0).
  • the detection direction ⁇ is sampled in steps of (180/N)° so that the N directions do not overlap in the range of 0 to 180°, that is, the scan angle range of 180°. Sampling is done.
  • sampling of detection directions in this embodiment is performed in non-overlapping N+k directions (k is an integer equal to or greater than 0) in a scan angle range ⁇ of less than (180° ⁇ 180°/number of projections), for example, ( ⁇ /( N+k-1)) in degrees.
  • Sampling intervals increments in the detection direction
  • the detection directions are sampled in N+k directions instead of N directions.
  • a system matrix is determined based on the number of detection elements (N+j), the number of projections (N+k), and the scan angle range ⁇ .
  • the system matrix determined here has (N+k) ⁇ (N+j) elements in the column direction corresponding to the detection direction and the detection element, and N ⁇ N elements in the row direction corresponding to the number of pixels of the tomographic image. .
  • the elements of the (j+k) ⁇ N+j ⁇ k row vector with N ⁇ N elements at the appropriate positions of the matrix are extracted from the system matrix W to produce a square system matrix W SQ .
  • the square system matrix W SQ has an inverse if the elements of the row vector to be removed are appropriate.
  • a discrete Radon inverse transform matrix which is the inverse matrix W SQ -1 , is calculated from the square system matrix W SQ .
  • any appropriate known method such as the cofactor method or the sweep method can be used.
  • ⁇ Reconstructed image generation processing> On the other hand, to obtain a tomographic image from the object 2, the object 2 is arranged between the irradiation device 10 and the detection device 20 having each detection element 22. FIG. While irradiating the beam 4 from the irradiation device 10, detection is performed in non-overlapping N+k directions.
  • this beam 4 continues to irradiate continuously, and while the directions (detection directions) of the object 2, irradiation device 10, and detection device 20 move smoothly, the signal from the detection element 22 of the detection device 20 is electrically sampled.
  • a typical example of illumination and detection in non-overlapping N+k directions is sampling at ( ⁇ /(N+k ⁇ 1))° steps.
  • the obtained signals can be organized into data that can generate a sinogram, if processing such as correction and scale conversion for the physical characteristics of the equipment and rearrangement of the arrangement order are performed as necessary.
  • the process of obtaining this sinogram is the irradiation and detection process, and is the discrete Radon transform.
  • the first vector X' has (N+k).times.(N+j) elements since it has sinogram values.
  • Equation ( 4 ) a second vector X SQ ', which is a column vector with N ⁇ N elements, can be obtained. Since the second vector X SQ ' satisfies the relationship of Equation (4) with the square system matrix W SQ , the equation ( The relationship 5) is established, and the third vector X SQ that gives a reconstructed image is obtained by discrete Radon inversion processing. A reconstructed image can be visualized by a devectorization process of returning the arrangement of the third vectors X SQ to an image of N pixels ⁇ N pixels by the inverse operation of vectorization. A tomographic image for one slice is thus reconstructed.
  • the processing of this embodiment is divided into two systems, the above-described discrete Radon inversion matrix generation processing and reconstructed image generation processing.
  • a user who normally takes reconstructed images performs the reconstructed image generation process described above, from illumination and detection to reconstructed image generation.
  • the user does not necessarily have to execute the discrete Radon inverse transform matrix generation process from the system matrix generation corresponding to the above-described operating parameters to the discrete Radon inverse transform matrix generation.
  • Discrete Radon inverse transform matrix generation processing is required in situations such as manufacturing, installation, and maintenance of a tomographic imaging system, which is the kind of processing handled by device manufacturers.
  • the reconstructed image generation processing is the permutation of the order of matrices and vectors and the calculation of matrix multiplication, and the processing such as iterative approximation, for which the amount of calculation processing becomes a problem, is not included in the range.
  • the reconstructed image generation processing used in the user's normal imaging is the elements of the row vector to be thinned out and the discrete Radon inversion matrix W SQ -1 prepared by the discrete Radon inversion matrix generation processing. It can be executed with a sufficiently light calculation load just by using it, and it can be said that this point brings high practicality to this embodiment. This is because the computational burden of the most frequent processing of acquiring a tomographic image is extremely small.
  • the method of the present embodiment is extremely advantageous for applications in which imaging of many slices is repeated to obtain a three-dimensional volume image.
  • the reconstruction image generation process is suitable for computational processing on many-core type processors such as GPGPU (General-purpose computing on graphics processing unit). This is because an improvement in the processing speed of GPGPU (General-purpose computing on graphics processing unit).
  • system matrix W the square system matrix W SQ , and the inverse matrix W SQ -1 of the system matrix W SQ , which reflect the geometric aspect, are included only in the discrete Radon inverse transform matrix generation process, and are used for measurement.
  • the separation of the reconstruction image generation process from the reconstruction image generation process essentially separates the problem of convergence when using iteration in algebraic methods and the noise aspect of the measurement. For this reason, in this embodiment, there is also the advantage that a phenomenon caused by an approximate solution that is unrelated to noise is less likely to occur, and that various noise reduction processes that are generally employed in image processing are likely to exhibit their intended effects.
  • the total number of rows in system matrix W is (N+k) ⁇ (N+j). Of these, the number of rows left to make the matrix a square matrix is N ⁇ N, and the total number of row vectors to be thinned out is (j+k) ⁇ N+j ⁇ k. Determine either the sequence of N ⁇ N numbers to be left or (j+k) ⁇ N+j ⁇ k numbers to be thinned out from the integer range [1, (N+k) ⁇ (N+j)] of indices (row numbers) that identify rows do. A randomized approach is used for this determination. That is, either N ⁇ N pieces to be left or (j+k) ⁇ N+j ⁇ k pieces to be thinned out is randomly determined by numerically generated random numbers.
  • a square system matrix W SQ can be obtained from the system matrix W by corresponding decimation processing.
  • the rank of the obtained square system matrix W SQ is calculated.
  • the calculated rank is the same value as N ⁇ N when the square system matrix W SQ has an inverse matrix, and a value smaller than N ⁇ N when it does not have an inverse matrix.
  • the rank value can determine whether the square system matrix W SQ can have an inverse. If the rank does not match N ⁇ N, another random number is used to select the decimated number sequence ⁇ m ⁇ by a random algorithm, calculate the rank, and determine whether or not it has an inverse matrix.
  • the matrix decimation process can be performed at a practical speed.
  • the row vector corresponding to p 1 is (1010), while the row vectors corresponding to p 2 , p 5 and p 9 are the same row vector (1010). Therefore, the row vectors corresponding to p 2 , p 5 and p 9 are thinned out.
  • the row vector corresponding to p 3 is (0101), but the row vector corresponding to p 8 and p 12 is the same row vector (0101). Therefore, row vectors corresponding to p 8 and p 12 are thinned out.
  • the row vector corresponding to p 6 is (1110), but the row vector corresponding to p 10 and p 14 is the same row vector (1110). Therefore, row vectors corresponding to p 10 and p 14 are thinned out.
  • row vector corresponding to p 7 is (0111), but the row vector corresponding to p 11 and p 15 is the same row vector (0111). Therefore, row vectors corresponding to p 11 and p 15 are thinned out.
  • the square system matrix W SQ composed of the elements of the row vectors corresponding to p 2 , p 3 , p 10 , p 11 as described above, has an inverse matrix, but as described above, corresponding to p 1 Since the elements of the row vector are elements of the same row vector as the element of the row vector corresponding to p 2 , the square system matrix W SQ composed of p 1 , p 3 , p 10 and p 11 is p 2 , It is identical to the square system matrix W SQ composed of the elements of the row vectors corresponding to p 3 , p 10 and p 11 and has an inverse.
  • the reconstructed image obtained from p 2 , p 3 , p 10 and p 11 and p 1 , p 3 , p 10 and p 11 By superimposing the reconstructed image obtained from , the SN ratio of the reconstructed image is improved and the image quality is improved. This is because it is generally said that the greater the number of measurements, the higher the measurement accuracy, and the SN ratio becomes n 1/2 times when measurements are performed n times.
  • the row vector elements corresponding to p 5 and p 9 are the same row vector elements as the row vector element corresponding to p 2 .
  • FIG. 1B is a diagram showing the overall configuration of the tomographic imaging system according to the embodiment of the present invention.
  • the tomography system 1 includes an irradiation device 10 , a detection device 20 , an operating parameter calculation device 30 , a discrete Radon inversion matrix generation device 40 and a reconstructed image generation device 50 .
  • the irradiation device 10, the detection device 20, the operating parameter calculation device 30, the discrete Radon inversion matrix generation device 40, and the reconstructed image generation device 50 may not be configured as one physical device, and the irradiation device 10 and A part or all of the detection device 20, the motion parameter calculation device 30, the discrete Radon inversion matrix generation device 40, and the reconstructed image generation device 50 may be configured as separate physical devices.
  • the devices configured as separate physical devices may be connected to each other via a wired/wireless computer network.
  • Each of the irradiation device 10, the detection device 20, the operating parameter calculation device 30, the discrete Radon inversion matrix generation device 40, and the reconstructed image generation device 50 does not need to be configured as one physical device, and can be configured as a plurality of physical devices. It may be configured from a typical device.
  • FIG. 8 is a block diagram showing the functional configuration of the operating parameter calculation device according to the embodiment of the present invention.
  • the operating parameter calculation device 30 includes an operating parameter interpreting unit 301 , an operating parameter receiving unit 303 , a system matrix determining unit 305 , a square system matrix generating unit 307 , an inverse matrix determining unit 309 , and an operating parameter calculation result output unit 311 .
  • the operation parameter interpreter 301 specifies the value, and if the input data is a set of values and/or a range of values of the operation parameter, , to specify one of the values.
  • the operational parameter reception unit 303 receives the values of one or two designated operational parameters out of the number of detection elements, the number of projections, and the scan angle range less than (180°-180°/number of projections).
  • the system matrix determination unit 305 sets the specified one or two operating parameter values and the operating parameters other than the operating parameters corresponding to the specified one or two operating parameter values to predetermined values. Based on this, determine the system matrix.
  • the inverse matrix determination unit 309 determines whether or not an inverse matrix exists for the determined square system matrix.
  • the operating parameter calculation result output unit 311 outputs predetermined values of operating parameters other than the operating parameters corresponding to the specified one or two operating parameter values in the form of display, voice output, or the like.
  • FIG. 9 is a block diagram showing the functional configuration of the discrete Radon inversion matrix generation device and reconstructed image generation device according to the embodiment of the present invention.
  • the discrete Radon inversion matrix generation device 40 includes an operating parameter reception unit 403 , a system matrix determination unit 405 , a square system matrix generation unit 407 , an inverse matrix determination unit 409 , a discrete Radon inversion matrix generation unit 411 and a storage unit 413 .
  • the operating parameter accepting unit 403 accepts the values of the three designated operating parameters.
  • the system matrix determination unit 405 determines the system matrix based on the values of the three designated operation parameters.
  • the inverse matrix determination unit 409 determines whether or not an inverse matrix exists for the determined square system matrix.
  • the discrete Radon inverse transform matrix generator 411 generates a discrete Radon inverse transform matrix W SQ ⁇ 1 that is an inverse matrix of the square system matrix W SQ .
  • the storage unit 413 determines that the elements of the row vectors that make up the square system matrix W SQ are the same as the elements of the row vector that were removed in step S207, the storage unit 413 makes up the square system matrix W SQ .
  • the reconstructed image generation device 50 includes a control unit 501, a sinogram generation unit 503, a first vector generation unit 505, a vector decimation unit 507, a third vector generation unit 509, a devectorization unit 511, a different measurement value second vector generation unit 513, A composite reconstructed image generation unit 515 and a storage unit 517 are provided.
  • the control unit 501 controls the irradiation device 10 and the detection device 20 to perform detection by the detection device 20 in non-overlapping N+j directions while irradiating the object 2 with the beam 4 from the irradiation device 10 .
  • the sinogram generation unit 503 performs processing such as correction of the physical characteristics of the device, conversion of scale, and rearrangement of the arrangement order, etc., as necessary, for the signals obtained from the detection elements 22 of the detection device 20, and generates a sinogram. do.
  • the first vector generation unit 505 generates a first vector X', which is a column vector obtained by rearranging the elements from the generated sinogram in an order matching the column direction of the system matrix W.
  • a vector decimation unit 507 thins out the elements of the first vector X' corresponding to the elements of the row vector removed when generating the square system matrix W SQ to generate a second vector X SQ '.
  • the third vector generation unit 509 uses the discrete Radon inversion matrix W SQ -1 generated by the discrete Radon inversion matrix generation unit 411 and the second vector X SQ ′ to generate the discrete Radon A third vector X SQ that provides a reconstructed image is generated by inverse transform processing.
  • the devectorization unit 511 generates a reconfigured image by devectorization processing of returning the arrangement of the third vectors X SQ to an image of N pixels ⁇ N pixels by the inverse operation of vectorization.
  • Different measurement value second vector generation section 513 configures square system matrix W SQ by square system matrix generation section 407 based on the index (row number) stored in storage section 413 by square system matrix generation section 407.
  • the index (row number) of the row vector constituting the square system matrix W SQ when it is determined that the row vector elements that are the same as the row vector elements removed by the square system matrix generation unit 407 exist among the row vector elements.
  • a composite reconstructed image generation unit 515 generates a composite reconstructed image with improved image quality by superimposing a plurality of generated reconstructed images.
  • the storage unit 517 stores the reconstructed image generated in step S319.
  • FIG. 10 is a diagram showing the hardware configuration of the operating parameter calculation device 30 according to the embodiment of the present invention.
  • the operating parameter calculation device 30 includes a CPU 30a, a RAM 30b, a ROM 30c, an external memory 30d, an input section 30e, an output section 30f, and a communication section 30g.
  • the RAM 30b, ROM 30c, external memory 30d, input section 30e, output section 30f, and communication section 30g are connected to the CPU 30a via a system bus 30h.
  • the CPU 30a comprehensively controls each device connected to the system bus 30h.
  • the ROM 30c and the external memory 30d store the BIOS and OS, which are control programs for the CPU 30a, and various programs and data necessary for realizing the functions executed by the computer.
  • the RAM 30b functions as the main memory and work area of the CPU.
  • the CPU 30a loads necessary programs and the like from the ROM 30c and the external memory 30d to the RAM 30b and executes the loaded programs to implement various operations.
  • the external memory 30d is composed of, for example, flash memory, hard disk, DVD-RAM, USB memory, and the like.
  • the input unit 30e accepts operation instructions and the like from users and the like.
  • the input unit 30e includes input devices such as input buttons, a keyboard, a pointing device, a wireless remote controller, and a microphone.
  • the output unit 30f outputs data processed by the CPU 30a and data stored in the RAM 30b, ROM 30c, and external memory 30d.
  • the output unit 30f is composed of an output device such as a CRT display, LCD, organic EL panel, printer, speaker, or the like.
  • the communication unit 30g is an interface for connecting and communicating with external devices via a network or directly.
  • the communication unit 30g is composed of interfaces such as a serial interface and a LAN interface, for example.
  • the hardware configurations of the discrete Radon inversion matrix generation device 40 and the reconstructed image generation device 50 are the same.
  • 11A and 11B are flow charts of an example of operation parameter calculation processing according to an embodiment of the present invention.
  • the user sets a value, a set of values, or a range of values for one or two of the number of detector elements, the number of projections, and the scan angle range less than (180°-180°/number of projections) as the operating parameters.
  • Input through the input unit 30e of the calculation device 30 (S101). If the input data is an operation parameter value, the operation parameter interpretation unit 301 of the operation parameter calculation device 30 designates the value, and the input data is a set of values and/or a range of values of the operation parameter. , one of the values is specified (S103). If the input data is a set of values and/or a range of values for an operating parameter, the following processing is repeated until the set of values and/or a range of values for the input operating parameter are exhausted. Therefore, the case where the input data is the value of the operation parameter will be described.
  • the operation parameter value reception unit 303 of the operation parameter calculation device 30 receives the values of the designated one or two operation parameters (S105).
  • the system matrix determination unit 305 determines predetermined values (initial values ), the system matrix W is determined (S107).
  • the system matrix W corresponds to each of the detection directions of the (N+k) directions, where the number of detection elements and the number of projections are (N+j) and (N+k) (where j and k are integers equal to or greater than 0).
  • (N+k) having a weight value indicating the contribution of each pixel of the two-dimensional tomographic image of N pixels ⁇ N pixels to the path of the part of the wave or particle heading to each of the (N+j) detecting elements 22 (N+k) x(N+j) columns or row vectors having N x N elements are arranged in the row or column direction so as to correspond to the pixels of the two-dimensional tomographic image.
  • duplicate row vector elements are removed from the row vector elements that make up the determined system matrix W (S109).
  • a square system matrix W SQ is generated by removing elements of the column vector (S111). Specifically, for example, for the index (row number) that identifies the row of the elements of the remaining row vector, either the index of the element of the row vector to be left or the index of the element of the row vector to be thinned out is numerically generated. Random numbers are used to determine the decimation number sequence ⁇ m ⁇ . Then, the square system matrix W SQ is generated by removing the elements of the row vector corresponding to the determined thinning sequence ⁇ m ⁇ .
  • the inverse matrix determination unit 309 determines whether or not an inverse matrix exists for the generated square system matrix W SQ (S113). Specifically, for example, the rank of the square system matrix W SQ is calculated, and whether or not the calculated rank is N ⁇ N determines whether or not the inverse matrix exists.
  • the operation parameter calculation result output unit 311 displays predetermined values of operation parameters other than the operation parameters corresponding to the specified one or two operation parameter values, outputs audio, etc. (S115).
  • step S107 If it is determined that the inverse matrix does not exist, the process returns to step S107, and for different predetermined values of the operation parameters other than the operation parameter corresponding to the specified one or two operation parameter values, do the steps.
  • the inverse matrix determination unit 309 returns to step S107 as necessary, and determines different predetermined values for the operation parameters other than the operation parameter corresponding to the specified one or two operation parameter values after step S107. Steps may be repeated. For example, even if it is determined that the inverse matrix exists in step S113, by repeating, the operating parameters other than the operating parameter corresponding to the specified one or two operating parameter values that have the inverse matrix range can be obtained.
  • ⁇ Generation Processing of Inverse Discrete Radon Transform Matrix> 12A and 12B are flowcharts of an example of a process for generating an inverse discrete Radon transform matrix according to an embodiment of the present invention.
  • the user can determine the number of detection elements N+j, the number of projections N+k, (180°-180 °/number of projections) are specified via the input section of the discrete Radon inversion matrix generator 40 (S201).
  • the operation parameter value reception unit 403 of the discrete Radon inversion matrix generation device 40 receives the values of the three designated operation parameters (S203).
  • the system matrix determination unit 405 determines the system matrix W based on the values of the three designated operation parameters (S205).
  • the square system matrix W SQ is generated by removing the elements of the column vector (S209). Then, if it is determined that among the elements of the row vectors that make up the square system matrix W SQ , there are elements that are the same as the elements of the row vector that were removed in step S207, then the elements of the row vectors that make up the square system matrix W SQ are An index (row number) specifying a row corresponding to the element and the removed element of the same row vector is stored in the storage unit 413 (S211).
  • the inverse matrix determination unit 409 determines whether or not an inverse matrix exists for the generated square system matrix W SQ (S213). If it is determined that an inverse matrix exists (S213-Yes), the discrete Radon inverse transform matrix generator 411 generates a discrete Radon inverse transform matrix W SQ -1 that is a square system matrix W SQ inverse matrix (S215). If it is determined that the inverse matrix does not exist (S213-No), return to step S209 and generate a square system matrix W SQ based on another random number (S217).
  • ⁇ Reconstructed image generation processing> 13A and 13B are flowcharts of an example of reconstructed image generation processing according to an embodiment of the present invention.
  • the object 2 is arranged between the irradiation device 10 and the detection device 20 (S301).
  • the control unit 501 of the reconstructed image generation device 50 controls the irradiation device 10 and the detection device 20 to irradiate the target object 2 with the beam 4 from the irradiation device 10, and detect by the detection device 20 in non-overlapping N+j directions. (S303).
  • the sinogram generation unit 503 of the reconstructed image generation device 50 processes the signals obtained from the detection elements 22 of the detection device 20, such as correcting the physical characteristics of the device and converting the scale, and organizes the arrangement order as necessary. etc. to generate a sinogram (S305).
  • the process of generating this sinogram is the irradiation and detection process, and is the discrete Radon transform.
  • the first vector generation unit 505 of the reconstructed image generation device 50 generates a first vector X′, which is a column vector obtained by rearranging the elements from the generated sinogram in an order matching the column direction of the system matrix W. (S307).
  • the first vector X' has (N+k).times.(N+j) elements since it has sinogram values. From irradiation and detection S301 to vectorization S307, quantization by analog-to-digital conversion of signals, processing of digital signals, and processing by temporarily storing data as necessary can be performed.
  • the sinogram it is sufficient for the sinogram to be obtained as long as the data that can generate it is obtained, and it is not necessary to display it, and it is not necessary to go through explicit data that can be said to be a sinogram. If the first vector corresponding to the vectorized sinogram is obtained, arbitrary processing (for example, processing by simulation) that performs equivalent arbitrary processing from irradiation and detection S301 to vectorization S307 is also included in the present embodiment. ing.
  • the vector decimation unit 507 of the reconstructed image generation device 50 thins out (vector decimation process) to generate the second vector X SQ ' (S309).
  • vector decimation processing it is possible to obtain the second vector X SQ ', which is a column vector having N ⁇ N elements that can satisfy the relationship of equation (4) with the square system matrix W SQ . .
  • the third vector generation unit 509 of the reconstructed image generation device 50 generates a discrete A third vector X SQ that provides a reconstructed image is generated by Radon inversion processing (S311).
  • the devectorization unit 511 of the reconstructed image generation device 50 generates a reconstructed image by devectorization processing of returning the arrangement of the third vectors X SQ to an image of N pixels ⁇ N pixels by the inverse operation of vectorization (S313). A tomographic image for one slice is thus reconstructed.
  • step S211 the separate measurement value second vector generation unit 513 of the reconstructed image generation device 50 generates row vectors forming the square system matrix W SQ based on the index (row number) stored in the storage unit 411 in step S211.
  • elements of the second vector X SQ ' corresponding to the elements of the row vectors forming the square system matrix W SQ when it is determined that the elements of the row vector removed in step S207 exist. are replaced with the elements of the second vector X SQ ' corresponding to the removed elements of the same row vector to generate another measured value second vector (S315).
  • the third vector generation unit 509 generates the third vector X SQ using the different measurement value second vector (S317). Then, the devectorization unit 511 performs devectorization processing on the generated third vector X SQ to generate a reconstructed image and stores it in the storage unit 517 (S319).
  • steps S315 to S319 are performed when it is determined in step S211 that among the elements of the row vectors forming the square system matrix W SQ , there are elements that are the same as the elements of the row vector removed in step S207. Repeat a desired number of times for different combinations of the row vectors that make up SQ and elements of the same removed row vectors (S321).
  • the synthesized reconstructed image generating unit 515 of the reconstructed image generating device 50 superimposes the generated reconstructed images to generate a synthesized reconstructed image with improved image quality (S323).
  • FIG. 14 is a diagram (referred to as an original image) showing verification data (8-bit numerical phantom) used for verification.
  • 15A, 16A, and 17A show the cases where the number of detection elements is made larger than the resolution, the number of projections is made larger than the resolution, and the number of detection elements and the number of projections are made larger than the resolution, respectively.
  • FIG. 10 is a diagram showing an example of a sinogram corresponding to . In each figure, numerical values indicating the beginning and end of pixel numbers representing pixel coordinates are shown on the upper and left sides.
  • areas with strong absorption are displayed brightly, and areas with weak absorption are displayed darkly.
  • the operation of irradiation and detection by the tomographic imaging apparatus 10 was also simulated, and numerical data corresponding to the measured values were obtained by calculation.
  • the sinogram representing this weak values of detection intensity are displayed brightly, and strong values are displayed darkly.
  • the sinograms of Figures 15A, 16A, and 17A were obtained based on the discrete Radon transforms for each measurement operation. Specifically, the following three patterns were confirmed.
  • the sinogram of FIG. 16A is displayed as an image with a size of 64 ⁇ 65 pixels.
  • the sinogram of FIG. 17A is displayed as an image with a size of 128 ⁇ 67 pixels.
  • addition of detection elements is performed by reducing the size of the detection elements so that the area of the detection element group does not change.
  • 15A, 16A, and 17A are for the operation of this embodiment, but are also used for conventional FBP and ML-EM processing.
  • the conventional method it is not particularly necessary to add detection elements or increase the number of projections, and a reconstructed image can be obtained by adopting a sinogram (not shown) acquired with a size of N ⁇ N pixels.
  • the sinogram with the size of N ⁇ N pixels has less information than the sinogram obtained by the present embodiment.
  • a sinogram with a size of (N+j) ⁇ (N+k) pixels according to the present embodiment is used for comparison between methods using images reconstructed from sinograms having the same amount of information.
  • the number of pixels in the detection element direction is the resolution.
  • a sinogram was used in which the size of the image was reduced by averaging the pixel values to be identical.
  • step S301 to S305 are simulated by computer simulation as described above
  • the process employed to determine the decimation sequence ⁇ m ⁇ for the matrix decimation process is the randomized algorithm described above.
  • the rank did not match (reach) N ⁇ N (64 ⁇ 64), and the process returned to step S209 and retried.
  • the ratio of retry processing actually required is (1) about 1 to 2 times on average for 10 algorithm operations when the number of detection elements is greater than the resolution, and (2) the number of projections is greater than the resolution. (3) When the number of detection elements and the number of projections are larger than the resolution, the average is about 3287 times for 10 algorithm operations. It takes about 10 to 20 minutes when the number of detection elements is larger than the resolution, and (2) the number of projections is larger than the resolution. 3) It took about 8 to 24 hours when the number of detection elements and projections was larger than the resolution. If the addition of the number of projections is as small as about 1 and the measurement angle range is small, the total number of rows does not increase significantly.
  • the rank condition is no longer satisfied, so the number of retries to select rows that satisfy the rank condition may increase compared to a matrix with a large number of additional projections and additional detection elements.
  • the calculation processing time is not determined only by the number of rows, but also changes depending on the dependency of each row. Especially when the scan angle range is narrow, it is considered that the independence of each row is reduced, so such an increase or decrease in processing time is a theoretically assumed behavior.
  • FIGS. 15B, 16B, and 17B are reconstructed images obtained in respective examples of patterns (1), (2), and (3) above.
  • 15C, 16C, and 17C are reconstructed images obtained using the Shepp-Logan filter in the conventional FBP method based on the sinograms of FIGS. 15A, 16A, and 17A.
  • 15D, 16D, and 17D are based on the sinograms of FIGS. 15A, 16A, and 17A, and iterative approximation calculation (ML-EM method) is performed 10 times as a conventional algebraic reconstruction method, and the local solution is 0. It is a reconstructed image obtained up to the number of times that does not occur.
  • 15E, 16E, and 17E are reconstructed images obtained without using a filter in the FBP method as a reference, based on the sinograms of FIGS. 15A, 16A, and 17A.
  • FIGS. 15E, 16E and 17E show the ISNR values of the reconstructed images in FIGS. 15B to 15D, FIGS. 16B to 16D, and FIGS. , Tables 4 and 5.
  • a scan angle range smaller than the “180° scan angle range” in the conventional expression is defined as “(180°-180°/number of projections) less than the scan angle range” in this specification and claims. It is expressed as "angle range”.
  • the specified scan angle range was described as being less than (180°-180°/number of projections), but the specified scan angle range is any suitable predetermined value less than 180° (for example, 170° or less, 165° or less, etc.).
  • the scan angle range is less than (180°-180°/number of projections), but if the number of additional detection elements j is 1 or more and/or the number of additional projections k is 1 or more, Even with a scan angle range of 180°, it is possible to generate a reconstructed image backed by an algebraically exact solution of the discrete Radon inverse transform with less computational resources.
  • the method of determining the contribution of the pixel to the beam 4 to the detection element 22 described with reference to FIG. 1B can also be employed in this embodiment.
  • the operation of the tomographic imaging apparatus 100 has been described as an example, but the present embodiment is not necessarily limited to employing the irradiation apparatus 10 .
  • the present embodiment can be similarly applied to other techniques for performing tomographic imaging of an object without using an irradiation device, such as SPECT (single photon emission tomography).
  • SPECT single photon emission tomography
  • a computer-implemented method for determining operating parameters of a tomographic imaging system that reconstructs an N pixel by N pixel tomographic image by an inverse discrete Radon transform over a scan angle range of less than (180°-180°/number of projections) and
  • a target object arranged in a detection range of a detection device having a plurality of detection elements arranged in at least one line is directed to the detection device by an irradiation device, and waves to be detected by the detection device are detected.
  • the waves or particles after passing through each part of the object are detected by each of the detection elements.
  • the step of generating the square system matrix first removes duplicate row or column vector elements from the determined system matrix, and selects a predetermined total of (j+k) ⁇ N+j ⁇ k row or column vectors The method of paragraph (1), wherein the step of generating a square system matrix by removing elements.
  • the step of generating the square system matrix first removes duplicate row vector or column vector elements from the determined system matrix, and then randomly selects row vector or column vector elements for the remaining row vector or column vector elements.
  • Method according to clause (2) the step of removing vectors to produce a square system matrix by removing elements of a predetermined total of (j+k) ⁇ N+j ⁇ k row or column vectors.
  • a computer-readable recording medium recording the program according to item (4).
  • a wave to be detected by the detection device is directed toward the detection device by the irradiation device with respect to an object arranged in a detection range of the detection device having a plurality of detection elements arranged in at least one row.
  • the waves or particles after passing through each part of the object are detected by each of the detection elements.
  • the square system matrix generator first removes duplicate row vector or column vector elements from the determined system matrix, and generates a predetermined total of (j+k) ⁇ N+j ⁇ k row vector or column vector elements. Apparatus according to clause (6) for generating a square system matrix by elimination.
  • the square system matrix generation unit first removes duplicate row vector or column vector elements from the determined system matrix, and then randomly generates a row vector or column vector for the remaining row vector or column vector elements.
  • the present disclosure is applicable for any tomographic imaging technique in which waves or particles are transmitted through the object itself, e.g. , SPECT (Single Photon Emission Computed Tomography), Optical Tomography, Optical CT apparatus, and Electron Beam Tomography (TEMT).
  • SPECT Single Photon Emission Computed Tomography
  • Optical Tomography Optical CT apparatus
  • ETT Electron Beam Tomography

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Health & Medical Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Medical Informatics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Optics & Photonics (AREA)
  • Surgery (AREA)
  • High Energy & Nuclear Physics (AREA)
  • Veterinary Medicine (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Public Health (AREA)
  • Pathology (AREA)
  • Radiology & Medical Imaging (AREA)
  • Biomedical Technology (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Molecular Biology (AREA)
  • Biophysics (AREA)
  • Animal Behavior & Ethology (AREA)
  • General Health & Medical Sciences (AREA)
  • Algebra (AREA)
  • Computational Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Computing Systems (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Apparatus For Radiation Diagnosis (AREA)
  • Image Processing (AREA)

Abstract

より小さな計算資源で、離散ラドン逆変換の代数的厳密解に裏付けられた再構成画像を生成することができるシステム、方法、プログラム、及びプログラムを記録した記録媒体を提供すること。　（１８０°－１８０°／投影数）未満のスキャン角度範囲で、離散逆ラドン変換により断層画像を再構成する断層撮像システムの動作パラメータを求めるためのコンピュータにより実行される方法であって、検出素子数、投影数、及び（１８０°－１８０°／投影数）未満のスキャン角度範囲のうちの指定された１つ又は２つの動作パラメータの値に基づいて、システム行列を決定するステップと、決定された前記システム行列から正方システム行列を生成するステップと、前記正方システム行列について、逆行列が存在するか否かを判定するステップと、前記各ステップを、必要に応じて繰り返して、前記正方システム行列の逆行列が存在する、前記指定された１つ又は２つ以外の動作パラメータに対応する動作パラメータの値を求めるステップとを含む方法。

Description

断層撮像のためのシステム、方法、プログラム、及びプログラムを記録した記録媒体
　本発明は、断層撮像のためのシステム、プログラム、及びプログラムを記録した記録媒体に関する。
　Ｘ線、ガンマ線、光波などの透過性の電磁波や粒子線、または地震波やその他の任意の透過性を示す波又は粒子を利用する断層撮像法（トモグラフィー）が実用化されている。その代表例が、医療、工業用途で用いられるＸ線ＣＴスキャナー、ＳＰＥＣＴ（Single Photon Emission Computed Tomography）、光トモグラフィー(Optical Tomography)、電子線トモグラフィー（ＴＥＭＴ）である。これらの撮像手法で対象物の断層画像が得られる数学的原理はラドン変換である。この場合、対象物内部の各部の減衰率を撮像情報として外部に取出す処理がラドン変換によって記述され、対象物の断層情報を外界に取り出した撮像情報から再構成する処理がラドン逆変換によって記述される。すなわち、対象物で減衰した透過後の波又は粒子の流れを照射方向又は検出方向を変えて各位置で検出する操作がラドン変換に対応する。逆にその検出した方向別位置別の強度情報から対象物内部の各部の減衰率を推定して画像を再構築する操作がラドン逆変換に対応する。連続座標や連続方位でのラドン変換／ラドン逆変換は、現実の画像では有限画素数を扱うためサンプリングされた座標や方位で実行される。サンプリングした座標や方位の下での画像再構成のためのラドン逆変換（離散ラドン逆変換）は多元の連立一次方程式を解くことつまり行列の演算に還元される。これらの数学的側面については標準的教科書に解説されている（例えば非特許文献１、２）。
　反復処理を行う逐次近似法（iterative reconstruction, ＩＲ）によってよりノイズやアーティファクトが軽減される手法も注目されている。その反復は、例えば１００回又はそれ以上繰返すことにより解を得ている。このように代数的手法に頼って多元の連立一次方程式の近似解を得ることによって再構成像を得る手法はＡＲＴ（Algebraic Reconstruction Technique）と呼ばれている。ＡＲＴは、例えば、Maximum Likelihood - Expectation Maximization (ML-EM）法、Additive - Simultaneous Iterative Reconstruction Techniques (ASIRT) 法、Multiplicative - Simultaneous Iterative Reconstruction Techniques (MSIRT) 法、と呼ばれる手法を含んでいる。従来のＡＲＴでは、計算量が膨大となり多量の計算資源を必要とする。また、イタレーションを重ねても、真の画像に対応する解に近づかず、局所解に捕獲されてしまうことも多い。他方、代数的手法にこだわらずに計算機での近似計算の効率を高め、現実的な計算資源で対処する手法も模索されてきた。この一例がＦＢＰ（Filtered Back Projection）である（例えば特許文献１）。ＦＢＰは代数的手法に比べ計算量が少ないことから現在Ｘ線ＣＴスキャナーにおいて最も普及している。
　断層撮像法の発展経過では、初期から今日まで撮像画素数増大の要請が継続しており、ＦＢＰの可能性が探求されてきた。しかしＦＢＰは、ＩＲに比して種々のアーティファクトが顕われやすく再構成像の品質では一般に不利である。計算能力の発展が続く今日では、再構成像の品質で本来優位なＡＲＴも再び注目されている（例えば特許文献２）。
国際公開第２０１４／０２１３４９号 国際公開第２０１３／１０５５８３号
Aviash C. Kak and Molcolm Slaney, "Principles of Computerized Tomographic Imaging," IEEE Press (1988) Stanley R. Deans, "The Radon Transform and Some of Its Applications," Krieger Publications (1983) Robert L. Siddon, "Fast calculation of the exact radiological　path for a three-dimensional CT array," Medical Physics, vol.12, no.2, pages 252-255 (1985) Sunnegardh, J. and Danielsson, P.-E. "A New Anti-Aliased Projection Operator for Iterative CT Reconstruction," 9 th International Meeting of Fully Three-Dimensional Image Reconstruction in Radiology and Nuclear Medicine, pages 125-127 (2009) Serhan O. Isikman, Waheb Bishara, Sam Mavandadi, Frank W. Yu, Steve Feng, Randy Lau, and Aydogan Ozcan"Lens-free optical tomographic microscope with a large imaging volume on a chip"PNAS May 3, 2011 108 (18) 7296-7301; https://doi.org/10.1073/pnas.1015638108 N. Kawase, M. Kato, H. Nishioka, H. Jinnai, Ultramicroscopy, 107, 8-15 (2007)"Transmission electron microtomography without the "missing wedge" for quantitative structural analysis"
　上述のように、再構成像の品質に関してはＡＲＴが優れているが、撮像画素数が増大した状況で従来のＡＲＴを採用しても、解くべき連立一次方程式の元（決定すべき変数）の数が増大してしまい膨大な計算資源が必要となり、その処理に長い時間を要する。
　また、断層撮像法においては、一般的に、スキャン角度範囲として１８０°必要であるが、電子線トモグラフィー（ＴＥＭＴ）や光学ＣＴ顕微鏡においては、スキャン角度範囲を１８０°取ることができず、例えば、非特許文献５に開示される光学ＣＴ顕微鏡では、スキャン角度範囲が１００°しか取れず、また、非特許文献６に開示されるＴＥＭＴでは、スキャン角度範囲が１６０°しか取れない。
　そこで、本発明は、より小さな計算資源で、離散ラドン逆変換の代数的厳密解に裏付けられた再構成画像を生成することができるシステム、方法、プログラム、及びプログラムを記録した記録媒体を提供することを目的の１つとする。
　また、本発明は、（１８０°－１８０°／投影数）未満のスキャン角度範囲で、より小さな計算資源で、離散ラドン逆変換の代数的厳密解に裏付けられた再構成画像を生成することができるシステム、方法、プログラム、及びプログラムを記録した記録媒体を提供することを目的の１つとする。
　本発明の１つの態様は、（１８０°－１８０°／投影数）未満のスキャン角度範囲で、離散逆ラドン変換により断層画像を再構成する断層撮像システムの動作パラメータを求めるためのコンピュータにより実行される方法であって、検出素子数、投影数、及び（１８０°－１８０°／投影数）未満のスキャン角度範囲のうちの指定された１つ又は２つの動作パラメータの値に基づいて、システム行列を決定するステップと、決定された前記システム行列から正方システム行列を生成するステップと、前記正方システム行列について、逆行列が存在するか否かを判定するステップと、前記各ステップを、必要に応じて繰り返して、前記正方システム行列の逆行列が存在する、前記指定された１つ又は２つ以外の動作パラメータに対応する動作パラメータの値を求めるステップとを含む方法を提供するものである。
　前記正方システム行列を生成するステップは、決定された前記システム行列から、まず重複する行ベクトル又は列ベクトルの要素を除去し、正方システム行列を生成するステップであるものとすることができる。
　前記正方システム行列を生成するステップは、決定された前記システム行列から、まず重複する行ベクトル又は列ベクトルの要素を除去し、次いで、残った行ベクトル又は列ベクトルの要素についてランダムに行ベクトル又は列ベクトルを除去して、正方システム行列を生成するステップであるものとすることができる。
　本発明の１つの態様は、前記方法をコンピュータに実行させるためのプログラムを提供するものである。
　本発明の１つの態様は、前記プログラムを記録したコンピュータ読み取り可能な記録媒体を提供するものである。
　本発明の１つの態様は、（１８０°－１８０°／投影数）未満のスキャン角度範囲で、離散逆ラドン変換により断層画像を再構成する断層撮像システムの動作パラメータを求めるための装置であって、検出素子数、投影数、及び（１８０°－１８０°／投影数）未満のスキャン角度範囲のうちの指定された１つ又は２つの動作パラメータの値に基づいてシステム行列を決定するシステム行列決定部と、前記システム行列から正方システム行列を生成する正方システム行列生成部と、前記正方システム行列について、逆行列の有無を判定する逆行列判定部とを備え、前記システム行列を決定し、前記正方システム行列を生成し、前記逆行列の有無を判定することを、必要に応じて繰り返して、前記正方システム行列の逆行列が存在する、前記指定された１つ又は２つ以外の動作パラメータに対応する動作パラメータの値を求める装置を提供するものである。
　前記正方システム行列生成部は、決定された前記システム行列から、まず重複する行ベクトル又は列ベクトルの要素を除去し、正方システム行列を生成するものとすることができる。
　前記正方システム行列生成部は、決定された前記システム行列から、まず重複する行ベクトル又は列ベクトルの要素を除去し、次いで、残った行ベクトル又は列ベクトルの要素についてランダムに行ベクトル又は列ベクトルを除去して、正方システム行列を生成するものとすることができる。
　本発明の１つの態様は、１８０°未満の所定の値のスキャン角度範囲で、離散逆ラドン変換により断層画像を再構成する断層撮像システムの動作パラメータを求めるためのコンピュータにより実行される方法であって、検出素子数、投影数、及び１８０°未満の所定の値のスキャン角度範囲のうちの指定された１つ又は２つの動作パラメータの値に基づいて、システム行列を決定するステップと、前記システム行列から正方システム行列を生成するステップと、前記正方システム行列について、逆行列が存在するか否かを判定するステップと、前記各ステップを、必要に応じて繰り返して、前記正方システム行列の逆行列が存在する、前記指定された１つ又は２つ以外の動作パラメータに対応する動作パラメータの値を求めるステップとを含む方法を提供するものである。
　本発明の１つの態様は、１８０°未満の所定の値のスキャン角度範囲で、離散逆ラドン変換により断層画像を再構成する断層撮像システムの動作パラメータを求めるための装置であって、検出素子数、投影数、及び１８０°未満の所定の値のスキャン角度範囲のうちの指定された１つ又は２つの動作パラメータの値に基づいてシステム行列を決定するシステム行列決定部と、前記システム行列から正方システム行列を生成する正方システム行列生成部と、前記正方システム行列について、逆行列の有無を判定する逆行列判定部とを備え、前記システム行列を決定し、前記正方システム行列を生成し、前記逆行列の有無を判定することを、必要に応じて繰り返して、前記正方システム行列の逆行列が存在する、前記指定された１つ又は２つ以外の動作パラメータに対応する動作パラメータの値を求める装置を提供するものである。
　本明細書及び特許請求の範囲においては、不明瞭にならない限り本開示の属する技術分野の慣用に従う用語法を採用することがある。波は電磁波、音波といった伝播による減衰を観念できる任意の波動を含み、光波（赤外、可視、紫外）、振動（地震など）を含む。粒子は、流れとなった粒子線において伝播による減衰を観念できる任意の粒子を含み、中性子線やその他の放射線も含む。電子や光子のように、量子力学的に波動性と粒子性を兼ね備えるものも、これらの波又は粒子のいずれか又は両方に含まれる。同様に、古典力学的に波動や粒子と捉えられるものも、これらの波又は粒子のいずれか又は両方に含まれる。これらの波又は粒子は、それらの人為的な放出源となる装置（照射装置）により生成されるものであっても、また、例えば宇宙線などの自然由来で生成されるものや、例えば対象物の各部から放射される放射線などの対象物の性質として生じるもの（本出願書面では「照射装置によらずに生じた波又は粒子」と呼ぶ）であっても本開示の実施に用いることができる。
　１８０°の角度範囲をスキャンする場合、０°のスキャン角度の撮像画像と１８０°のスキャン角度の撮像画像は同一となるため、慣用的に「１８０°のスキャン角度範囲」という表現が、スキャン開始角度０°からスキャン終了角度（１８０°－１８０°／投影数）のスキャン角度範囲を意味するものとして用いられている。そこで、慣用的な表現での「１８０°のスキャン角度範囲」よりも小さいスキャン角度範囲を、本明細書及び特許請求の範囲においては、「（１８０°－１８０°／投影数）未満のスキャン角度範囲」と表現している。
　上記構成を有する本発明によれば、より小さな計算資源で、離散ラドン逆変換の代数的厳密解に裏付けられた再構成画像を生成することができるシステム、受信装置、送信装置、方法、プログラム、及びプログラムを記録した記録媒体を提供することができる。
　また、上記構成を有する本発明によれば、（１８０°－１８０°／投影数）未満のスキャン角度範囲で、より小さな計算資源で、離散ラドン逆変換の代数的厳密解に裏付けられた再構成画像を生成することができるシステム、方法、プログラム、及びプログラムを記録した記録媒体を提供することができる。
本発明の実施形態の一例の断層撮像システムおいて、画像が取得される対象物と検出装置の関係を示す斜視図である。 本発明の実施形態の一例の断層撮像システムにおいて、画像が取得される対象物の切断面における平面配置を含む概略構成を説明するための模式図であり、また、本発明の実施形態に係る断層撮像 システムの全体構成を示す図である。 従来の離散ラドン逆変換の代数的手法の原理を説明する説明図である。 検出素子数を解像度よりも大きくした場合の一例の断層画像の画素、波又は粒子のビームの照射、検出素子の位置の関係を示す図である。 投影数を解像度よりも大きくした場合の一例の断層画像の画素、波又は粒子のビームの照射、検出素子の位置の関係を示す図である。 検出素子数及び投影数を解像度よりも大きくした場合の一例の断層画像の画素、波又は粒子のビームの照射、検出素子の位置の関係を示す図である。 １６までの解像度についての検出素子数及び投影数の解像度に対する追加数に対する、厳密解が存在する解像度とスキャン角度範囲を示す図である。 ６４までの解像度についての検出素子数及び投影数の解像度に対する追加数に対する、厳密解が存在する解像度とスキャン角度範囲を示す図である。 本発明の実施形態の行列デシメーション処理における効率的な行ベクトルの選択方法の一例を示す図である。 本発明の実施形態に係る動作パラメータ算出装置の機能構成を示すブロック図である。 本発明の実施形態に係る離散ラドン逆変換行列生成装置、再構成画像生成装置の機能構成を示すブロック図である。 本発明の実施形態に係る動作パラメータ算出装置３０のハードウエア構成を示す図である。 本発明の実施形態に係る動作パラメータ算出処理の一例のフローチャートである。 本発明の実施形態に係る動作パラメータ算出処理の一例のフローチャート 本発明の実施形態に係る離散ラドン逆変換行列の生成処理の一例のフローチャートである。 本発明の実施形態に係る離散ラドン逆変換行列の生成処理の一例のフローチャートである。 本発明の実施形態に係る再構成画像生成処理の一例のフローチャートである。 本発明の実施形態に係る再構成画像生成処理の一例のフローチャートである。 検証のために用いた検証用データ（８ｂｉｔ数値ファントム）を示す図である。 検出素子数を解像度よりも大きくした場合のサイノグラムの例を示す図である。 検出素子数を解像度よりも大きくした場合の実施例において得られた再構成画像である。 図１５Ａのサイノグラムに基づいて、従来のＦＢＰ法においてＳｈｅｐｐ－Ｌｏｇａｎフィルタを用いて得られた再構成画像である。 図１５Ａのサイノグラムに基づいて、従来の代数再構成法として逐次近似計算（ＭＬ－ＥＭ法）を１０ 回行って得られた再構成画像である。 図１５Ａのサイノグラムに基づいて、基準としてＦＢＰ法においてフィルタを用いず得られた再構成画像である。 投影数を解像度よりも大きくした場合のサイノグラムの例を示す図である。 投影数を解像度よりも大きくした場合の実施例において得られた再構成画像である。 図１６Ａのサイノグラムに基づいて、従来のＦＢＰ法においてＳｈｅｐｐ－Ｌｏｇａｎフィルタを用いて得られた再構成画像である。 図１６Ａのサイノグラムに基づいて、従来の代数再構成法として逐次近似計算（ＭＬ－ＥＭ法）を１０ 回行って得られた再構成画像である。 図１６Ａのサイノグラムに基づいて、基準としてＦＢＰ法においてフィルタを用いず得られた再構成画像である。 検出素子数及び投影数を解像度よりも大きくした場合のサイノグラムの例を示す図である。 検出素子数及び投影数を解像度よりも大きくした場合の実施例において得られた再構成画像である。 図１７Ａのサイノグラムに基づいて、従来のＦＢＰ法においてＳｈｅｐｐ－Ｌｏｇａｎフィルタを用いて得られた再構成画像である。 図１７Ａのサイノグラムに基づいて、従来の代数再構成法として逐次近似計算（ＭＬ－ＥＭ法）を１０ 回行って得られた再構成画像である。 図１７Ａのサイノグラムに基づいて、基準としてＦＢＰ法においてフィルタを用いず得られた再構成画像である。
（原理）
　まず、本開示の断層撮像の原理について、図１Ａ～図７を参照して説明する。本実施形態における例示の幾何学的配置を説明する。図１Ａは、本実施形態の一例の断層撮像システム１００において、画像が取得される対象物２と検出装置２０の関係を示す斜視図である。また、図１Ｂは、本実施形態の一例の断層撮像システム１００において、画像が取得される対象物２の切断面における平面配置を含む概略構成を説明するための模式図である。この空間の一平面に、断層画像を構成するＮ画素×Ｎ画素を持つ画素の群が定義され、この整数Ｎが解像度となる。従来の手法では、この整数Ｎは検出装置２０の検出素子２２の数と一致させている。波又は粒子のビーム４は、その平面に含まれており、照射装置１０から照射される。検出装置２０と照射装置１０は、相対配置が固定されたまま対象物２に対して回転することができる。このような制御を含む装置の動作のために、断層撮像システム１００は、後述の制御部５０１を備えている。この回転は、上記画素の群や対象物２に対する相対的なものである。このため必ずしも検出装置２０と照射装置１０のみが回転されるとは限らず、検出装置２０と照射装置１０が設置床面（図示しない）に対し固定されていて対象物２が、画素の群とともに回転される場合もある。
　照射装置１０により波又は粒子のビーム４を照射し検出装置２０によりそのビーム４を対象物２の反対側で検出する動作は、従来の手法では検出方向θのスキャン範囲を画素数Ｎによって分割した刻みで検出方向についてサンプリングされる。この検出方向についてのサンプリング数を以下では「投影数」とも呼ぶ。図１Ｂのようなビーム４が平行なものでは、検出方向θのスキャン範囲は通常０°～１８０°とされ、すなわち、スキャン角度範囲φが１８０°とされ、従来の手法では検出装置２０からの検出信号が（１８０／Ｎ）°の刻みでサンプリングされる。なお、本実施形態を含めて、断層撮像システム１００の検出方向θの典型的なスキャン動作は、一定速度で増加又は減少させる定速回転をしながら、検出装置２０の値を取得するタイミングを電気的に制御することによりサンプリングが行われる。ただし、本実施形態はこのようなものに限定されるものではない。
　図２は、従来の離散ラドン逆変換の代数的手法の原理を説明する説明図である。この図は波又は粒子を照射し、透過後の強度データを用いて断層画像を再構築する処理の本質を示す初等的なものである。図示される縦２×横２に並ぶ正方形領域は、解像度Ｎ＝２に対応させてそれぞれが断層画像の画素を意図したもので、ａ1～ａ4の４つの値を持ちうる。これらの値は、各画素の位置での対象物２によるビームの減衰の程度を示し、決定されるべき値である。複数の画素を通り、解像度Ｎ＝２に対応させて１８０°／Ｎつまり９０°向きを変えて描かれた破線矢印は波又は粒子のビームの照射と検出素子の位置に対応している。矢印が指し示す先に記す数値は、検出装置の検出素子がその矢印の先に位置している際に得られる値（例えば減衰量）を意味しており、ここでは通った画素の値の合計が数値で示されている。つまり、この矢印が波又は粒子のビームを表し、矢印が通る画素と矢印の向きで照射及び検出（以下「照射」とのみ記す）を示している。矢印の先の合計値は、その矢印の先に検出装置の検出素子がある場合の、各検出素子で検出された減衰成分に対応している。矢印の指し示す値つまり検出装置から得られる値のみによってａ1～ａ4の４つの値を決定することが、検出装置の検出素子の示す値に基づいて離散ラドン逆変換を実行することにより断層画像を再構成する処理に対応する。
　図２に示された例では、ａ1～ａ4のうち矢印それぞれがつないでいるものが線形和を与えるため、次の連立方程式に書き直せる。

Figure JPOXMLDOC01-appb-I000001
ここで、ｘ線ＣＴの場合について例示すれば、ｐ1～ｐ4は図２のｘ線の各画素及び各検出方向での測定値であり、ａ1～ａ4はｘ線吸収係数である。この連立方程式は一意に値を定めることはできない。ここで、

Figure JPOXMLDOC01-appb-I000002
のように任意の数Ａを用いてシフトさせた新たなａ′1～ａ′4を代入しても式（１）は不変である。つまり、式（１）は元（決定すべき変数）の数が４であるのに対し、矢印の示す検出によって導き出せる連立方程式の実質的な数は３である。この場合、ａ１～ａ４の４つの値はいずれも不定（indeterminate）となり値を特定できない。
　この関係を別の数学的関係で整理すると、上述式（１）は、
Figure JPOXMLDOC01-appb-M000003
 のように、より一般に行列の形式に表現される。ｗnmはＮ2行Ｎ2列の正方行列である行列Ｗ（「システム行列」と呼ぶ）の要素である。一般に、行列Ｗが逆行列Ｗ-1を持つか否かは、階数の値による判定等の種々の数学的手法により判定することができる。すなわち、行列Ｗが逆行列を持つことは、Ｎ2行Ｎ2列の正方行列Ｗの階数（ｒａｎｋ（ｗnm））の値がＮ×Ｎに達することと同値である。行列Ｗが逆行列を持てば式（３）に対応する連立方程式を一意に解くことができる。上記式（３）の行列では、階数を計算するとｒａｎｋ（ｗnm）＝３＜４であり、逆行列を持つために必要な値に達していない。
　一方、行列Ｗが逆行列Ｗ-1を持つ場合は、ｐ1～ｐ4の値からａ1～ａ4の値を求めることができる。すなわち、測定値であるｐ1～ｐ4にこの逆行列を作用させるという代数的に厳密な手法のみによって、ａ1～ａ4を一意に決定することが可能となる。この違いが示唆するのは、照射の段階で適切な工夫をすれば、現実の断層画像の再構築処理のためのより一般の形をもつ多元連立一次方程式でも代数的手法によって解を求められる余地があること、すなわち、代数的手法によって不定ではない厳密解が得られる余地があることである。本発明者は、スキャン角度範囲が（１８０°－１８０°／投影数）未満であっても、厳密解を得ることができる場合があることを見出し、さらに、厳密解を得るための照射の段階での適切な工夫として、検出素子数、検出方向のサンプル数である投影数、スキャン角度範囲の３つの動作パラメータを調整すればよいことを見出した。また、この場合、検出素子数、投影数といったサンプル数を断層画像の解像度よりも大きくすると、厳密解の存在確率が大きくなり、またサンプル数を解像度よりも大きくすればするほど、厳密解の存在確率が上がると共に、厳密解が存在する範囲が大きくなり、スキャン角度範囲が小さい方向に厳密解が存在する範囲が拡がることを見出した。
　図１Ｂを再び参照する。ビーム４は上記座標にて全ての画素範囲をカバーする広がりを持ち、各検出素子２２に向かう部分はある幅τを持っている。当該検出素子２２に対するビーム４についてのその画素のもつ寄与は、典型的にはビーム４のうち、検出素子２２の一つに向かうものがもつ幅τと各画素（縦横それぞれδ）の重なりである。図１Ｂでは１つの画素についてこの重なりの様子をハッチで示している。Ｎ画素×Ｎ画素の全ての画素について、検出方向θ別、検出素子２２の別に寄与が定まる。この寄与を数値として保持しているのがシステム行列Ｗである。システム行列Ｗを用いると、照射と画像再構成のうち、照射を表す離散ラドン変換は、行列の積により、
　Ｘ′＝　Ｗ　Ｘ　　　　　　（４）
と表現される。ここで、システム行列Ｗは、総画素数に応じた空間についてのサンプリング点の並びを行方向に、検出方向の刻みと検出素子とに応じた照射及び検出についてのサンプリング点の並びを列方向にとるのが通常である。従来の手法を図１Ｂの画素数などの設定に適用した場合には、システム行列Ｗが行方向にも列方向にもＮ×Ｎ個の要素を持つ正方行列となりその要素の総数はＮ4となる。これは、画素がＮ画素×Ｎ画素、検出方向θのスキャン範囲がＮサンプル、検出素子２２の数がＮ素子であるためである。なお、画質を高めるためにＮの値は時代とともに増大しており、典型例では１０２４すなわち１０3程度とされ、従来のものでは２０００程度が最大値であるが、本実施形態の場合にもＮの値が得に限定されるものではない。Ｎが１０3程度のものではシステム行列は行方向列方向ともに１０6程度の規模、要素数は１０12程度となる。
　式（４）のＸは各画素の位置での減衰（吸収、散乱も含む）を表す像である減衰像を列ベクトルとして表現したものである。実際の空間での減衰像は図１Ｂに示したＮ画素×Ｎ画素の二次元画像であるが、システム行列Ｗとの演算のために、列方向を総画素数Ｎ×Ｎになるようにベクトル化したものがＸである。式（４）において減衰像Ｘにシステム行列Ｗが作用して得られる左辺の列ベクトルＸ′は、サイノグラムと呼ばれる像に対応している。システム行列Ｗにて列方向を検出方向のサンプリング点と検出素子の並びとしたことに対応して、列ベクトルＸ′も検出方向のサンプリング点と検出素子それぞれについて得られる検出強度を示している。なお、サイノグラムは検出素子の軸と検出方向刻みの軸との２次元に列ベクトルＸ′の要素を並べ替えたものである。以上のように、式（４）は照射と検出を表現したものである。図２に示した図では、列ベクトルＸはａ1～ａ4の各値に、また列ベクトルＸ′は各矢印の先に付した値にそれぞれ対応し、システム行列Ｗは、各矢印と各画素の対応関係を反映するものといえる。
　断層画像の再構成とは、式（４）が実現していることを前提として、システム行列Ｗとサイノグラムの列ベクトルＸ′から列ベクトルＸを算出し、列ベクトルＸから減衰像を得る処理である。つまり断層画像の再構成は形式的に次式により表現される。
　Ｘ＝　Ｗ-1　Ｘ′　　　　　　（５）
ここで、Ｗ-1はＷの逆行列であり、式（５）は式（４）の両辺にＷ-1を作用させて左右辺を入れ替えたものである。実際には、正方行列であるＷが常に逆行列Ｗ-1をもつとは限らない。図２でａ1～ａ4の各値が定まらないのは、それらの照射と検出の仕方ではＷ-1が存在しないことに対応する。なお上述した図１Ｂに基づく説明は平行ビームを照射するものであるが、本実施形態はファンビームやコーンビームといった平行ではないビームを利用するような照射装置１０や検出装置２０を採用する場合にも当業者には明らかなわずかな変更によって適用することができる。また、図１Ｂの画素数などの設定以外においても、本実施形態の手法を適用することができる。なお、式（４）、（５）を、数学的に等価な他の形式により表現することもでき、たとえば、システム行列Ｗのために、上述したものとは異なり、総画素数に応じた空間についてのサンプリング点の並びを列方向に、検出方向の刻みと検出素子とに応じた照射及び検出についてのサンプリング点の並びを行方向にとる。そしてそれに対応させて、式（４）のＸを各画素の位置に応じて行ベクトルとして表現することもできる。その場合式（４）、（５）に相当する式は、それぞれの各辺に転置操作を施したものとなる。この転置操作では行列やベクトルの行と列を入れ替え、ベクトルに作用させる行列が左からではなく右から作用させる順序に変わるという表現上の影響がある。ベクトルが「列又は行ベクトル」のように択一的に記載され、対応する行列の要素の並びが「列方向又は行方向」のように択一的に記載されているときは、これらの組み合わせの選択は、それぞれの択一的記載における同順のものを組み合わせて選ぶこととする。この表現の選択は本出願の本質に影響を与えないため、本願においては特に断りの無い場合には式（４）のようにサイノグラムをベクトル化したＸ等を列ベクトルにより表現し、システム行列等の行列は左から作用させる表記にもとづいて説明する。
　本発明者は、スキャン角度範囲が（１８０°－１８０°／投影数）未満であっても、厳密解を得ることができる場合があることを見出し、さらに、上で述べた、厳密解を得るための照射の段階での適切な工夫として、検出素子数、検出方向のサンプル数である投影数、スキャン角度範囲の３つの動作パラメータを調整すればよいことを見出した。また、この場合、検出素子数、投影数といったサンプル数を断層画像の解像度よりも大きくすると、厳密解の存在確率が大きくなり、またサンプル数を解像度よりも大きくすればするほど、厳密解の存在確率が上がると共に、厳密解が存在する範囲が大きくなり、スキャン角度範囲が小さい方向に厳密解が存在する範囲が拡がることを見出した。以下、その例として、検出素子数を解像度よりも大きくした場合、投影数を解像度よりも大きくした場合、検出素子数及び投影数を解像度よりも大きくした場合について説明する。
（１）検出素子数を解像度よりも大きくした場合
　例として、解像度が２、すなわち断層画像が２画素×２画素であり、検出素子数が４、投影数が２、スキャン角度範囲が２２°の場合について説明する。
　この場合、幾何学的配置は、図３に示されるようなものとなるので、各検出素子により検出される値ｐ1～ｐ8は、各画素の位置での対象物２によるビームの減衰の程度ａ1～ａ4を用いて、式（６）のように記述される。

Figure JPOXMLDOC01-appb-I000004
そして、これを行列の形式で表現すると式（７）のようになり、システム行列Ｗは式（８）のようになる。

Figure JPOXMLDOC01-appb-I000005

Figure JPOXMLDOC01-appb-I000006
ここで、システム行列Ｗは、正方行列でないため、式（７）において逆行列を左からかけてａ1～ａ4を求めることができない。そこで、システム行列Ｗを構成する８個の行ベクトルの要素のうちの４個の行ベクトルの要素を間引くことによって、システム行列Ｗを正方化する。この正方化された行列を正方システム行列ＷSQと呼ぶ。
　この場合、８個の行ベクトルの要素から４個の行ベクトルの要素を間引くことによって、４個の行ベクトルの要素を選択することになるので、組み合わせは全部で8Ｃ4＝７０通りあるが、正方システム行列ＷSQが逆行列を持つ組み合わせは９通りである。そのうちの１つが、点線で囲まれた行ベクトルを選択した場合であり、この場合、得られた正方システム行列ＷSQに対して式（９）が成り立つ。
Figure JPOXMLDOC01-appb-I000007
そして、逆行列ＷSQ -1を左からかけることによって、式（１０）のように、ａ1～ａ4を求めることができ、断層画像が得られる。

Figure JPOXMLDOC01-appb-I000008
（２）投影数を解像度よりも大きくした場合
　例として、解像度が２、すなわち断層画像が２画素×２画素であり、検出素子数が２、投影数が６、スキャン角度範囲が６５．２°の場合について説明する。
　この場合、幾何学的配置は、図４に示されるようなものとなるので、各検出素子により検出される値ｐ1～ｐ8は、各画素の位置での対象物２によるビームの減衰の程度ａ1～ａ4を用いて、式（１１）のように記述される。

Figure JPOXMLDOC01-appb-I000009
そして、これを行列の形式で表現すると式（１２）のようになり、システム行列Ｗは式（１３）のようになる。
Figure JPOXMLDOC01-appb-I000010

Figure JPOXMLDOC01-appb-I000011
　この場合、１２個の行ベクトルの要素から８個の行ベクトルの要素を間引くことによって、４個の行ベクトルの要素を選択することになるので、組み合わせは全部で12Ｃ4＝４９５通りあるが、正方システム行列ＷSQが逆行列を持つ組み合わせは１６０通りである。そのうちの１つが、点線で囲まれた行ベクトルを選択した場合であり、この場合、得られた正方システム行列ＷSQに対して式（１４）が成り立つ。
Figure JPOXMLDOC01-appb-I000012
そして、逆行列ＷSQ -1を左からかけることによって、式（１５）のように、ａ1～ａ4を求めることができ、断層画像が得られる。

Figure JPOXMLDOC01-appb-I000013
（３）検出素子数及び投影数を解像度よりも大きくした場合
　例として、解像度が２、すなわち断層画像が２画素×２画素であり、検出素子数が４、投影数が４、スキャン角度範囲が５４°の場合について説明する。
　この場合、幾何学的配置は、図５に示されるようなものとなるので、各検出素子により検出される値ｐ1～ｐ16は、各画素の位置での対象物２によるビームの減衰の程度ａ1～ａ4を用いて、式（１６）のように記述される。

Figure JPOXMLDOC01-appb-I000014
そして、これを行列の形式で表現すると式（１７）のようになり、システム行列Ｗは式（１８）のようになる。
Figure JPOXMLDOC01-appb-I000015

Figure JPOXMLDOC01-appb-I000016
　この場合、１６個の行ベクトルから１２個の行ベクトルの要素を間引くことによって、４個の行ベクトルの要素を選択することになるので、組み合わせは全部で16Ｃ4＝１８２０通りあるが、正方システム行列ＷSQが逆行列を持つ組み合わせは４６５通りである。そのうちの１つが、点線で囲まれた行ベクトルを選択した場合であり、この場合、得られた正方システム行列ＷSQに対して式（１９）が成り立つ。
Figure JPOXMLDOC01-appb-I000017
そして、逆行列ＷSQ -1を左からかけることによって、式（２０）のように、ａ1～ａ4を求めることができ、断層画像が得られる。

Figure JPOXMLDOC01-appb-I000018
　以上は、検出素子数及び／又は投影数を解像度よりも大きくした場合であったが、上述のように、検出素子数及び投影数が解像度と同じであっても、（１８０°－１８０°／投影数）未満のスキャン角度範囲に対して厳密解を有する場合があることを本発明者は見出した。そこで、検出素子数及び投影数の解像度に対する追加数に対する、厳密解が存在する解像度とスキャン角度範囲を、数値計算によって調べた。結果を図６Ａ及び図６Ｂに示す。図６Ａは、１６までの解像度についての検出素子数及び投影数の解像度に対する追加数に対する、厳密解が存在する解像度とスキャン角度範囲を示す図であり、図６Ｂは、６４までの解像度についての検出素子数及び投影数の解像度に対する追加数に対する、厳密解が存在する解像度とスキャン角度範囲を示す図である。図６Ｂについては、検出素子数及び投影数の両方を追加した場合の図を省略している。図６Ａ及び図６Ｂは、後述の（計算機シミュレーションによる検証（実施例））と同様の計算実験条件において、追加検出素子数、追加投影数、スキャン角度範囲、解像度を変えて得られたものである。図６Ａ及び図６Ｂから分かるように、スキャン角度範囲が（１８０°－１８０°／投影数）未満で、検出素子数、投影数、スキャン角度範囲を適切に選択することによって、厳密解を得ることができる。また、検出素子数、投影数といったサンプル数を断層画像の解像度よりも大きくすると、厳密解の存在確率が大きくなり、またサンプル数を解像度よりも大きくすればするほど、厳密解の存在確率が上がると共に、厳密解が存在する範囲が大きくなり、スキャン角度範囲が小さい方向に厳密解が存在する範囲が拡がることが分かる。したがって、例えば、検出素子数と投影数を適切に選択すれば、所望の（１８０°－１８０°／投影数）未満のスキャン角度範囲で厳密解を得ることができる。例えば、上述の非特許文献６の例では、スキャン角度範囲を１６０°以下にする必要があるが、１６０°以下のスキャン角度範囲に対して、厳密解を得ることができる適切な検出素子数と投影数を探索すれば、探索された検出素子数と投影数で測定を行うことにより厳密解に基づく再構成画像を得ることができる。
　ここで、検出素子を追加する場合には、（１）検出素子のサイズを小さくして検出素子を追加する場合（典型的には、検出素子群の領域が変わらないように検出素子のサイズを小さくして検出素子を追加する場合）、（２）検出素子のサイズを変えずに検出素子を追加する場合、（３）検出素子のサイズを小さくして検出素子を追加する場合が考えられる。（１）～（３）のいずれの場合であっても、後述の数学的な処理に鑑みて、追加する検出素子数が同じであれば、逆行列をもつ正方システム行列を得ることにより厳密解を得る点での作用効果は同じである。したがって、図６Ａ及び図６Ｂの結果は、検出素子群と計測対象物の相対的な大小といった物理的な長さの情報によらず、システム行列の構造のみに着目した評価であるが、厳密解を得る条件の決定には充分である。
　また、本実施形態により得られる再構成画像は、厳密解に基づく品質の高い画像であるので、入力された断層撮像画像から再構成画像を予測する機械学習システムのための学習用データとして有用である。
　以上を一般化して説明すると、以下のようになる。
＜離散ラドン逆変換行列の生成処理＞
　断層撮像でＮ画素×Ｎ画素（Ｎは正の整数）の二次元画像を再構成する場合、従来の手法では、検出装置２０の検出素子２２の数はＮ個とされるが、本実施形態においては、（Ｎ＋ｊ）個（ｊは０以上の整数）とされる。検出方向θについてのサンプリングは、従来の手法では、上述したように、０～１８０°の範囲、すなわち１８０°のスキャン角度範囲において重複しないＮ方向となるように、（１８０／Ｎ）°刻みでサンプリングが行われる。これに対し本実施形態の検出方向のサンプリングは、（１８０°－１８０°／投影数）未満のスキャン角度範囲φにおいて重複しないＮ＋ｋ方向（ｋは０以上の整数）でなされ、例えば（φ／（Ｎ＋ｋ－１））°刻みで行われる。サンプリングの間隔（検出方向の刻み）は、φ／（Ｎ＋ｋ－１）°といった等間隔でなく、異なる間隔としてもよい。このように、Ｎ個に代えて（Ｎ＋ｊ）個の検出素子で、Ｎ方向に代えてＮ＋ｋ方向で検出方向のサンプリングを行う。
　この検出素子数（Ｎ＋ｊ）、投影数（Ｎ＋ｋ）、スキャン角度範囲φに基づいて、システム行列を決定する。ここで決定されるシステム行列は、検出方向と検出素子それぞれに対応させて列方向に（Ｎ＋ｋ）×（Ｎ＋ｊ）、断層画像の画素数に対応させて行方向にＮ×Ｎだけの要素をもつ。そして、逆行列を算出することができるようにするために、システム行列Ｗから行列の適切な位置にある（ｊ＋ｋ）×Ｎ＋ｊ×ｋ個の、Ｎ×Ｎ個の要素をもつ行ベクトルの要素を除去して、正方システム行列ＷSQを生成する。除去される行ベクトルの要素が適切であれば、正方システム行列ＷSQは逆行列を持つ。正方システム行列ＷSQからその逆行列ＷSQ -1である離散ラドン逆変換行列を算出する。ここで逆行列ＷSQ -1の算出に当たっては、例えば余因子法や掃き出し法など任意の適切な周知の手法を用いることができる。
＜再構成画像生成処理＞
　一方、対象物２からの断層画像を得るには、照射装置１０と各検出素子２２をもつ検出装置２０との間に、対象物２を配置する。照射装置１０からのビーム４を照射しながら、重複しないＮ＋ｋ方向で検出が行われる。典型的に、このビーム４は連続して照射し続けていて、対象物２と照射装置１０及び検出装置２０の向き（検出方向）はなめらかに動きつつ、検出装置２０の検出素子２２からの信号が電気的にサンプリングされる。重複しないＮ＋ｋ方向での照射及び検出の典型例は、（φ／（Ｎ＋ｋ－１））°刻みのタイミングでのサンプリングである。得られる信号は、必要に応じて機器の物理的特性について補正やスケールの換算などの処理と配列順序の整理などをすれば、サイノグラムを生成しうるデータに整理される。このサイノグラムを得る処理が、照射及び検出の処理であり、離散ラドン変換である。このサイノグラムからは、システム行列Ｗ列方向の並びに合わせるような順序で要素を並べ替えた列ベクトルを得ることができるが、このベクトルを本実施形態では第１ベクトルと呼ぶ。第１ベクトルＸ′は、サイノグラムの値を持つので（Ｎ＋ｋ）×（Ｎ＋ｊ）要素を持つ。
　正方システム行列ＷSQを生成した際に除去された行ベクトルの要素に対応する、第１ベクトルＸ′の要素を間引くベクトルデシメーション処理を行うことにより、正方システム行列ＷSQとの間で式（４）の関係を満たすことができるＮ×Ｎ要素をもつ列ベクトルである第２ベクトルＸSQ′を得ることができる。第２ベクトルＸSQ′は正方システム行列ＷSQとの間で式（４）の関係を満たすため、正方システム行列ＷSQの逆行列である離散ラドン逆変換行列ＷSQ -1を用いて式（５）の関係が成り立ち、離散ラドン逆変換処理によって再構成像を与える第３ベクトルＸSQが得られる。この第３ベクトルＸSQの並びをベクトル化の逆操作によってＮ画素×Ｎ画素の画像に戻すデベクタライズ処理により再構成像を画像化することができる。これで１スライス分の断層画像が再構築される。
　本実施形態の処理は、上述の離散ラドン逆変換行列生成処理と再構成画像生成処理の２系統に分かれている。通常再構成像を撮像するユーザは、照射及び検出から再構成画像生成に至る上述の再構成画像生成処理を実行する。つまりユーザは、上述の動作パラメータに対応するシステム行列生成から離散ラドン逆変換行列生成に至る離散ラドン逆変換行列の生成処理を必ずしも実行する必要はない。離散ラドン逆変換行列生成処理が必要となるのは、断層撮像システムを作製したり、設置やメンテナンスしたりする場面で、いわば機器メーカーが対応するような処理である。間引くべき行ベクトルの要素及び離散ラドン逆変換行列ＷSQ -1が得られていれば、再構成像を撮像しようとするユーザは自らの測定のためにはそれらを利用するだけでよい。
　特筆すべきは、再構成画像生成処理が行列やベクトルの順序の入れ替えや行列積の演算であり、例えば逐次近似といった計算処理量が問題となる処理がその範囲に含まれていないことである。本実施形態において、ユーザの通常の撮像において利用される再構成画像生成処理は、離散ラドン逆変換行列生成処理により準備された、間引くべき行ベクトルの要素及び離散ラドン逆変換行列ＷSQ -1を使うだけで十分に軽い計算負担で実行可能であり、この点は、本実施形態に高い実用性をもたらすものといえる。断層画像を取得するという最も頻度の高い処理の計算負担がごく少ないためである。特に、三次元ボリューム像を得るために多数のスライスの撮像を繰り返す用途には本実施形態の手法が極めて有利である。その繰り返しに離散ラドン逆変換行列生成処理が不要なことに加え、再構成画像生成処理は、ＧＰＧＰＵ（General-purpose computing on graphics processing unit）などのメニーコアタイプのプロセッサーでの計算処理に適するため、今後の処理速度の向上も期待できるからである。
　また、幾何学的な側面を反映するシステム行列Ｗや正方システム行列ＷSQ、システム行列ＷSQの逆行列ＷSQ -1が、離散ラドン逆変換行列生成処理のみに含まれており、測定のための再構成画像生成処理と分離されていることから、代数的手法に反復を利用したときの収束の問題と測定のためのノイズの側面とが本質的に分離されている。このため、本実施形態では、ノイズとは無関係な近似解に起因した現象が生じにくく、画像処理で一般に採用される各種のノイズ軽減処理が目的通りに効果を発揮しやすいという利点もある。
＜間引くべき行ベクトルの要素の具体的な決定手法＞
　上述したように、離散ラドン逆変換行列ＷSQ -1が正方システム行列ＷSQの逆行列として求まることは、間引くべき行ベクトルの要素の適切さに依存している。その手法は、結果として逆行列を持つ正方システム行列を生成することができれば、限定はされるものではないが、ここでは、階数での判定を利用する乱択アルゴリズム（randomized algorithm）について以下説明する。
　システム行列Ｗの行の総数は（Ｎ＋ｋ）×（Ｎ＋ｊ）である。このうち、行列を正方行列にするために残す行数はＮ×Ｎであり、間引く行ベクトルの総数は（ｊ＋ｋ）×Ｎ＋ｊ×ｋである。行を特定するインデックス（行番号）の整数の範囲［１，（Ｎ＋ｋ）×（Ｎ＋ｊ）］のうちから、残すＮ×Ｎ個か間引く（ｊ＋ｋ）×Ｎ＋ｊ×ｋ個かいずれかの数列を決定する。この決定のために乱択手法を用いる。つまり、残すＮ×Ｎ個か間引く（ｊ＋ｋ）×Ｎ＋ｊ×ｋ個かのいずれかを数値生成された乱数でランダムに決定する。これにより、最終的に間引き数列｛ｍ｝を決定する。間引き数列｛ｍ｝が決まれば、対応した間引き処理によってシステム行列Ｗから正方システム行列ＷSQを得ることができる。次に、求めた正方システム行列ＷSQの階数（rank）を算出する。算出した階数は、正方システム行列ＷSQが逆行列を持つ場合はＮ×Ｎと同じ値となり、逆行列を有しない場合はＮ×Ｎよりも小さい値となる。このため、階数の値によって、正方システム行列ＷSQが逆行列を持ちうるかどうかを判定することができる。階数がＮ×Ｎと一致しなければ、再び別の乱数に基づいてランダムアルゴリズムによって間引き数列｛ｍ｝の選択し、階数を算出して逆行列を有するか否かの判定をすればよい。このような乱択アルゴリズムを採用すれば、実用的な速度で行列デシメーション処理を実行することができる。
　ここで、行ベクトルの異なるデシメーションによって、同じ行ベクトルの組み合わせが選択された場合は、正方システム行列ＷSQは同一になる。本発明者はこの点に着目して、以下のような効率的な行ベクトルの選択方法を着想した。
　以下、上記の「（３）検出素子数及び投影数を解像度よりも大きくした場合」のデシメーションを例として説明する。
　図７を参照して、ｐ1に対応する行ベクトルは、(１０１０)であるが、ｐ2、ｐ5、ｐ9に対応する行ベクトルが同じ行ベクトル(１０１０)である。そこで、ｐ2、ｐ5、ｐ9に対応する行ベクトルを間引く。
　次に、ｐ3に対応する行ベクトルは、(０１０１)であるが、ｐ8、ｐ12に対応する行ベクトルが同じ行ベクトル(０１０１)である。そこで、ｐ8、ｐ12に対応する行ベクトルを間引く。
　続いて、ｐ6に対応する行ベクトルは、(１１１０)であるが、ｐ10、ｐ14に対応する行ベクトルが同じ行ベクトル(１１１０)である。そこで、ｐ10、ｐ14に対応する行ベクトルを間引く。
　さらに、ｐ7に対応する行ベクトルは、(０１１１)であるが、ｐ11、ｐ15に対応する行ベクトルが同じ行ベクトル(０１１１)である。そこで、ｐ11、ｐ15に対応する行ベクトルを間引く。
　以上の間引き処理によって、残った行ベクトルは、互いに異なる行ベクトルとなる。残った行ベクトルは６個であるから、正方システム行列ＷSQを生成するための組み合わせの数は、16Ｃ4＝１８２０通りから6Ｃ4＝１５通りまで減らすことができる。このようにして、重複する行ベクトルをまず間引く（除去する）ことによって、選択する候補となる行ベクトルの数を減らすことができ、組み合わせの数を減らすことができるので、逆行列の有無の判定のための計算量を減らすことができる。
＜間引かれた重複する行ベクトルを利用した画質向上方法＞
　また、行ベクトルの異なるデシメーションによって、同じ行ベクトルの組み合わせが選択された場合は、正方システム行列ＷSQは同一になる点に着目して、本発明者は、以下のような効率的な再構成画像の画質向上方法を着想した。
　例えば、上述のようにｐ2、ｐ3、ｐ10、ｐ11に対応する行ベクトルの要素から構成される正方システム行列ＷSQは逆行列を有するが、上述のように、ｐ1に対応する行ベクトルの要素は、ｐ2に対応する行ベクトルの要素と同じ行ベクトルの要素であるから、ｐ1、ｐ3、ｐ10、ｐ11から構成される正方システム行列ＷSQは、ｐ2、ｐ3、ｐ10、ｐ11に対応する行ベクトルの要素から構成される正方システム行列ＷSQと同一であり、逆行列を有する。一方、ｐ1とｐ2は、異なる検出素子によって独立に測定されているから、ｐ2、ｐ3、ｐ10、ｐ11から求められる再構成画像とｐ1、ｐ3、ｐ10、ｐ11から求められる再構成画像とを重ね合わせることにより、再構成画像のＳＮ比が向上し、画質が向上する。なぜならば、一般に、測定回数が多いほど測定精度は向上し、ｎ回測定を行うと、ＳＮ比はｎ1/2倍になるといわれているからである。同様に、上述のように、ｐ5、ｐ9に対応する行ベクトルの要素は、ｐ2に対応する行ベクトルの要素と同じ行ベクトルの要素であるから、結局、ｐ2、ｐ3、ｐ10、ｐ11から求められる再構成画像、ｐ1、ｐ3、ｐ10、ｐ11から求められる再構成画像、ｐ5、ｐ3、ｐ10、ｐ11から求められる再構成画像、ｐ9、ｐ3、ｐ10、ｐ11から求められる再構成画像とを重ね合わせることにより、再構成画像のＳＮ比がさらに向上し、画質が向上する。
　このようにして、（１）逆行列を持つ組み合わせを早く発見するためのデシメーション法と、（２）そのデシメーション法で間引いた行ベクトルの要素を含む逆行列の組み合わせに基づく再構成画像を用いて精度を上げる手法を組み合わせた手法により、最も計算量が少なく、効率的に、精度の高い再構成画像を得ることができる。
（本発明の具体的な実施形態）
　以下、本発明の具体的な実施形態について図面を参照して説明する。上記の（原理）の項と重複する説明は省略する。
　図１Ｂは、本発明の実施形態に係る断層撮像システムの全体構成を示す図である。断層撮像システム１は、照射装置１０、検出装置２０、動作パラメータ算出装置３０、離散ラドン逆変換行列生成装置４０、再構成画像生成装置５０を備える。照射装置１０、検出装置２０、動作パラメータ算出装置３０、離散ラドン逆変換行列生成装置４０、再構成画像生成装置５０は、１つの物理的な装置として構成されていなくてもよく、照射装置１０及び検出装置２０、動作パラメータ算出装置３０、離散ラドン逆変換行列生成装置４０、再構成画像生成装置５０は、その一部又は全部のそれぞれが別個の物理的な装置として構成されてもよい。その場合、別個な物理的装置として構成された装置は互いに有線／無線のコンピュータネットワークを介して接続されてもよい。照射装置１０、検出装置２０、動作パラメータ算出装置３０、離散ラドン逆変換行列生成装置４０、再構成画像生成装置５０の各々は、１つの物理的な装置として構成される必要はなく、複数の物理的な装置から構成されてもよい。
　図８は、本発明の実施形態に係る動作パラメータ算出装置の機能構成を示すブロック図である。
　動作パラメータ算出装置３０は、動作パラメータ解釈部３０１、動作パラメータ受付部３０３、システム行列決定部３０５、正方システム行列生成部３０７、逆行列判定部３０９、動作パラメータ算出結果出力部３１１を備える。
　動作パラメータ解釈部３０１は、ユーザにより入力されたデータが動作パラメータの値であった場合、その値を指定し、入力されたデータが動作パラメータの値の組及び／又は値の範囲であった場合、その中の１つの値を指定する。
　動作パラメータ受付部３０３は、検出素子数、投影数、及び（１８０°－１８０°／投影数）未満のスキャン角度範囲のうちの指定された１つ又は２つの動作パラメータの値を受け付ける。
　システム行列決定部３０５は、前記指定された１つ又は２つの動作パラメータの値、及び前記指定された１つ又は２つの動作パラメータの値に対応する動作パラメータ以外の動作パラメータについては所定の値に基づいて、システム行列を決定する。
　正方システム行列生成部３０７は、ｊ＝ｋ＝０以外の場合に、決定された前記システム行列から、所定の計（ｊ＋ｋ）×Ｎ＋ｊ×ｋ個の行ベクトル又は列ベクトルの要素を除去することにより正方システム行列を生成する。
　逆行列判定部３０９は、決定された正方システム行列について、逆行列が存在するか否かを判定する。
　動作パラメータ算出結果出力部３１１は、指定された１つ又は２つの動作パラメータの値に対応する動作パラメータ以外の動作パラメータの所定の値を表示、音声出力等の形態で出力する。
　図９は、本発明の実施形態に係る離散ラドン逆変換行列生成装置、再構成画像生成装置の機能構成を示すブロック図である。
　離散ラドン逆変換行列生成装置４０は、動作パラメータ受付部４０３、システム行列決定部４０５、正方システム行列生成部４０７、逆行列判定部４０９、離散ラドン逆変換行列生成部４１１、記憶部４１３を備える。
　動作パラメータ受付部４０３は、指定された３つの動作パラメータの値を受け付ける。
　システム行列決定部４０５は、指定された３つの動作パラメータの値に基づいて、システム行列を決定する。
　正方システム行列生成部４０７は、ｊ＝ｋ＝０以外の場合に、決定された前記システム行列から、所定の計（ｊ＋ｋ）×Ｎ＋ｊ×ｋ個の行ベクトル又は列ベクトルの要素を除去することにより正方システム行列を生成する。
　逆行列判定部４０９は、決定された正方システム行列について、逆行列が存在するか否かを判定する。
　離散ラドン逆変換行列生成部４１１は、正方システム行列ＷSQの逆行列である離散ラドン逆変換行列ＷSQ -1を生成する。
　記憶部４１３は、正方システム行列ＷSQを構成する行ベクトルの要素のうちで、ステップＳ２０７において除去された行ベクトルの要素と同じものが存在すると判定した場合、その正方システム行列ＷSQを構成する行ベクトルの要素及び除去された同じ行ベクトルの要素に対応する行を特定するインデックス（行番号）を記憶する。
　再構成画像生成装置５０は、制御部５０１、サイノグラム生成部５０３、第１ベクトル生成部５０５、ベクトルデシメーション部５０７、第３ベクトル生成部５０９、デベクタライズ部５１１、別測定値第２ベクトル生成部５１３、合成再構成画像生成部５１５、記憶部５１７を備える。
　制御部５０１は、照射装置１０及び検出装置２０を制御して、照射装置１０からのビーム４を対象物２に照射しながら、重複しないＮ＋ｊ方向で検出装置２０によって検出を行う。
　サイノグラム生成部５０３は、検出装置２０の各検出素子２２から得られる信号について、必要に応じて機器の物理的特性について補正やスケールの換算などの処理と配列順序の整理などを行い、サイノグラムを生成する。
　第１ベクトル生成部５０５は、生成されたサイノグラムから、システム行列Ｗの列方向の並びに合わせるような順序で要素を並べ替えた列ベクトルである第１ベクトルＸ′を生成する。
　ベクトルデシメーション部５０７は、正方システム行列ＷSQを生成した際に除去された行ベクトルの要素に対応する、第１ベクトルＸ′の要素を間引き、第２ベクトルＸSQ′を生成する。
　第３ベクトル生成部５０９は、離散ラドン逆変換行列生成部４１１により生成された離散ラドン逆変換行列ＷSQ -1と第２ベクトルＸSQ′を用いて、式（５）に基づいて、離散ラドン逆変換処理によって再構成像を与える第３ベクトルＸSQを生成する。
　デベクタライズ部５１１は、第３ベクトルＸSQの並びをベクトル化の逆操作によってＮ画素×Ｎ画素の画像に戻すデベクタライズ処理により再構成画像を生成する。
　別測定値第２ベクトル生成部５１３は、正方システム行列生成部４０７によって記憶部４１３に記憶されたインデックス（行番号）に基づいて、正方システム行列生成部４０７によって、正方システム行列ＷSQを構成する行ベクトルの要素のうちで、正方システム行列生成部４０７により除去された行ベクトルの要素と同じものが存在すると判定された場合の正方システム行列ＷSQを構成する行ベクトルのインデックス（行番号）に対応する、第２ベクトルＸSQ′の要素の一部又は全部を、除去された同じ行ベクトルの要素に対応する第２ベクトルＸSQ′の要素に置き換えて、別測定値第２ベクトルを生成する。
　合成再構成画像生成部５１５は、生成された複数の再構成画像を重ね合わせて、画質が向上した合成再構成画像を生成する。
　記憶部５１７は、ステップＳ３１９で生成された再構成画像を記憶する。
　図１０は、本発明の実施形態に係る動作パラメータ算出装置３０のハードウエア構成を示す図である。動作パラメータ算出装置３０は、ＣＰＵ３０ａ、ＲＡＭ３０ｂ、ＲＯＭ３０ｃ、外部メモリ３０ｄ、入力部３０ｅ、出力部３０ｆ、通信部３０ｇを含む。ＲＡＭ３０ｂ、ＲＯＭ３０ｃ、外部メモリ３０ｄ、入力部３０ｅ、出力部３０ｆ、通信部３０ｇは、システムバス３０ｈを介して、ＣＰＵ３０ａに接続されている。
　ＣＰＵ３０ａは、システムバス３０ｈに接続される各デバイスを統括的に制御する。
　ＲＯＭ３０ｃや外部メモリ３０ｄには、ＣＰＵ３０ａの制御プログラムであるＢＩＯＳやＯＳ、コンピュータが実行する機能を実現するために必要な各種プログラムやデータ等が記憶されている。
　ＲＡＭ３０ｂは、ＣＰＵの主メモリや作業領域等として機能する。ＣＰＵ３０ａは、処理の実行に際して必要なプログラム等をＲＯＭ３０ｃや外部メモリ３０ｄからＲＡＭ３０ｂにロードして、ロードしたプログラムを実行することで各種動作を実現する。
　外部メモリ３０ｄは、例えば、フラッシュメモリ、ハードディスク、ＤＶＤ－ＲＡＭ、ＵＳＢメモリ等から構成される。
　入力部３０ｅは、ユーザ等からの操作指示等を受け付ける。入力部３０ｅは、例えば、入力ボタン、キーボード、ポインティングデバイス、ワイヤレスリモコン、マイクロフォン等の入力デバイスから構成される。
　出力部３０ｆは、ＣＰＵ３０ａで処理されるデータや、ＲＡＭ３０ｂ、ＲＯＭ３０ｃや外部メモリ３０ｄに記憶されるデータを出力する。出力部３０ｆは、例えば、ＣＲＴディスプレイ、ＬＣＤ、有機ＥＬパネル、プリンタ、スピーカ等の出力デバイスから構成される。
　通信部３０ｇは、ネットワークを介して又は直接、外部機器と接続・通信するためのインタフェースである。通信部３０ｇは、例えばシリアルインタフェース、ＬＡＮインタフェース等のインタフェースから構成される。
　離散ラドン逆変換行列生成装置４０及び再構成画像生成装置５０のハードウエア構成も同様である。
　図８に示される動作パラメータ算出装置３０の各部は、ＲＯＭや外部メモリに記憶された各種プログラムが、ＣＰＵ、ＲＡＭ、ＲＯＭ、外部メモリ、入力部、出力部、通信部等を資源として使用することで実現される。
　以上のシステム構成を前提に、本発明の実施形態に係る断層撮像システムの断層撮像処理の例を、図１～図１７Ｅ等を参照して、以下に説明する。上記の（原理）の項と重複する説明は省略する。
＜動作パラメータ算出処理＞
　図１１Ａ及び図１１Ｂは、本発明の実施形態に係る動作パラメータ算出処理の一例のフローチャートである。
　ユーザが、検出素子数、投影数、及び（１８０°－１８０°／投影数）未満のスキャン角度範囲のうちの１つ又は２つの動作パラメータの値、値の組、又は値の範囲を動作パラメータ算出装置３０の入力部３０ｅを介して入力する（Ｓ１０１）。動作パラメータ算出装置３０の動作パラメータ解釈部３０１は、入力されたデータが動作パラメータの値であった場合、その値を指定し、入力されたデータが動作パラメータの値の組及び／又は値の範囲であった場合、その中の１つの値を指定する（Ｓ１０３）。入力されたデータが動作パラメータの値の組及び／又は値の範囲であった場合、入力された動作パラメータの値の組及び／又は値の範囲が尽きるまで以下の処理を繰り返すが、以下、簡単のため、入力されたデータが動作パラメータの値であった場合について説明する。動作パラメータ算出装置３０の動作パラメータ値受付部３０３は、指定された１つ又は２つの動作パラメータの値を受け付ける（Ｓ１０５）。
　システム行列決定部３０５は、指定された１つ又は２つの動作パラメータの値、及び指定された１つ又は２つの動作パラメータの値に対応する動作パラメータ以外の動作パラメータについては所定の値（初期値）に基づいて、システム行列Ｗを決定する（Ｓ１０７）。システム行列Ｗは、指定された又は所定の検出素子数、投影数を（Ｎ＋ｊ）、（Ｎ＋ｋ）としたとき（ｊ、ｋは０以上の整数）、（Ｎ＋ｋ）方向の前記検出方向それぞれに対応させて、波又は粒子のうち（Ｎ＋ｊ）個の検出素子２２のそれぞれに向かう部分の行路に対するＮ画素×Ｎ画素の二次元断層画像の各画素の寄与を示す重み値を要素にもつ（Ｎ＋ｋ）×（Ｎ＋ｊ）個の、Ｎ×Ｎ個の要素をもつ列又は行ベクトルを、二次元断層画像の画素に対応させて行又は列方向に並べたものである。
　正方システム行列生成部３０７は、ｊ＝ｋ＝０以外の場合に、決定されたシステム行列Ｗから、所定の計（ｊ＋ｋ）×Ｎ＋ｊ×ｋ個の行ベクトルの要素を除去することにより正方システム行列ＷSQを生成する。
　このために、まず、決定されたシステム行列Ｗを構成する行ベクトルの要素の中で、重複する行ベクトルの要素を除去する（Ｓ１０９）。
　次いで、残った行ベクトルの要素について、残す行ベクトルの要素か、間引く行ベクトルの要素かのいずれかを数値生成された乱数でランダムに決定し、計（ｊ＋ｋ）×Ｎ＋ｊ×ｋ個の行ベクトル又は列ベクトルの要素を除去することにより正方システム行列ＷSQを生成する（Ｓ１１１）。具体的には、例えば、残った行ベクトルの要素の行を特定するインデックス（行番号）に対して、残す行ベクトルの要素のインデックスか、間引く行ベクトルの要素のインデックスかのいずれかを数値生成された乱数でランダムに決定し、間引き数列｛ｍ｝を決定する。そして、決定された間引き数列｛ｍ｝に対応する行ベクトルの要素を除去することにより正方システム行列ＷSQを生成する。
　逆行列判定部３０９は、生成された正方システム行列ＷSQについて、逆行列が存在するか否かを判定する（Ｓ１１３）。具体的には、例えば、正方システム行列ＷSQの階数を算出し、算出された階数がＮ×Ｎか否かで、逆行列が存在するか否かを判定する。
　逆行列が存在すると判定された場合、動作パラメータ算出結果出力部３１１は、指定された１つ又は２つの動作パラメータの値に対応する動作パラメータ以外の動作パラメータの所定の値を表示、音声出力等の形態で出力する（Ｓ１１５）。
　逆行列が存在しないと判定された場合、ステップＳ１０７に戻り、指定された１つ又は２つの動作パラメータの値に対応する動作パラメータ以外の動作パラメータについて異なる所定の値に対して、ステップＳ１０７以降のステップを行う。
　逆行列判定部３０９は、必要に応じてステップＳ１０７に戻り、指定された１つ又は２つの動作パラメータの値に対応する動作パラメータ以外の動作パラメータについて異なる所定の値に対して、ステップＳ１０７以降のステップを繰り返してもよい。例えば、ステップＳ１１３で逆行列が存在すると判定された場合でも、繰り返しを行うことによって、逆行列を持つような、指定された１つ又は２つの動作パラメータの値に対応する動作パラメータ以外の動作パラメータの範囲を求めることができる。
＜離散ラドン逆変換行列の生成処理＞
　図１２Ａ及び図１２Ｂは、本発明の実施形態に係る離散ラドン逆変換行列の生成処理の一例のフローチャートである。
　ユーザが、上記の動作パラメータ算出処理によって求められた動作パラメータに基づいて、又は上記の動作パラメータ算出処理によって求められた動作パラメータとは無関係に検出素子数Ｎ＋ｊ、投影数Ｎ＋ｋ、（１８０°－１８０°／投影数）未満のスキャン角度範囲φの３つの動作パラメータを、離散ラドン逆変換行列生成装置４０の入力部を介して指定する（Ｓ２０１）。離散ラドン逆変換行列生成装置４０の動作パラメータ値受付部４０３は、指定された３つの動作パラメータの値を受け付ける（Ｓ２０３）。
　システム行列決定部４０５は、指定された３つの動作パラメータの値に基づいて、システム行列Ｗを決定する（Ｓ２０５）。
　正方システム行列生成部４０７は、ステップＳ１０９～Ｓ１１１と同様にして、ｊ＝ｋ＝０以外の場合に、決定されたシステム行列Ｗから、所定の計（ｊ＋ｋ）×Ｎ＋ｊ×ｋ個の行ベクトルの要素を除去することにより正方システム行列ＷSQを生成する。このために、まず、決定されたシステム行列Ｗを構成する行ベクトルの要素の中で、重複する行ベクトルの要素を除去する（Ｓ２０７）。次いで、残った行ベクトルの要素について、残す行ベクトルの要素か、間引く行ベクトルの要素かのいずれかを数値生成された乱数でランダムに決定し、計（ｊ＋ｋ）×Ｎ＋ｊ×ｋ個の行ベクトル又は列ベクトルの要素を除去することにより正方システム行列ＷSQを生成する（Ｓ２０９）。そして、正方システム行列ＷSQを構成する行ベクトルの要素のうちで、ステップＳ２０７において除去された行ベクトルの要素と同じものが存在すると判定した場合、その正方システム行列ＷSQを構成する行ベクトルの要素及び除去された同じ行ベクトルの要素に対応する行を特定するインデックス（行番号）を記憶部４１３に記憶する（Ｓ２１１）。
　逆行列判定部４０９は、ステップＳ１１３と同様に、生成された正方システム行列ＷSQについて、逆行列が存在するか否かを判定する（Ｓ２１３）。逆行列が存在すると判定された場合（Ｓ２１３－Ｙｅｓ）、離散ラドン逆変換行列生成部４１１は、正方システム行列ＷSQ逆行列である離散ラドン逆変換行列ＷSQ -1を生成する（Ｓ２１５）。逆行列が存在しないと判定された場合（Ｓ２１３－Ｎｏ）、ステップＳ２０９に戻って、別の乱数に基づいて正方システム行列ＷSQを生成する（Ｓ２１７）。
＜再構成画像生成処理＞
　図１３Ａ及び図１３Ｂは、本発明の実施形態に係る再構成画像生成処理の一例のフローチャートである。
　照射装置１０と検出装置２０との間に、対象物２を配置する（Ｓ３０１）。再構成画像生成装置５０の制御部５０１は、照射装置１０及び検出装置２０を制御して、照射装置１０からのビーム４を対象物２に照射しながら、重複しないＮ＋ｊ方向で検出装置２０によって検出を行う（Ｓ３０３）。
　再構成画像生成装置５０のサイノグラム生成部５０３は、検出装置２０の各検出素子２２から得られる信号について、必要に応じて機器の物理的特性について補正やスケールの換算などの処理と配列順序の整理などを行い、サイノグラムを生成する（Ｓ３０５）。このサイノグラムを生成する処理が、照射及び検出の処理であり、離散ラドン変換である。
　再構成画像生成装置５０の第１ベクトル生成部５０５は、生成されたサイノグラムから、システム行列Ｗの列方向の並びに合わせるような順序で要素を並べ替えた列ベクトルである第１ベクトルＸ′を生成する（Ｓ３０７）。第１ベクトルＸ′は、サイノグラムの値を持つので（Ｎ＋ｋ）×（Ｎ＋ｊ）要素を持つ。なお、照射及び検出Ｓ３０１からベクトル化Ｓ３０７までは、信号のアナログデジタル変換による量子化、デジタル信号の処理や、必要に応じて一時的にデータを保存して処理することもできる。また、サイノグラムは、それを生成しうるデータが得られれば十分であり、必ずしもその表示を行う必要はなく、サイノグラムといえる明示的なデータを経由する必要もない。サイノグラムをベクトル化したものに相当する第１ベクトルが得られれば、照射及び検出Ｓ３０１以降ベクトル化Ｓ３０７まで等価な任意の処理を行う任意の処理（例えば、シミュレーションによる処理）も本実施形態に含まれている。
　再構成画像生成装置５０のベクトルデシメーション部５０７は、正方システム行列ＷSQを生成した際に除去された行ベクトルのインデックス（行番号）に対応する、第１ベクトルＸ′の要素を間引き（ベクトルデシメーション処理）、第２ベクトルＸSQ′を生成する（Ｓ３０９）。このベクトルデシメーション処理を行うことにより、正方システム行列ＷSQとの間で式（４）の関係を満たすことができるＮ×Ｎ要素をもつ列ベクトルである第２ベクトルＸSQ′を得ることができる。
　第２ベクトルＸSQ′は正方システム行列ＷSQとの間で式（４）の関係を満たすため、正方システム行列ＷSQの逆行列である離散ラドン逆変換行列ＷSQ -1を用いて式（５）の関係が成り立つ。再構成画像生成装置５０の第３ベクトル生成部５０９は、ステップＳ２１３５で生成された離散ラドン逆変換行列ＷSQ -1と第２ベクトルＸSQ′を用いて、式（５）に基づいて、離散ラドン逆変換処理によって再構成像を与える第３ベクトルＸSQを生成する（Ｓ３１１）。
　再構成画像生成装置５０のデベクタライズ部５１１は、この第３ベクトルＸSQの並びをベクトル化の逆操作によってＮ画素×Ｎ画素の画像に戻すデベクタライズ処理により再構成画像を生成する（Ｓ３１３）。これで１スライス分の断層画像が再構成される。
　再構成画像生成装置５０の別測定値第２ベクトル生成部５１３は、ステップＳ２１１で記憶部４１１に記憶されたインデックス（行番号）に基づいて、ステップＳ２１１において正方システム行列ＷSQを構成する行ベクトルの要素のうちで、ステップＳ２０７において除去された行ベクトルの要素と同じものが存在すると判定した場合の正方システム行列ＷSQを構成する行ベクトルの要素に対応する、第２ベクトルＸSQ′の要素の一部又は全部を、除去された同じ行ベクトルの要素に対応する第２ベクトルＸSQ′の要素に置き換えて、別測定値第２ベクトルを生成する（Ｓ３１５）。
　第３ベクトル生成部５０９は、別測定値第２ベクトルを用いて第３ベクトルＸSQを生成する（Ｓ３１７）。そして、デベクタライズ部５１１は、生成された第３ベクトルＸSQにデベクタライズ処理を行い、再構成画像を生成し、記憶部５１７に記憶する（Ｓ３１９）。
　このステップＳ３１５～Ｓ３１９を、ステップＳ２１１において正方システム行列ＷSQを構成する行ベクトルの要素のうちで、ステップＳ２０７において除去された行ベクトルの要素と同じものが存在すると判定した場合の正方システム行列ＷSQを構成する行ベクトルと、除去された同じ行ベクトルの要素との異なる組み合わせに対して所望の回数繰り返す（Ｓ３２１）。再構成画像生成装置５０の合成再構成画像生成部５１５は、生成された複数の再構成画像を重ね合わせて、画質が向上した合成再構成画像を生成する（Ｓ３２３）。
（計算機シミュレーションによる検証（実施例））
　本実施形態による画像再構成手法の有効性を計算機シミュレーションにより確認した。逆ラドン変換を含む処理全般は、対比のための従来の手法のもの（比較例）も含めて、数式処理ソフトウエアMathematica（Wolfram Research, Inc., Champaign, Illinois）及び画像処理ソフトウエアＩｍａｇｅＪ（アメリカ国立衛生研究所）を用いて実現した。照射、検出、サイノグラム生成のシミュレーション処理も含め、各実施例及び比較例のシミュレーション処理を、Ｉｎｔｅｌ社製Ｘｅｏｎ　Ｅ５プロセッサー（３．５ＧＨｚ、６―Ｃｏｒｅタイプ）のＣＰＵと６４ＧＢの主記憶装置とを搭載したコンピュータ、又はＡｐｐｌｅ社製Ｍ１プロセッサー（８－Ｃｏｒｅタイプ）のＣＰＵと１６ＧＢの主記憶装置とを搭載したコンピュータのいずれかを用いて行った。図１４は検証のために用いた検証用データ（８ｂｉｔ数値ファントム）を示す図（原画像という）である。また、図１５Ａ、図１６Ａ、図１７Ａは、上述の検出素子数を解像度よりも大きくした場合、投影数を解像度よりも大きくした場合、検出素子数及び投影数を解像度よりも大きくした場合、それぞれに対応するサイノグラムの例を示す図である。各図には画素の座標を表す画素番号の最初及び最後を示す数値を上辺及び左辺に示した。原画像（数値ファントム）は、計算機シミュレーションにおいて対象物のある断面における減衰率構造を模擬するために人為的に与える数値データであり、画素が空間位置に対応した解像度がＮの画像、本実施例では、Ｎ＝６４つまり６４×６４の画像として表現できる。図１４では吸収の強い領域を明るく、弱い領域を暗く表示している。計算機シミュレーションでは、断層撮像装置１０による照射及び検出の動作も模擬し、計算によって測定値に相当する数値データを得た。これを表すサイノグラムでは、検出強度の弱い値を明るく、強い値を暗く表示している。図１５Ａ、図１６Ａ、図１７Ａのサイノグラムは、各測定動作にあたる離散ラドン変換に基づき得られたものである。具体的には、以下の３パターンで確認した。
（１）検出素子数を解像度よりも大きくした場合の例として、Ｎ＝６４、追加検出素子数ｋ＝１２８として、（Ｎ＋ｋ）＝１９２個の検出素子２２を持つ検出装置２０によって、追加投影数ｊ＝０を採用して、スキャン角度範囲９８．４°において、（Ｎ＋ｊ）＝６４方向からの投影分にて離散ラドン変換を実行した。これを反映するため、図１５Ａのサイノグラムは１９２×６４素のサイズをもつ画像として表示している。
（２）投影数を解像度よりも大きくした場合の例として、Ｎ＝６４、追加検出素子数ｋ＝０として、（Ｎ＋ｋ）＝６４個の検出素子２２を持つ検出装置２０によって、追加投影数ｊ＝１を採用して、スキャン角度範囲１１８．１５°において、（Ｎ＋ｊ）＝６５方向からの投影分にて離散ラドン変換を実行した。これを反映するため、図１６Ａのサイノグラムは６４×６５画素のサイズをもつ画像として表示している。
（３）検出素子数及び投影数を解像度よりも大きくした場合の例として、Ｎ＝６４、追加検出素子数ｋ＝６４として、（Ｎ＋ｋ）＝１２８個の検出素子２２を持つ検出装置２０によって、追加投影数ｊ＝３を採用して、スキャン角度範囲９３．３°において、（Ｎ＋ｊ）＝６７方向からの投影分にて離散ラドン変換を実行した。これを反映するため、図１７Ａのサイノグラムは１２８×６７画素のサイズをもつ画像として表示している。
なお、本実施例では、検出素子の追加は、検出素子群の領域が変わらないように検出素子のサイズを小さくすることにより行っている。図１５Ａ、図１６Ａ、図１７Ａのサイノグラムは本実施形態の動作のためのものであるが、従来のＦＢＰ法やＭＬ－ＥＭ法の処理のためにも利用している。従来の手法では検出素子の追加や投影数の追加をすることは特段必要ではなくＮ×Ｎ画素のサイズで取得したサイノグラム（図示しない）を採用して再構成画像を得ることができる。しかし、このＮ×Ｎ画素のサイズのサイノグラムは本実施形態によって得られるサイノグラムと比較して、情報量が小さくなる。このため、ここでは、同じ情報量を持つサイノグラムから再構成した画像による、手法間の比較とするために、本実施形態による（Ｎ＋ｊ）×（Ｎ＋ｋ）画素のサイズのサイノグラムを採用している。ただし、従来のＦＢＰ法やＭＬ－ＥＭ法の処理においては、再構成される画像のサイズを、実施例により得られた画像のサイズと同一にするために、検出素子方向の画素数が解像度と同一になるように画素値の平均をとることにより画像のサイズを縮小したサイノグラムを用いた。
　本実施形態の検証として、上記の３つのパターンについて、上記の検出素子数（Ｎ＋ｊ）、投影数（Ｎ＋ｋ）、スキャン角度範囲φに対して、上記のステップＳ２０１～Ｓ２０９、Ｓ２１３～Ｓ２１５及びＳ３０１～Ｓ３１３の処理（ステップＳ３０１～Ｓ３０５については上述のように計算機シミュレーションによる模擬）を行い、再構成画像を取得した。行列デシメーション処理の間引き数列｛ｍ｝を決定するために採用した処理は、上述の乱択アルゴリズムである。当該アルゴリズムに基づき間引き数列｛ｍ｝を決定する処理を実際に多数回行ったところ、ステップＳ２１３において、階数がＮ×Ｎ（６４×６４）に一致せず（達せず）、ステップＳ２０９に戻るリトライ処理が必要となることがあった。実際に必要となったリトライ処理の割合は、（１）検出素子数を解像度よりも大きくした場合において１０回のアルゴリズム動作に対する平均で１～２回程度、（２）投影数を解像度よりも大きくした場合において１０回のアルゴリズム動作に対する平均で３５回程度、（３）検出素子数及び投影数を解像度よりも大きくした場合において１０回のアルゴリズム動作に対する平均で３２８７回程度、であった。実施例の確認に用いた計算機で階数の条件を満たす正方行列を生成するのにそれぞれ（１）検出素子数を解像度よりも大きくした場合において１０～２０分程度、（２）投影数を解像度よりも大きくした場合において２０～３０分程度、３）検出素子数及び投影数を解像度よりも大きくした場合において８～２４時間程度かかった。なお、投影数の追加が１程度と小さいく、かつ計測角度範囲が小さい場合には、全体の行数は大きく増加しないが、その中から行の組を選び出す際に、１つでも従属な行が含まれると階数の条件を満たさなくなるため、追加投影数、追加検出素子数が多い場合の行列に比べて、階数の条件を満たす行を選び出す際のリトライの回数が増える場合があった。このことは、計算処理時間が行数の大小のみで決まらず、各行の従属性によっても変化することを示している。特にスキャン角度範囲が狭い場合に、各行の独立性が低下すると考えられることから、このような処理時間の増減は理論的に想定される振る舞いとなっている。
　図１５Ｂ、図１６Ｂ、図１７Ｂは、上記のパターン（１）、（２）、（３）のそれぞれの実施例において得られた再構成画像である。図１５Ｃ、図１６Ｃ、図１７Ｃは、図１５Ａ、図１６Ａ、図１７Ａのサイノグラムに基づいて、従来のＦＢＰ法においてＳｈｅｐｐ－Ｌｏｇａｎフィルタを用いて得られた再構成画像である。図１５Ｄ、図１６Ｄ、図１７Ｄは、図１５Ａ、図１６Ａ、図１７Ａのサイノグラムに基づいて、従来の代数再構成法として逐次近似計算（ＭＬ－ＥＭ法）を１０回行って局所解が０とならない回数までとして得られた再構成画像である。また、図１５Ｅ、図１６Ｅ、図１７Ｅは、図１５Ａ、図１６Ａ、図１７Ａのサイノグラムに基づいて、基準としてＦＢＰ法においてフィルタを用いず得られた再構成画像である。
　図１４と図１５Ｂ、図１４と図１６Ｂ、図１４と図１７Ｂを各々比較すれば明らかなように、本実施形態においては再現性が極めて高い再構成画像が得られた。特に、各従来法による再構成画像にはスキャン角度範囲を制限したことによる影響が大きく現れているが、本実施形態による再構成画像にはこれが見られない。このような高い再現性は、リトライ処理の有無にかかわらず乱択アルゴリズムにより得られたどの間引き数列｛ｍ｝を採用した場合でも同様であった。これは、従来のＳｈｅｐｐ－Ｌｏｇａｎフィルタを用いたＦＢＰ法（図１５Ｃ、図１６Ｃ、図１７Ｃ）、及びＭＬ－ＥＭ法（図１５Ｄ、図１６Ｄ、図１７Ｄ）とは異なり、離散ラドン逆変換行列ＷSQ -1が代数的な厳密解となることに裏付けられている本実施形態の利点であり、それが画像自体の高い再現性を通じ確認された。
　次に、再構成像の再現性の比較を数値的指標に基づいて実施した。具体的には、各再構成画像の再現性について式（２１）のＩＳＮＲ（Improvement in Signal-to-Noise ratio）を用いて評価した。
Figure JPOXMLDOC01-appb-M000019
 ここで、
Figure JPOXMLDOC01-appb-M000020
 は、順に、図１４に示す原画像、対象とする図１５Ｂ～図１５Ｄ、図１６Ｂ～図１６Ｄ、図１７Ｂ～図１７Ｄのいずれかの再構成画像、及び単純な逆ラドン変換（フィルタ無しの逆投影）の再構成画像について、それぞれをベクトル表示したものである。単純な逆ラドン変換の再構成画像は、比較のための基準であり、その画像は、上述のように図１５Ｅ、図１６Ｅ、図１７Ｅに示されている。式（２１）を用いて単純な逆ラドン変換の再構成画像を基準に求めた図１５Ｂ～図１５Ｄ、図１６Ｂ～図１６Ｄ、図１７Ｂ～図１７Ｄそれぞれにおける再構成画像のＩＳＮＲの値を表３、表４、表５に示す。
Figure JPOXMLDOC01-appb-T000021
Figure JPOXMLDOC01-appb-T000022
Figure JPOXMLDOC01-appb-T000023
 ここで、本実施例のＩＳＮＲ値は、８ｂｉｔのデータに対して厳密に解けたことから、分母が０となった為に無限大に発散した。各ＩＳＮＲの値から、ＦＢＰ法や逐次近似再構成法を用いる従来例に比べ、本実施形態による再構成画像の再現性が大幅に高まっていることが数値的にも確認された。
　上記において、慣用的な表現での「１８０°のスキャン角度範囲」よりも小さいスキャン角度範囲を、本明細書及び特許請求の範囲においては、「（１８０°－１８０°／投影数）未満のスキャン角度範囲」と表現しているとした。上記実施形態においては、指定されるスキャン角度範囲を（１８０°－１８０°／投影数）未満として説明したが、指定されるスキャン角度範囲が１８０°未満の任意の適切な所定の値以下（例えば１７０°以下、１６５°以下等）に制限されたものとしてもよい。
　また、上記実施形態においては、スキャン角度範囲を（１８０°－１８０°／投影数）未満として説明したが、追加検出素子数ｊが１以上及び／又は追加投影数ｋが１以上の場合は、スキャン角度範囲が１８０°であっても、より小さな計算資源で、離散ラドン逆変換の代数的厳密解に裏付けられた再構成画像を生成することができる。
　上記実施形態においては、図１Ｂを参照して説明した検出素子２２に対するビーム４についてのその画素のもつ寄与を寄与の有無として１か０で決定する手法を採用した。しかしながら、図１Ｂでハッチを付した部分の面積を寄与とする手法に加え、他の実用に付されている手法は本実施形態においても採用できる。例えば、面積に代え、掃く長さに応じた重みを持たせたり、その掃く長さを効率よく評価するシドンのアルゴリズム（Siddon's algorithm）を採用して上記寄与を決定することも有利である（非特許文献３）。さらに、さらなる高画質化を行うために、アンチエイリアス処理できるSunnegardh and Danielsonの方法を利用することも有利である（非特許文献４）。
　また、上記実施形態においては、一例の断層撮像装置１００の動作により説明したが、本実施形態は、必ずしも照射装置１０を採用するものには限定されない。ＳＰＥＣＴ（単一光子放射断層撮影）など、照射装置を用いることなく対象物の断層撮影を行う他の手法においても、本実施形態を同様に適用することができる。
　以上、本発明について、例示のためにいくつかの実施形態に関して説明してきたが、本発明はこれに限定されるものでなく、本発明の範囲及び精神から逸脱することなく、形態及び詳細について、様々な変形及び修正を行うことができることは、当業者に明らかであろう。
　以上の説明に関連して、更に以下の各項を開示する。
（１）
　（１８０°－１８０°／投影数）未満のスキャン角度範囲で、離散逆ラドン変換によりＮ画素×Ｎ画素の断層画像を再構成する断層撮像システムの動作パラメータを求めるためのコンピュータにより実行される方法であって、
　前記断層撮像システムは、少なくとも一列に並んだ複数の検出素子をもつ検出装置の検出範囲に配置された対象物に対して、照射装置により前記検出装置に向けて前記検出装置の検出対象となる波または粒子を照射しながら、または照射装置によらずに生じた前記波または粒子を対象物の各部に透過させながら、前記対象物の各部を透過した後の前記波または粒子を前記検出素子それぞれにより受けて前記検出素子別の強度値を得る検出動作を、スキャン角度範囲にわたって、前記対象物からみた前記波または粒子の相対的な検出方向における互いに重複しない、投影数の方向にわたって前記検出方向を変更して繰り返して得られた前記検出素子別の強度値を取得するものであり、
　検出素子数、投影数、及び（１８０°－１８０°／投影数）未満のスキャン角度範囲のうちの指定された１つ又は２つの動作パラメータの値を受け付けるステップと、
　前記指定された１つ又は２つの動作パラメータの値、及び前記指定された１つ又は２つの動作パラメータの値に対応する動作パラメータ以外の動作パラメータについては所定の値に基づいて、システム行列を決定し、ここで、前記システム行列は、指定された又は所定の検出素子数、投影数を（Ｎ＋ｊ）、（Ｎ＋ｋ）としたとき（ｊは０以上の整数、ｋは０以上の整数）、（Ｎ＋ｋ）方向の前記検出方向それぞれに対応させて、前記波または粒子のうち（Ｎ＋ｊ）個の前記検出素子のそれぞれに向かう部分の行路に対するＮ画素×Ｎ画素の二次元断層画像の各画素の寄与を示す重み値を要素にもつ（Ｎ＋ｋ）×（Ｎ＋ｊ）個の、Ｎ×Ｎ個の要素をもつ列または行ベクトルを、前記二次元断層画像の画素に対応させて行または列方向に並べたものであるステップと、
　ｊ＝ｋ＝０以外の場合に、決定された前記システム行列から、所定の計（ｊ＋ｋ）×Ｎ＋ｊ×ｋ個の行ベクトルまたは列ベクトルの要素を除去することにより正方システム行列を生成するステップと、
　決定された前記正方システム行列について、逆行列が存在するか否かを判定するステップと、
　前記システム行列を決定するステップ、前記正方システム行列を生成するステップ、及び前記逆行列が存在するか否かを判定するステップを、必要に応じて、前記指定された１つ又は２つの動作パラメータの値に対応する動作パラメータ以外の動作パラメータについて異なる所定の値に対して繰り返して、前記正方システム行列の逆行列が存在する、前記指定された１つ又は２つ以外の動作パラメータに対応する動作パラメータの値を求めるステップと、
を含む方法。
（２）
　前記正方システム行列を生成するステップは、決定された前記システム行列から、まず重複する行ベクトルまたは列ベクトルの要素を除去し、所定の計（ｊ＋ｋ）×Ｎ＋ｊ×ｋ個の行ベクトルまたは列ベクトルの要素を除去することにより正方システム行列を生成するステップである項（１）に記載の方法。
（３）
　前記正方システム行列を生成するステップは、決定された前記システム行列から、まず重複する行ベクトルまたは列ベクトルの要素を除去し、次いで、残った行ベクトルまたは列ベクトルの要素についてランダムに行ベクトルまたは列ベクトルを除去して、所定の計（ｊ＋ｋ）×Ｎ＋ｊ×ｋ個の行ベクトルまたは列ベクトルの要素を除去することにより正方システム行列を生成するステップである項（２）に記載の方法。
（４）
　項（１）～（３）に記載の方法をコンピュータに実行させるためのプログラム。
（５）
　項（４）に記載のプログラムを記録したコンピュータ読み取り可能な記録媒体。
（６）
　（１８０°－１８０°／投影数）未満のスキャン角度範囲で、離散逆ラドン変換によりＮ画素×Ｎ画素の断層画像を再構成する断層撮像システムの動作パラメータを求めるための装置であって、
　前記断層撮像システムは、少なくとも一列に並んだ複数の検出素子をもつ検出装置の検出範囲に配置された対象物に対して、照射装置により前記検出装置に向けて前記検出装置の検出対象となる波または粒子を照射しながら、または照射装置によらずに生じた前記波または粒子を対象物の各部に透過させながら、前記対象物の各部を透過した後の前記波または粒子を前記検出素子それぞれにより受けて前記検出素子別の強度値を得る検出動作を、スキャン角度範囲にわたって、前記対象物からみた前記波または粒子の相対的な検出方向における互いに重複しない、投影数の方向にわたって前記検出方向を変更して繰り返して得られた前記検出素子別の強度値を取得するものであり、
　検出素子数、投影数、及び（１８０°－１８０°／投影数）未満のスキャン角度範囲のうちの指定された１つ又は２つの動作パラメータの値を受け付ける動作パラメータ値受付部と、
　前記指定された１つ又は２つの動作パラメータの値、及び前記指定された１つ又は２つの動作パラメータの値に対応する動作パラメータ以外の動作パラメータについては所定の値に基づいて、システム行列を決定し、ここで、前記システム行列は、指定された又は所定の検出素子数、投影数を（Ｎ＋ｊ）、（Ｎ＋ｋ）としたとき（ｊは０以上の整数、ｋは０以上の整数）、（Ｎ＋ｋ）方向の前記検出方向それぞれに対応させて、前記波または粒子のうち（Ｎ＋ｊ）個の前記検出素子のそれぞれに向かう部分の行路に対するＮ画素×Ｎ画素の二次元断層画像の各画素の寄与を示す重み値を要素にもつ（Ｎ＋ｋ）×（Ｎ＋ｊ）個の、Ｎ×Ｎ個の要素をもつ列または行ベクトルを、前記二次元断層画像の画素に対応させて行または列方向に並べたものであるシステム行列決定部と、
　ｊ＝ｋ＝０以外の場合に、決定された前記システム行列から、所定の計（ｊ＋ｋ）×Ｎ＋ｊ×ｋ個の行ベクトルまたは列ベクトルの要素を除去することにより正方システム行列を生成する正方システム行列生成部と、
　決定された前記正方システム行列について、逆行列が存在するか否かを判定する逆行列判定部と、
を備え、
　前記システム行列を決定すること、前記正方システム行列を生成すること、及び前記逆行列が存在するか否かを判定することを、必要に応じて、前記指定された１つ又は２つの動作パラメータの値に対応する動作パラメータ以外の動作パラメータについて異なる所定の値に対して繰り返して、前記正方システム行列の逆行列が存在する、前記指定された１つ又は２つ以外の動作パラメータに対応する動作パラメータの値を求める装置。
（７）
　前記正方システム行列生成部は、決定された前記システム行列から、まず重複する行ベクトルまたは列ベクトルの要素を除去し、所定の計（ｊ＋ｋ）×Ｎ＋ｊ×ｋ個の行ベクトルまたは列ベクトルの要素を除去することにより正方システム行列を生成する項（６）に記載の装置。
（８）
　前記正方システム行列生成部は、決定された前記システム行列から、まず重複する行ベクトルまたは列ベクトルの要素を除去し、次いで、残った行ベクトルまたは列ベクトルの要素についてランダムに行ベクトルまたは列ベクトルを除去して、所定の計（ｊ＋ｋ）×Ｎ＋ｊ×ｋ個の行ベクトルまたは列ベクトルの要素を除去することにより正方システム行列を生成する項（６）に記載の装置。
　本開示は、波又は粒子が例えば対象物に照射されたりすることによって、対象物自体を波又は粒子が透過するような任意の断層撮像手法のために利用可能であり、例えば、Ｘ線ＣＴスキャナー、ＳＰＥＣＴ（Single Photon Emission Computed Tomography）、光トモグラフィー(Optical Tomography)、光学ＣＴ装置、電子線トモグラフィー（ＴＥＭＴ）が挙げられる。
２　対象物
４　ビーム
１００　断層撮像システム
１０　照射装置
２０　検出装置
２２　検出素子
３０　動作パラメータ算出装置
３０１　動作パラメータ解釈部
３０３　動作パラメータ受付部
３０５　システム行列決定部
３０７　正方システム行列生成部
３０９　逆行列判定部
３１１　逆行列判定部
３０ａ　ＣＰＵ
３０ｂ　ＲＡＭ
３０ｃ　ＲＯＭ
３０ｄ　外部メモリ
３０ｅ　入力部
３０ｆ　出力部
３０ｇ　通信部
３０ｈ　システムバス
４０　離散ラドン逆変換行列生成装置
４０３　動作パラメータ受付部
４０５　システム行列決定部
４０７　正方システム行列生成部
４０９　逆行列判定部
４１１　離散ラドン逆変換行列生成部
４１３　記憶部
５０　再構成画像生成装置
５０１　制御部
５０３　サイノグラム生成部
５０５　第１ベクトル生成部
５０７　ベクトルデシメーション部
５０９　第３ベクトル生成部
５１１　デベクタライズ部
５１３　別測定値第２ベクトル生成部
５１５　合成再構成画像生成部
５１７　記憶部

Claims (10)

  1. 　（１８０°－１８０°／投影数）未満のスキャン角度範囲で、離散逆ラドン変換により断層画像を再構成する断層撮像システムの動作パラメータを求めるためのコンピュータにより実行される方法であって、
    　検出素子数、投影数、及び（１８０°－１８０°／投影数）未満のスキャン角度範囲のうちの指定された１つ又は２つの動作パラメータの値に基づいて、システム行列を決定するステップと、
    　前記システム行列から正方システム行列を生成するステップと、
    　前記正方システム行列について、逆行列が存在するか否かを判定するステップと、
    　前記各ステップを、必要に応じて繰り返して、前記正方システム行列の逆行列が存在する、前記指定された１つ又は２つ以外の動作パラメータに対応する動作パラメータの値を求めるステップとを含む方法。
  2. 　前記正方システム行列を生成するステップは、決定された前記システム行列から、まず重複する行ベクトル又は列ベクトルの要素を除去し、正方システム行列を生成するステップである請求項１に記載の方法。
  3. 　前記正方システム行列を生成するステップは、決定された前記システム行列から、まず重複する行ベクトル又は列ベクトルの要素を除去し、次いで、残った行ベクトル又は列ベクトルの要素についてランダムに行ベクトル又は列ベクトルを除去して、正方システム行列を生成するステップである請求項２に記載の方法。
  4. 　請求項１～３に記載の方法をコンピュータに実行させるためのプログラム。
  5. 　請求項４に記載のプログラムを記録したコンピュータ読み取り可能な記録媒体。
  6. 　（１８０°－１８０°／投影数）未満のスキャン角度範囲で、離散逆ラドン変換により断層画像を再構成する断層撮像システムの動作パラメータを求めるための装置であって、
    　検出素子数、投影数、及び（１８０°－１８０°／投影数）未満のスキャン角度範囲のうちの指定された１つ又は２つの動作パラメータの値に基づいてシステム行列を決定するシステム行列決定部と、
    　前記システム行列から正方システム行列を生成する正方システム行列生成部と、
    　前記正方システム行列について、逆行列の有無を判定する逆行列判定部と、
    を備え、前記システム行列を決定し、前記正方システム行列を生成し、前記逆行列の有無を判定することを、必要に応じて繰り返して、前記正方システム行列の逆行列が存在する、前記指定された１つ又は２つ以外の動作パラメータに対応する動作パラメータの値を求める装置。
  7. 　前記正方システム行列生成部は、決定された前記システム行列から、まず重複する行ベクトル又は列ベクトルの要素を除去し、正方システム行列を生成する請求項６に記載の装置。
  8. 　前記正方システム行列生成部は、決定された前記システム行列から、まず重複する行ベクトル又は列ベクトルの要素を除去し、次いで、残った行ベクトル又は列ベクトルの要素についてランダムに行ベクトル又は列ベクトルを除去して、正方システム行列を生成する請求項６に記載の装置。
  9. 　１８０°未満の所定の値のスキャン角度範囲で、離散逆ラドン変換により断層画像を再構成する断層撮像システムの動作パラメータを求めるためのコンピュータにより実行される方法であって、
    　検出素子数、投影数、及び１８０°未満の所定の値のスキャン角度範囲のうちの指定された１つ又は２つの動作パラメータの値に基づいて、システム行列を決定するステップと、
    　前記システム行列から正方システム行列を生成するステップと、
    　前記正方システム行列について、逆行列が存在するか否かを判定するステップと、
    　前記各ステップを、必要に応じて繰り返して、前記正方システム行列の逆行列が存在する、前記指定された１つ又は２つ以外の動作パラメータに対応する動作パラメータの値を求めるステップとを含む方法。
  10. 　１８０°未満の所定の値のスキャン角度範囲で、離散逆ラドン変換により断層画像を再構成する断層撮像システムの動作パラメータを求めるための装置であって、
    　検出素子数、投影数、及び１８０°未満の所定の値のスキャン角度範囲のうちの指定された１つ又は２つの動作パラメータの値に基づいてシステム行列を決定するシステム行列決定部と、
    　前記システム行列から正方システム行列を生成する正方システム行列生成部と、
    　前記正方システム行列について、逆行列の有無を判定する逆行列判定部と、
    を備え、前記システム行列を決定し、前記正方システム行列を生成し、前記逆行列の有無を判定することを、必要に応じて繰り返して、前記正方システム行列の逆行列が存在する、前記指定された１つ又は２つ以外の動作パラメータに対応する動作パラメータの値を求める装置。
PCT/JP2022/002247 2021-01-21 2022-01-21 断層撮像のためのシステム、方法、プログラム、及びプログラムを記録した記録媒体 WO2022158575A1 (ja)

Priority Applications (3)

Application Number Priority Date Filing Date Title
JP2022576763A JPWO2022158575A1 (ja) 2021-01-21 2022-01-21
EP22742694.7A EP4283343A1 (en) 2021-01-21 2022-01-21 System, method, and program for tomographic imaging, and recording medium in which program is recorded
US18/262,230 US20240153162A1 (en) 2021-01-21 2022-01-21 System, Method, and Program for Tomographic Imaging, and Recording Medium in which Program is Recorded

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2021-008321 2021-01-21
JP2021008321 2021-01-21

Publications (1)

Publication Number Publication Date
WO2022158575A1 true WO2022158575A1 (ja) 2022-07-28

Family

ID=82548766

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2022/002247 WO2022158575A1 (ja) 2021-01-21 2022-01-21 断層撮像のためのシステム、方法、プログラム、及びプログラムを記録した記録媒体

Country Status (4)

Country Link
US (1) US20240153162A1 (ja)
EP (1) EP4283343A1 (ja)
JP (1) JPWO2022158575A1 (ja)
WO (1) WO2022158575A1 (ja)

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20090123048A1 (en) * 2007-05-09 2009-05-14 Jean-Daniel Leroux Image Reconstruction Methods Based on Block Circulant System Matrices
WO2013105583A1 (ja) 2012-01-10 2013-07-18 株式会社 東芝 逐次近似法を用いたｘ線コンピュータ断層撮影装置（ｘ線ｃｔ装置）
WO2014021349A1 (ja) 2012-07-30 2014-02-06 株式会社 東芝 Ｘ線コンピュータ断層撮影装置及び画像再構成方法
CN104899827A (zh) * 2015-05-26 2015-09-09 大连理工大学 基于固定分辨率条件下的离散Radon投影和Mojette投影转换方法
JP2015226753A (ja) * 2014-06-02 2015-12-17 株式会社東芝 Ｘ線ｃｔ装置及びその制御方法
WO2019230740A1 (ja) * 2018-05-28 2019-12-05 国立研究開発法人理化学研究所 オーバーサンプリングによる断層画像データの取得方法、取得装置、および制御プログラム
WO2019230741A1 (ja) * 2018-05-28 2019-12-05 国立研究開発法人理化学研究所 角度オフセットによる断層画像データの取得方法、取得装置、および制御プログラム

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20090123048A1 (en) * 2007-05-09 2009-05-14 Jean-Daniel Leroux Image Reconstruction Methods Based on Block Circulant System Matrices
WO2013105583A1 (ja) 2012-01-10 2013-07-18 株式会社 東芝 逐次近似法を用いたｘ線コンピュータ断層撮影装置（ｘ線ｃｔ装置）
WO2014021349A1 (ja) 2012-07-30 2014-02-06 株式会社 東芝 Ｘ線コンピュータ断層撮影装置及び画像再構成方法
JP2015226753A (ja) * 2014-06-02 2015-12-17 株式会社東芝 Ｘ線ｃｔ装置及びその制御方法
CN104899827A (zh) * 2015-05-26 2015-09-09 大连理工大学 基于固定分辨率条件下的离散Radon投影和Mojette投影转换方法
WO2019230740A1 (ja) * 2018-05-28 2019-12-05 国立研究開発法人理化学研究所 オーバーサンプリングによる断層画像データの取得方法、取得装置、および制御プログラム
WO2019230741A1 (ja) * 2018-05-28 2019-12-05 国立研究開発法人理化学研究所 角度オフセットによる断層画像データの取得方法、取得装置、および制御プログラム

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
AVIASH C. KAKMOLCOLM SLANEY: "Principles of Computerized Tomographic Imaging", 1988, IEEE PRESS
N. KAWASEM. KATOH. NISHIOKAH. JINNAI: "Transmission electron microtomography without the ''missing wedge'' for quantitative structural analysis", ULTRAMICROSCOPY, vol. 107, 2007, pages 8 - 15, XP005781314, DOI: 10.1016/j.ultramic.2006.04.007
ROBERT L. SIDDON: "Fast calculation of the exact radiological path for a three-dimensional CT array", MEDICAL PHYSICS, vol. 12, no. 2, 1985, pages 252 - 255, XP000577395, DOI: 10.1118/1.595715
SERHAN O. ISIKAMWAHEB BISHARASAM MAVANDADIFRANK W. YUSTEVE FENGRANDY LAUAYDOGAN OZCAN: "Lens-free optical tomographic microscope with a large imaging volume on a chip", PNAS, vol. 108, no. 18, 3 May 2011 (2011-05-03), pages 7296 - 7301, XP055121935, Retrieved from the Internet <URL:https://di.org/10.1073/pnas.1015638108> DOI: 10.1073/pnas.1015638108
STANLEY R. DEANS: "The Radon Transform and Some of its Applications", 1983, KRIEGER PUBLICATIONS
SUNNEGARDH, J.DANIELSSON, P.-E.: "A New Anti-Aliased Projection Operator for Iterative CT Reconstruction", 9TH INTERNATIONAL MEETING OF FULLY THREE-DIMENSIONAL IMAGE RECONSTRUCTION IN RADIOLOGY AND NUCLEAR MEDICINE, 2009, pages 125 - 127

Also Published As

Publication number Publication date
EP4283343A1 (en) 2023-11-29
US20240153162A1 (en) 2024-05-09
JPWO2022158575A1 (ja) 2022-07-28

Similar Documents

Publication Publication Date Title
RU2709437C1 (ru) Способ обработки изображений, устройство обработки изображений и носитель данных
US9672638B2 (en) Spectral X-ray computed tomography reconstruction using a vectorial total variation
US8135186B2 (en) Method and system for image reconstruction
US7840053B2 (en) System and methods for tomography image reconstruction
Jia et al. GPU-based fast low-dose cone beam CT reconstruction via total variation
EP2862516B1 (en) CT imaging methods and systems
US8571287B2 (en) System and method for iterative image reconstruction
JP7236110B2 (ja) オーバーサンプリングによる断層画像データの取得方法、取得装置、および制御プログラム
RU2510080C2 (ru) Устройство для обработки изображения, способ обработки изображения и среда долговременного хранения информации
JP2020036877A (ja) 反復的画像再構成フレームワーク
JP2020516345A (ja) 深層学習に基づくトモグラフィ再構成
WO2020237873A1 (zh) 基于神经网络的螺旋ct图像重建方法和设备及存储介质
KR20100133950A (ko) 동적인 제약들에 따른 오브젝트 주변의 사용을 통한 단층 촬영에 있어서의 양 감소 및 이미지 강화
Shi et al. Review of CT image reconstruction open source toolkits
JP7273272B2 (ja) 角度オフセットによる断層画像データの取得方法、取得装置、および制御プログラム
Wang et al. Sparse-view cone-beam CT reconstruction by bar-by-bar neural FDK algorithm
Hou et al. Analysis of compressed sensing based CT reconstruction with low radiation
WO2022158575A1 (ja) 断層撮像のためのシステム、方法、プログラム、及びプログラムを記録した記録媒体
Yan et al. Comparison of heuristic and deterministic algorithms in neutron coded imaging reconstruction
Mohan et al. 4D model-based iterative reconstruction from interlaced views
Buzmakov et al. Efficient and effective regularised ART for computed tomography
JP5493072B2 (ja) Ｃｔ装置、ｃｔ装置における画像再構成方法、及び電子回路部品
Iskender et al. Scatter correction in X-ray CT by physics-inspired deep learning
Saha et al. Multi-axial CT reconstruction from few view projections
Hou et al. A compressed sensing approach to low-radiation CT reconstruction

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: 22742694

Country of ref document: EP

Kind code of ref document: A1

ENP Entry into the national phase

Ref document number: 2022576763

Country of ref document: JP

Kind code of ref document: A

NENP Non-entry into the national phase

Ref country code: DE

ENP Entry into the national phase

Ref document number: 2022742694

Country of ref document: EP

Effective date: 20230821