WO2012050149A1 - Ｘ線ｃｔ装置及び画像再構成方法 - Google Patents

Ｘ線ｃｔ装置及び画像再構成方法 Download PDF

Info

Publication number
WO2012050149A1
WO2012050149A1 PCT/JP2011/073473 JP2011073473W WO2012050149A1 WO 2012050149 A1 WO2012050149 A1 WO 2012050149A1 JP 2011073473 W JP2011073473 W JP 2011073473W WO 2012050149 A1 WO2012050149 A1 WO 2012050149A1
Authority
WO
WIPO (PCT)
Prior art keywords
projection data
ray
matrix
successive approximation
extended
Prior art date
Legal status (The legal status 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 status listed.)
Ceased
Application number
PCT/JP2011/073473
Other languages
English (en)
French (fr)
Inventor
高橋　悠
後藤　大雅
廣川　浩一
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Hitachi Healthcare Manufacturing Ltd
Original Assignee
Hitachi Medical Corp
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 Hitachi Medical Corp filed Critical Hitachi Medical Corp
Priority to CN201180049466.5A priority Critical patent/CN103153192B/zh
Priority to US13/824,697 priority patent/US9123156B2/en
Priority to JP2012538707A priority patent/JP5828841B2/ja
Publication of WO2012050149A1 publication Critical patent/WO2012050149A1/ja
Anticipated expiration legal-status Critical
Ceased legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T12/00—Tomographic reconstruction from projections
    • G06T12/20—Inverse problem, i.e. transformations from projection space into object space
    • 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]
    • A61B6/032—Transmission 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/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/52—Devices using data or image processing specially adapted for radiation diagnosis
    • A61B6/5258—Devices using data or image processing specially adapted for radiation diagnosis involving detection or reduction of artifacts or noise
    • A61B6/5264—Devices using data or image processing specially adapted for radiation diagnosis involving detection or reduction of artifacts or noise due to motion
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2211/00—Image generation
    • G06T2211/40—Computed tomography
    • G06T2211/424—Iterative

Definitions

  • the present invention relates to an X-ray CT apparatus or the like for reconstructing an image by a successive approximation method.
  • X-ray CT equipment irradiates a subject with fan beam (fan beam) or cone beam (cone or pyramid beam) X-rays, measures the X-rays transmitted through the subject with an X-ray detector, This is an apparatus for obtaining a tomographic image of a subject by reconstructing measurement data from multiple directions.
  • the image reconstruction method in the X-ray CT apparatus is roughly divided into an analysis method and a successive approximation method.
  • the analysis method is a method of solving a problem analytically based on the projected cutting plane theorem.
  • the successive approximation method is a method of mathematically modeling an observation system that leads to acquisition of projection data and estimating the best image by an iterative method based on the mathematical model.
  • the advantage of the analysis method is that the reconstructed image can be obtained directly from the actual projection data, so that the amount of calculation is overwhelmingly small.
  • the advantage of the successive approximation method is that the physical process leading to the acquisition of projection data and the statistical fluctuations included in the actual projection data can be considered as mathematical models and statistical models, respectively. Beam artifacts) and quantum noise on the image can be reduced.
  • the Feldkamp method which is an analysis method, or a method obtained by improving the Feldkamp method is mainly used because of the small amount of calculation.
  • the practical use of the successive approximation method has begun to be studied.
  • the successive approximation method is a method in which an evaluation index of an image is set in advance and the image is sequentially updated so that an evaluation value obtained by quantifying the evaluation index takes a maximum value or a minimum value.
  • the evaluation index the contradiction between the forward projection data obtained by converting the image into projection data during the update process and the actual projection data, probabilistic likelihood, etc. are used, and the function for calculating the evaluation value is called the evaluation function .
  • Non-Patent Document 1 proposes a successive approximation method using a penalized weighted square error function as an evaluation function. As a method proposed so far, as proposed in Non-Patent Document 1, it is common to operate matrices that are transposed in the forward projection process and the backprojection process.
  • Non-Patent Document 2 proposes a method of performing sequential updating by applying view direction weights to back projection processing.
  • view direction weighted back projection process the back projection process using the view direction weight.
  • the view direction weighted backprojection process itself is a technique proposed in Non-Patent Document 3 and used in the analysis method.
  • the view direction weighted backprojection processing has the following advantages. (1) The redundancy of projection data can be eliminated. (2) Time resolution can be improved.
  • the update formula for performing the sequential update includes a relaxation coefficient related to the convergence speed and stability of the successive approximation method. In order for the successive approximation method to converge stably, it is necessary to set the relaxation coefficient within a specific range.
  • Non-patent document 2 describes that the relaxation coefficient is determined empirically.
  • Non-Patent Document 4 proposes a method of calculating a relaxation coefficient using a power method.
  • the power method is an iterative method for obtaining the maximum eigenvalue of an arbitrary matrix.
  • Patent Document 1 describes a method of performing back projection processing after generating actual column data and virtual channel data by expanding actual projection data in the analysis method. If the method of Patent Document 1 can be applied to the successive approximation method, the restriction on the region can be relaxed.
  • Non-Patent Document 2 As in the method of Non-Patent Document 2, by using the view direction weighted back projection processing by the successive approximation method, it is possible to reduce the deterioration of the image quality due to the movement of the subject. However, the method of Non-Patent Document 2 has a problem caused by the operator empirically determining the relaxation coefficient.
  • the sufficient condition of the relaxation coefficient for convergence of the successive approximation method changes depending on the bed feeding speed and the shooting FOV (Field Of View) as shooting conditions. It is very cumbersome and time consuming for the operator to determine the relaxation coefficient empirically according to these many imaging conditions.
  • the relaxation coefficient when the relaxation coefficient is set to a small value, the convergence condition of the successive approximation process is satisfied regardless of the shooting conditions due to the nature of the update matrix. However, if the relaxation coefficient is set to a too small value, the convergence of the successive approximation process is delayed, and the image quality of the reconstructed image is also deteriorated.
  • the relaxation coefficient is determined using the power method
  • the relaxation coefficient is automatically calculated according to the shooting conditions, so that the successive approximation method itself converges stably.
  • the method of Non-Patent Document 4 requires the sequential calculation of a large-scale matrix operation, and the calculation time Will increase.
  • Non-Patent Document 2 and Non-Patent Document 4 it is necessary to set an optimal value for the relaxation coefficient. That is, if the relaxation coefficient is not set to an optimum value, there remains a problem that the calculation time increases.
  • the present invention has been made in view of the above-described problems, and a first object thereof is to provide an X-ray CT apparatus and the like that reconstructs an image by a successive approximation method that can be stably converged and executed at high speed. That is.
  • the second object is to provide an X-ray CT apparatus and the like that can suppress image quality degradation when the successive approximation method is applied to data with body movement.
  • a third object is to provide an X-ray CT apparatus or the like that can suppress image quality degradation when the successive approximation method is applied under imaging conditions in which data loss occurs during spiral scanning or axial scanning with a high bed moving speed. It is.
  • the first invention is an X-ray CT apparatus that reconstructs a tomographic image of a subject by a successive approximation method in which matrices that are not transposed are applied to each other in forward projection processing and backprojection processing.
  • the actual projection data using an imaging unit that acquires actual projection data of the subject based on the imaging conditions, and an update formula of the successive approximation method that includes the imaging conditions and a relaxation coefficient that determines convergence of calculation.
  • an arithmetic unit that reconstructs the tomographic image by performing successive approximation on the X-ray CT apparatus, wherein the relaxation coefficient is analytically calculated based on the imaging conditions .
  • a second invention is an image reconstruction method for reconstructing a tomographic image of a subject by a successive approximation method in which matrices that are not transposed are applied to each other in forward projection processing and backprojection processing.
  • an X-ray CT apparatus that reconstructs an image by a successive approximation method that converges stably and can be executed at high speed.
  • an X-ray CT apparatus or the like that can suppress image quality degradation when the successive approximation method is applied to data with body movement.
  • an X-ray CT apparatus and the like that can suppress image quality deterioration when the successive approximation method is applied under imaging conditions in which data loss occurs during spiral scanning or axial scanning with a high bed moving speed.
  • the X-ray CT apparatus 1 includes a scanner 2 (imaging unit) on which an X-ray tube 11 and a detector 12 are mounted, a bed 4 on which a subject 10 is placed, and data obtained from the detector 12. And an input device 6 such as a mouse, a trackball, a keyboard, and a touch panel, a display device 7 that displays a reconstructed image, and the like.
  • a scanner 2 imaging unit
  • an input device 6 such as a mouse, a trackball, a keyboard, and a touch panel
  • a display device 7 that displays a reconstructed image, and the like.
  • the operator inputs shooting conditions and reconstruction parameters via the input device 6.
  • the imaging conditions are, for example, a bed feeding speed, a tube current, a tube voltage, and a slice position.
  • the reconstruction parameter is, for example, a region of interest, a reconstructed image size, a backprojection phase width, a reconstruction filter function, or the like.
  • the X-ray CT apparatus 1 is roughly composed of a scanner 2, an operation unit 3, and a bed 4.
  • the scanner 2 includes an X-ray tube 11 (X-ray generator), a detector 12, a collimator 13, a drive device 14, a central controller 15, an X-ray controller 16, a high voltage generator 17, a scanner controller 18, and a bed control.
  • the apparatus 19 includes a bed movement measuring apparatus 20, a collimator control apparatus 21, a preamplifier 22, an A / D converter 23, and the like.
  • the central controller 15 inputs imaging conditions and reconstruction parameters from the input device 6 in the operation unit 3, and sends control signals necessary for imaging to the collimator controller 21, the X-ray controller 16, the scanner controller 18, and the bed control. Transmit to device 19.
  • the collimator control device 21 controls the position of the collimator 13 based on the control signal.
  • the X-ray control device 16 controls the high voltage generator 17 based on the control signal.
  • the high voltage generator 17 applies a tube voltage and a tube current to the X-ray tube 11 (X-ray generator).
  • X-ray generator In the X-ray tube 11, electrons with energy corresponding to the applied tube voltage are emitted from the cathode, and the emitted electrons collide with the target (anode), whereby X-rays with energy corresponding to the electron energy are Is irradiated.
  • the scanner control device 18 controls the drive device 14 based on the control signal.
  • the driving device 14 circulates around the subject 10 around a gantry portion on which the X-ray tube 11, the detector 12, the preamplifier 22, and the like are mounted.
  • the bed control device 19 controls the bed 4 based on the control signal.
  • the X-ray irradiated from the X-ray tube 11 is limited in the irradiation region by the collimator 13, is absorbed (attenuated) according to the X-ray attenuation coefficient in each tissue in the subject 10, passes through the subject 10, It is detected by a detector 12 arranged at a position facing the tube 11.
  • the detector 12 includes a plurality of detection elements arranged in a two-dimensional direction (a channel direction and a column direction perpendicular to the channel direction). X-rays received by each detection element are converted into projection data (hereinafter referred to as “actual projection data”).
  • the X-rays detected by the detector 12 are converted into current, amplified by the preamplifier 22, converted into digital data by the A / D converter 23, LOG converted, calibrated, and used as actual projection data. Input to the arithmetic unit 5.
  • the computing device 5 includes a reconstruction computing device 31, an image processing device 32, and the like.
  • the input / output device 9 includes an input device 6, a display device 7, a storage device 8 (storage unit), and the like.
  • the reconstruction calculation device 31 performs an image reconstruction process using the actual projection data, and generates a reconstructed image.
  • the reconstruction calculation device 31 generates filter correction projection data by superimposing a reconstruction filter on the actual projection data of each view, and weights the filter correction projection data in the view direction (hereinafter referred to as “view direction weight”). .) Is weighted and the back projection process is performed to form a non-destructive tomographic image as a distribution map of the X-ray attenuation coefficient inside the subject 10.
  • the reconstruction calculation device 31 stores the generated reconstruction image in the storage device 8. In addition, the reconstruction calculation device 31 displays a reconstruction image as a CT image on the display device 7. Alternatively, the image processing device 32 performs image processing on the reconstructed image stored in the storage device 8, and displays the reconstructed image subjected to the image processing on the display device 7 as a CT image.
  • the X-ray CT apparatus 1 is a multi-slice CT that uses a detector 12 in which detector elements are arranged in a two-dimensional direction, and a single that uses a detector 12 in which detector elements are arranged in one row, that is, in a one-dimensional direction (channel direction only). Broadly divided into slice CT.
  • multi-slice CT an X-ray beam spreading in a cone shape or a pyramid shape is irradiated from an X-ray tube 11 as an X-ray source in accordance with the detector 12.
  • an X-ray beam spreading in a fan shape is emitted from the X-ray tube 11.
  • X-ray irradiation is performed while the gantry section circulates around the subject 10 placed on the bed 4 (however, scanogram imaging is excluded). ).
  • imaging in which the bed 4 is fixed during imaging and the X-ray tube 11 circulates around the subject 10 in a circular orbit is called an axial scan or the like.
  • imaging in which the bed 4 moves and the X-ray tube 11 circulates around the subject 10 in a spiral trajectory is called a spiral scan.
  • the bed controller 20 keeps the bed 4 stationary in the case of the axial scan. Further, in the case of the helical scan, the bed control device 20 translates the bed 4 in the body axis direction in accordance with the bed feeding speed as an imaging condition input via the input device 6.
  • the X-ray CT apparatus 1 is, for example, a multi-slice CT. Further, the scan method of the X-ray CT apparatus 1 is, for example, a rotate-rotate method (third generation).
  • Non-Patent Document 1 a successive approximation method using a penalized weighted square error function as an evaluation function.
  • x (k) is a vector (image vector) representing an image in the k-th sequential update
  • y is a vector representing actual projection data.
  • A is a matrix that associates an image with projection data, and is referred to as a system matrix because it represents the characteristics of the imaging system via the mathematical model described above.
  • a T represents a transposed matrix of A.
  • Ax (k) corresponds to processing (forward projection processing) for converting the image vector x (k) into a projection data vector.
  • a T ( ⁇ ) corresponds to processing (back projection processing) for converting the projection data vector in parentheses into an image vector.
  • D is a diagonal matrix having a weighting factor weighted to the difference value between the actual projection data and the forward projection data as a diagonal component.
  • a value corresponding to the number of detected photons is used as the weighting factor.
  • ⁇ is an optional parameter that adjusts the strength of the penalty term.
  • the penalty term has a role of suppressing the high-frequency enhancement effect of the image by applying the successive approximation method.
  • S is a matrix having the inverse of each element of the vector s represented by the following equation as a diagonal component.
  • c is a vector having the same number of elements as the image vector and having 1 as a value for all elements.
  • R in Equation (1) and R ′ in Equation (2) are linear operators of the first and second derivatives of the penalty term, and both are conversion matrices from image vectors to image vectors.
  • Non-Patent Document 2 for example.
  • the view direction weighted back projection process itself is proposed in Non-Patent Document 3.
  • Non-Patent Document 2 is expressed by the following equation.
  • ⁇ is a relaxation coefficient related to the convergence speed and stability of the successive approximation method.
  • equation (5) is modified as follows.
  • the matrix Q in the equation (7) is called an update matrix.
  • the update matrix is a matrix that determines an image update amount and an update direction per one-time sequential update process based on the update formula.
  • I represents a unit matrix.
  • the power eigenmethod is used to calculate the maximum eigenvalue of the update matrix, thereby calculating the spectrum radius ⁇ (Q) shown in Equation (8), and the relaxation coefficient. ⁇ is set.
  • the power method is not a method for solving the problem analytically, and it is necessary to sequentially calculate a large-scale matrix operation, which increases the calculation time.
  • FIG. 4 (a) shows the arrangement of the detector 12 and the X-ray tube 11 in the imaging section. Moreover, scanFOV41 in the imaging
  • the scanFOV41 is a region where an image is formed based on the projection cut surface theorem.
  • the scanFOV41 in the photographing section is almost circular.
  • FIG. 4 (b) shows the arrangement of the detector 12 and the X-ray tube 11 in the body axis direction. Moreover, scanFOV41 in the body axis direction is shown.
  • the scanFOV 41 in the body axis direction is a polygon such as a triangle, a quadrangle, or a pentagon, although it varies depending on the bed moving speed.
  • Fig. 5 shows the difference in scanFOV41 in the body axis direction due to the difference in bed movement speed. Comparing FIG. 5 (a) and FIG. 5 (b), it can be seen that the faster the bed moving speed, the narrower the width of the scanFOV 41 in the body axis direction. If the width of scanFOV41 in the body axis direction is narrower than the specified value, it means that the 180-degree view necessary to reconstruct a tomographic image has not been acquired, which limits the area where images can be created. Means that
  • the area in which an image can be created is also limited in the peripheral visual field in the body axis direction during axial scanning.
  • a relaxation coefficient value for stably converging the successive approximation method is analytically calculated. Then, an image reconstruction process is performed by a successive approximation method that performs a view direction weighted backprojection process.
  • the view direction weighted backprojection process uses the scan characteristics of the X-ray CT apparatus 1 shown below.
  • FIG. 6A shows the transmission path of the target view and the transmission path of the opposite phase view in the channel direction.
  • FIG. 6B shows the transmission path of the target view and the transmission path of the opposite phase view in the column direction.
  • data having such a relationship will be referred to as actual projection data of opposite phase.
  • the X-ray CT apparatus 1 collects data while rotating the X-ray tube 11 and the detector 12, an imaging time difference occurs in the view direction of the actual projection data. If the bed moving speed during imaging is constant, the imaging time difference is proportional to the distance between the body axis direction position of the X-ray tube 11 and the body axis direction position of the image.
  • FIG. 7 is a conceptual diagram of view direction weights.
  • the arithmetic device 5 of the X-ray CT apparatus 1 has, as the view direction weight, according to the distance between the body axis direction position of the slice including the pixel of interest and the body axis direction position of the X-ray tube 11, A weight is added so as to add 1 between the actual projection data in the opposite phase.
  • the arithmetic unit 5 of the X-ray CT apparatus 1 calculates the view direction weight by a predetermined number of views on both sides centering on the view in the body axis direction position of the X-ray tube 11 that is closest to the target pixel in the body axis direction. 0 weights except for the view.
  • a view to which a weight other than 0 is attached is referred to as a backprojection view
  • the number of backprojection views is referred to as a backprojection view number.
  • the optimal number of backprojection views is determined by the helical pitch (value determined by the bed feeding speed and the gantry rotation speed as imaging conditions) and the length of the detector 12 in the row direction.
  • the view direction weight By defining the view direction weight in this way, the time resolution can be improved as compared with the case where the view direction weight is not used. This is because the time phases of a plurality of views that contribute to the image reconstruction process are substantially matched by the view direction weight.
  • the definition of view direction weights may use a common backprojection view number for all pixels and determine a central view for each slice. Further, the definition of the view direction weight may determine the central view for each pixel as disclosed in “JP 2004-199163”. Further, the definition of the view direction weight may determine the number of backprojection views in each pixel.
  • the view direction weight is the actual projection data Since it is set for each element of the vector and the image vector, it is defined as w ij .
  • the elements of the system matrix A and a ij, the elements of the system matrix B multiplied by the view direction weight putting and b ij, b ij can be expressed as the following equation .
  • the back projection matrix representing the view direction weighted back projection processing is B T.
  • the system matrix A is a matrix that associates an image with projection data.
  • the computing device 5 of the X-ray CT apparatus 1 calculates each element a ij based on the apparatus specifications of the X-ray CT apparatus 1 and the imaging conditions input from the input device 6.
  • a calculation method of the system matrix A for example, the Joseph method proposed in “Joseph, PM (1982). An Improved Algorithm for Reprojecting Rays Through Pixel Images. IEEE Transactions on Medical Imaging MI-1, 192-196.” Is mentioned.
  • Examples of the apparatus specifications of the X-ray CT apparatus 1 used for the calculation of the system matrix A include the following. -Number of elements in the column direction of the detector 12-Number of elements in the channel direction of the detector 12-Size of one element in the column direction of the detector 12-Size of one element in the channel direction of the detector 12-Between the X-ray tube 11 and the rotation center Distance d sod The distance d sid between the X-ray tube 11 and the center of the detector 12 (detector center)
  • imaging conditions used for the calculation of the system matrix A include the following.
  • Beam pitch value obtained by dividing the bed feeding speed by the row direction detector size at the rotation center: “value corresponding to high speed” or “value corresponding to low speed”.
  • -Shooting rotation speed-Shooting start angle-Rotation speed-Number of shooting views per rotation-Shooting method “Axial scan” or "Helix scan”.
  • -Number of pixels (in the cross section): Number of pixels in the horizontal direction x Number of pixels in the vertical direction-Number of image slices-Pixel size (in the cross section) -Image slice thickness and image center position (within the cross section): “Matches the rotation center”. -Image slicing direction center position: “Same as the number of rotations of shooting 2.5.”
  • the system matrix A is different when the beam pitch, which is one of the imaging conditions, is “a value corresponding to high speed” and “a value corresponding to low speed”.
  • the relaxation coefficient ⁇ at which the successive approximation method converges stably is also different.
  • the spectrum radius of the update matrix is calculated by the power method, but there is a problem that the calculation time and the amount of memory used for the calculation are large. Therefore, the X-ray CT apparatus 1 according to the present invention calculates the operator norm of the update matrix instead of calculating the spectrum radius of the update matrix.
  • the matrix operator norm has the following properties.
  • equation (8) can be replaced as follows: Obviously, when the operator norm of the update matrix is smaller than 1 (when the equation (13) is satisfied), the spectral radius of the update matrix is smaller than 1 (the equation (8) is satisfied).
  • Non-Patent Document 3 an example of introducing a view direction weighted back projection process and a relaxation coefficient to the method proposed in Non-Patent Document 3 will be given.
  • the present invention is not limited to this example, and can be applied to any successive approximation method including an update matrix in which the spectrum radius can be adjusted by a relaxation coefficient.
  • S is a matrix having the inverse of each element of the vector s represented by the following equation as a diagonal component.
  • the matrix norm of the matrix on the left side of equation (17) can have any number of orders.
  • the arithmetic unit 5 of the X-ray CT apparatus 1 adds the elements of the matrix I- ⁇ (SB T DA + ⁇ SR) in the column direction, the maximum value in the row direction, operator norm
  • the arithmetic unit 5 of the X-ray CT apparatus 1 analytically calculates the relaxation coefficient ⁇ based on the imaging conditions.
  • the X-ray CT apparatus 1 can determine the relaxation coefficient ⁇ at high speed and easily compared to the iterative solution of the eigenvalue problem represented by the power method.
  • the relaxation coefficient ⁇ is determined based on the shooting conditions, the stability of the algorithm can be ensured as compared with the case where the relaxation coefficient ⁇ is determined empirically.
  • the scanner 2 acquires the actual projection data of the subject 10 based on the imaging conditions, and the operator instructs the image reconstruction via the input device 6, the arithmetic device 5 of the X-ray CT apparatus 1 Based on the imaging conditions input from 6, matrixes A, B, D, R, and R ′ in equation (17) are calculated (step 1).
  • the arithmetic unit 5 calculates each element of the matrix I ⁇ (SB T DA + ⁇ SR) (step 2).
  • the arithmetic unit 5 calculates the operator norm
  • the arithmetic unit 5 of the X-ray CT apparatus 1 analytically calculates the relaxation coefficient ⁇ based on the imaging conditions.
  • the arithmetic unit 5 of the X-ray CT apparatus 1 determines the relaxation coefficient ⁇ , it re-creates the tomographic image by performing successive approximation on the actual projection data using the update formula of the successive approximation method including the imaging conditions and the relaxation coefficient ⁇ .
  • the successive approximation method in the first embodiment applies matrices that are not transposed in the forward projection process and the backprojection process.
  • the update formula of the successive approximation method in the first embodiment is the formula (14).
  • the relaxation coefficient is calculated so as to satisfy Expression (17) defined based on the operator norm of the update matrix Q shown in Expression (16).
  • Expression (17) defined based on the operator norm of the update matrix Q shown in Expression (16).
  • (17) is an expression indicating that the operator norm of the update matrix Q is smaller than 1.
  • the arithmetic device 5 of the X-ray CT apparatus 1 calculates the value of the relaxation coefficient ⁇ applied to the equation (14) as the maximum value satisfying the conditional equation of the equation (17). That is, it is desirable for the arithmetic unit 5 to set as large a value as possible within the range that satisfies the conditional expression (17) as the value of the relaxation coefficient ⁇ applied to the expression (14).
  • the successive approximation method in the first embodiment performs back projection processing using view direction weights.
  • Arithmetic unit 5 calculates view direction weights with a weight of 0 except for a predetermined number of views on both sides around the view of the X-ray tube position closest to the target pixel in the body axis direction.
  • a back projection matrix B T that acts as a back projection process is calculated using the weight, and the successive approximation is performed on the actual projection data using the update formula of the formula (14) to which the back projection matrix B T is applied, and the subject 10 Reconstruct a tomographic image.
  • the arithmetic device 5 of the X-ray CT apparatus 1 calculates the relaxation coefficient ⁇ based on the imaging conditions, and uses the relaxation coefficient ⁇ calculated by itself to perform the actual projection. It was assumed that successive approximation to the data was performed.
  • the device that performs the calculation process of the relaxation coefficient ⁇ is not limited to the arithmetic device 5.
  • the arithmetic unit 5 of the X-ray CT apparatus 1 or other computer has a relaxation coefficient ⁇ shown in FIG. 8 for each of various imaging conditions (for example, imaging conditions used for the calculation of the system matrix A described above).
  • the calculation process is executed, and the calculation result is stored in the storage device 8 of the X-ray CT apparatus 1. That is, the storage device 8 stores the relaxation coefficient ⁇ for each photographing condition.
  • the arithmetic device 5 of the X-ray CT apparatus 1 uses the imaging condition input from the input device 6 as a search key, and the setting value of the relaxation coefficient ⁇ stored in the storage device 8 It is good also as what performs successive approximation with respect to real projection data using the setting value of the relaxation coefficient (alpha) acquired as a search result.
  • the image reconstruction process by the successive approximation method that performs the data expansion back projection process is executed using the relaxation coefficient ⁇ that is analytically calculated by the first embodiment.
  • Such an image reconstruction process is effective when it is necessary to image a long range in the body axis direction of the subject 10 within a predetermined time. For example, it is effective in cases such as emergency medicine where priority is given to shortening the shooting time over image quality. However, according to the present invention, it is possible to reduce the photographing time and to maintain the image quality at a certain level.
  • the projection data is expanded in the channel direction and the column direction, and the scanFOV 41 is expanded using virtual projection data (hereinafter referred to as “expansion projection data”). Is the purpose.
  • extended scanFOV the scanFOV41 that is expanded by the application of the data expansion type back projection process is referred to as “extended scanFOV”.
  • Fig. 9 (a) shows expansion in the channel direction.
  • the expansion channels 42a and 42b are projection data expanded in the channel direction.
  • the channel direction expansion scanFOV43 is obtained by expanding the scanFOV41 by applying the data expansion back projection process using the expansion channels 42a and 42b.
  • Fig. 9 (b) shows expansion in the column direction.
  • the extended column 44 is projection data extended in the column direction.
  • the column direction extended scanFOV 45 is obtained by extending the scanFOV 41 by using the extended column 44 and applying the data extended back projection process.
  • the arithmetic device 5 of the X-ray CT apparatus 1 performs a process for calculating the number of extended columns and the number of extended channels based on the imaging conditions input from the input device 6 (step 11).
  • the distance between the X-ray tube 11 and the rotation center is d sod .
  • the distance between the X-ray tube 11 and the center of the detector 12 is d sid .
  • FOM the FOM size is the size of the (Field Of Measurement) 51 and d phi.
  • the FOM size d phi as shown in FIG. 12, among the vertices of the reconstructed FOV53 (image region to be created) is twice the distance between the rotation center and the farthest vertex 54 from the center of rotation.
  • d dtc be the aperture width of the detection element in the body axis direction
  • r re be the number of columns of the detection element
  • r im be the number of columns to be expanded.
  • a straight line connecting the X-ray tube 11 and the detector center 52 is defined as the q axis. Since the length of one row of detector elements is constant, by calculating r im, expansion columns is obtained.
  • X-ray tube side plane a plane closer to the X-ray tube 11 (hereinafter referred to as “X-ray tube side plane”) out of two planes orthogonal to the q axis and in contact with the FOM 51.
  • d s the length of the X-ray tube side plane and the X-ray beam intersecting in the body axis direction.
  • F is the backprojection phase width, and in this embodiment, is a parameter for adjusting the trade-off between image quality degradation caused by estimation error of extended projection data and image quality degradation caused by loss of actual projection data. .
  • r im is calculated according to the following equation.
  • the arithmetic device 5 of the X-ray CT apparatus 1 calculates the number of columns r im to be expanded by the equation (20).
  • the arithmetic unit 5 may expand the channel so as to include the extended scanFOV43 in the imaging section, and can easily calculate the number of extended channels.
  • the arithmetic device 5 executes a process of creating extended projection data based on the number of extended columns and the number of extended channels calculated in step 11 (step 12).
  • extended projection data in the column direction is estimated by extrapolation from difference data between forward projection data and actual projection data.
  • extended projection data in the column direction is estimated by zero-order extrapolation using data that is located on the outermost side in the column direction among the data obtained as actual projection data.
  • the projection value ⁇ of the extended projection data is It is expressed as the following formula.
  • the computing device 5 of the X-ray CT apparatus 1 calculates the projection value ⁇ of the extended projection data by the equation (21).
  • the arithmetic device 5 executes the back projection process by combining the actual projection data of the subject 10 acquired by the scanner 2 and the extended projection data created in Step 12 (Step 13).
  • the present invention can also be applied to the axial scan.
  • the channel direction and the column direction are expanded, only one of the channel direction and the column direction may be expanded.
  • back projection processing is performed by using extended projection data obtained by extending real projection data in the channel direction and / or column direction.
  • the arithmetic device 5 of the X-ray CT apparatus 1 analytically calculates the relaxation coefficient ⁇ based on the imaging conditions input by the input device 6 as in the first embodiment. Then, the arithmetic device 5 calculates a conversion matrix P from the actual projection data to the extended projection data based on the imaging conditions input by the input device 6, and combines the actual projection data and the extended projection data using the conversion matrix P.
  • the successive approximation method by combining the first embodiment and the second embodiment uses the extended projection data obtained by extending the actual projection data in the channel direction and / or the column direction, and the view direction weight. By using this, back projection processing is performed.
  • the arithmetic device 5 of the X-ray CT apparatus 1 analytically calculates the relaxation coefficient ⁇ based on the imaging conditions input by the input device 6 as in the first embodiment.
  • the arithmetic device 5 calculates a conversion matrix P from the actual projection data to the extended projection data based on the imaging conditions input by the input device 6. Further, the computing device 5 calculates a view direction weight that gives a weight of 0 except for a predetermined number of views on both sides around the view of the X-ray tube position with the closest distance in the body axis direction to the target pixel, A back projection matrix B T to be operated as back projection processing is calculated using the view direction weight.
  • the arithmetic unit 5 calculates an extended backprojection matrix that acts as a backprojection process combining the actual projection data and the extended projection data using the transformation matrix P and the view direction weight, and this extended backprojection matrix is applied.
  • the tomographic image is reconstructed by performing successive approximation on the actual projection data and the extended projection data using the updating formula (14).

Landscapes

  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Apparatus For Radiation Diagnosis (AREA)
  • Health & Medical Sciences (AREA)
  • Mathematical Analysis (AREA)
  • Radiology & Medical Imaging (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Algebra (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Pulmonology (AREA)
  • Mathematical Physics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Chemical & Material Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Biochemistry (AREA)
  • General Health & Medical Sciences (AREA)
  • Immunology (AREA)
  • Pathology (AREA)

Abstract

　安定に収束し、かつ高速に実行可能な逐次近似法によって画像を再構成するX線CT装置等を提供するために、X線CT装置1の演算装置5は、入力装置6から入力される撮影条件に基づいて、行列A、B、D、R、R'を算出する(ステップ1)。次に、演算装置5は、行列I－α(SBTDA＋βSR)の各要素を算出する(ステップ2)。次に、演算装置5は、行列I－α(SBTDA＋βSR)の作用素ノルム||I－α(SBTDA＋βSR)||を算出する(ステップ3)。そして、演算装置5は、所定の条件式を満たすように、緩和係数αを決定する(ステップ4)。

Description

Ｘ線ＣＴ装置及び画像再構成方法
　本発明は、逐次近似法によって画像を再構成するX線CT装置等に関するものである。
　X線CT装置は、ファンビーム(扇形ビーム)もしくはコーンビーム(円錐形または角錐形ビーム)状のX線を被検体に照射し、被検体を透過したX線をX線検出器により計測し、多方向からの計測データを再構成することにより被検体の断層像を得る装置である。
　X線CT装置における画像再構成法は、解析法と逐次近似法に大別される。解析法は、投影切断面定理に基づき解析的に問題を解く方法である。逐次近似法は、投影データ取得に至る観測系を数学的にモデル化し、数学モデルに基づいて反復法により最良の画像を推定する方法である。
　両者を比較すると、解析法の利点としては、実投影データから直接的に再構成画像が得られるため、計算量が圧倒的に少ないことが挙げられる。一方、逐次近似法の利点としては、投影データの取得に至る物理的過程および実投影データに含まれる統計的な揺らぎをそれぞれ数学モデルや統計モデルとして考慮できるため、解析法において発生するアーチファクト(コーンビームアーチファクトなど)や画像上の量子ノイズを低減できることが挙げられる。
　従来、マルチスライスCTにおける画像再構成法としては、計算量が少ないことを採用理由として、解析法であるFeldkamp法、またはFeldkamp法を改良した手法が主に用いられている。しかし、近年のコンピュータの高性能化に伴い、逐次近似法も実用化が検討され始めている。
　逐次近似法は、画像の評価指標を事前に設定しておき、評価指標を数値化した評価値が最大値もしくは最小値をとるように画像を逐次更新する方法である。評価指標としては、更新の過程において画像を投影データに変換した順投影データと実投影データ間の矛盾や確率的な尤もらしさなどが用いられ、評価値を算出するための関数は評価関数と呼ばれる。
　非特許文献1では、罰則付き加重二乗誤差関数を評価関数として用いる逐次近似法が提案されている。これまでに提案されている手法としては、非特許文献1において提案されているように、順投影処理と逆投影処理とで互いに転置の関係にある行列を作用させることが一般的である。
　一方、少数ではあるが、順投影処理と逆投影処理とで互いに転置の関係ではない行列を作用させる逐次近似法も提案されている。非特許文献2では、ビュー方向重みを逆投影処理に適用して逐次更新を行う手法が提案されている。
　以降では、ビュー方向重みを使用する逆投影処理を「ビュー方向重み付き逆投影処理」と呼ぶことにする。
　ビュー方向重み付き逆投影処理自体は、非特許文献3に提案されており、解析法において使用されている技術である。ビュー方向重み付き逆投影処理は、次のような利点を持つ。
(1)投影データの冗長性を排除できる。
(2)時間分解能を向上できる。
　順投影処理と逆投影処理とで互いに転置の関係ではない行列を作用させる場合、逐次更新を行うための更新式には、逐次近似法の収束の速度と安定性に関わる緩和係数が含まれる。逐次近似法が安定に収束するためには、緩和係数を特定の範囲に設定する必要がある。
　非特許文献2では、緩和係数を経験的に決定すると記載されている。一方、非特許文献4では、べき乗法を用いて緩和係数を算出する手法が提案されている。ベキ乗法は、任意の行列の最大固有値を求める反復解法である。
　ところで、高速な寝台移動速度によってらせんスキャンを行った場合、列方向の実投影データが不足し、逐次近似法により画像が得られる領域が制限される問題がある。
　この問題に対して、特許文献1では、解析法において、実投影データを拡張して仮想列データおよび仮想チャネルデータを作成した後、逆投影処理を行う手法が記載されている。仮に、特許文献1の手法を逐次近似法に応用することができれば、領域の制限を緩和することができる。
　以降では、特許文献1の手法のように、実投影データを拡張して仮想列データおよび仮想チャネルデータを作成する処理を「データ拡張型逆投影処理」と呼ぶことにする。
特開2009-90139号公報
H. Erdogan et. al., "Ordered subsets algorithms for transmission tomography," Phys. Med. Biol., Vol.44, pp.2835-2851, 1999 J. Sunnegardh, "Combining analytical and iterative reconstruction in helicalcone-beam CT," Linkoping Studies in Science andTechnology Thesis No. 1301, 2007 S. Wesarg et. al., "Parker weights revisited," Med. Phys.Vol.29, No.3, pp372-378, March 2002 G. L. Zeng and G. T. Gullberg, "Unmatched projector/backprojector pairs in an iterativereconstruction algorithm," IEEE. Trans. Med. Imag,Vol.19, No.5, pp548-555, May 2000
　従来技術によって解決されておらず、本発明が解決しようとする課題は、以下の通りである。
　X線CT装置のスキャン撮影では、X線管とX線検出器を回転させながらデータを収集していくため、実投影データのビュー方向において撮影時間差が生じる。そのため、スキャン中に被検体が動いた場合、ビュー方向に位置情報が異なる実投影データが収集されることになる。このようなデータに対して、順投影処理と逆投影処理とで互いに転置の関係にある行列を用いて逐次近似法を行った場合、被検体の動きによる矛盾を反映した画像が最適な画像として再構成される。その結果、得られる画像には、動きに起因した画質劣化を伴うことから、何らかの補正方法が必要である。
　非特許文献2の手法のように、逐次近似法にてビュー方向重み付き逆投影処理を使用することによって、被検体の動きに起因する画質の劣化を低減することはできる。しかし、非特許文献2の手法では、操作者が経験的に緩和係数を決定することに起因する問題がある。
　逐次近似法が収束するための緩和係数の十分条件は、撮影条件としての寝台送り速度や撮影FOV(Field Of View)等に依存して変化する。これらの多数の撮影条件に応じて、操作者が経験的に緩和係数を決定することは非常に煩雑であり、時間がかかる。
　また、一般に、緩和係数を小さな値に設定すると、更新行列の性質により、撮影条件に因らず、逐次近似処理の収束条件は満たす。しかし、緩和係数をあまりに小さな値に設定してしまうと、逐次近似処理の収束が遅くなり、再構成される画像の画質も劣化する。
　また、非特許文献4の手法のように、ベキ乗法を使用して緩和係数を決定する場合、撮影条件に応じて緩和係数を自動的に算出するため、逐次近似法自体は安定に収束する。しかし、X線CT装置の画像再構成に使用する画像と実投影データの次元数が膨大であるため、非特許文献4の手法では、大規模な行列演算を逐次計算する必要が生じ、計算時間が増大する。
　更に、特許文献1の手法のデータ拡張型逆投影処理を、逐次近似法の逆投影処理に適用した場合、順投影処理と逆投影処理とで互いに転置の関係ではない行列を逆投影処理において作用させることになる。そうすると、非特許文献2及び非特許文献4と同様、緩和係数に最適な値を設定する必要がある。すなわち、緩和係数を最適な値に設定しなければ、計算時間が増大するという課題が残る。
　本発明は、前述した問題点に鑑みてなされたもので、第1の目的は、安定に収束し、かつ高速に実行可能な逐次近似法によって画像を再構成するX線CT装置等を提供することである。第2の目的は、体動のあるデータに対し逐次近似法を適用した場合の画質劣化を抑制できるX線CT装置等を提供することである。第3の目的は、寝台移動速度が高速のらせんスキャン時またはアキシャルスキャン時においてデータ欠損が生ずる撮影条件において、逐次近似法を適用した場合の画質劣化を抑制できるX線CT装置等を提供することである。
　前述した目的を達成するために第1の発明は、順投影処理と逆投影処理とで互いに転置ではない行列を作用させる逐次近似法によって被写体の断層像を再構成するX線CT装置であって、撮影条件に基づいて前記被写体の実投影データを取得する撮影部と、前記撮影条件、及び、計算の収束性を定める緩和係数を含む前記逐次近似法の更新式を用いて、前記実投影データに対する逐次近似を行うことによって前記断層像を再構成する演算部と、を具備し、前記緩和係数は、前記撮影条件に基づいて解析的に算出されることを特徴とするX線CT装置である。
　第2の発明は、順投影処理と逆投影処理とで互いに転置ではない行列を作用させる逐次近似法によって被写体の断層像を再構成する画像再構成方法であって、撮影条件に基づいて前記被写体の実投影データを取得するステップと、前記撮影条件、及び、計算の収束性を定める緩和係数を含む前記逐次近似法の更新式を用いて、前記実投影データに対する逐次近似を行うことによって前記断層像を再構成するステップと、を含み、前記緩和係数は、前記撮影条件に基づいて解析的に算出されることを特徴とする画像再構成方法である。
　本発明により、安定に収束し、かつ高速に実行可能な逐次近似法によって画像を再構成するX線CT装置等を提供することができる。また、体動のあるデータに対し逐次近似法を適用した場合の画質劣化を抑制できるX線CT装置等を提供することができる。また、寝台移動速度が高速のらせんスキャン時またはアキシャルスキャン時においてデータ欠損が生ずる撮影条件において、逐次近似法を適用した場合の画質劣化を抑制できるX線CT装置等を提供することができる。
X線CT装置1の全体外観図 X線CT装置1の構成図 アキシャルスキャンとらせんスキャンを説明する為の図 scanFOV41を説明する為の図 寝台移動速度の違いによるscanFOV41の違いを説明する為の図 注目ビューと対向位相ビューの相関を説明する為の図 ビュー方向重みを説明する為の図 緩和係数の算出処理を示すフローチャート 拡張scanFOVを説明する為の図 データ拡張型逆投影処理を示すフローチャート 拡張列数の算出処理を説明する為の図 FOMサイズを説明する為の図
　以下図面に基づいて、本発明の実施形態を詳細に説明する。最初に、図1～図3を参照しながら、全ての実施形態に共通するX線CT装置1の構成及びX線CT装置1による処理を説明する。
　図1に示すように、X線CT装置1は、X線管11や検出器12が搭載されるスキャナ2(撮影部)、被検体10を載置する寝台4、検出器12から得られるデータの処理を行う演算装置5(演算部)、マウス、トラックボール、キーボード、タッチパネルなどの入力装置6、及び再構成画像などを表示する表示装置7などを含む。
　操作者は、入力装置6を介して、撮影条件や再構成パラメータなどを入力する。撮影条件は、例えば、寝台送り速度、管電流、管電圧、スライス位置などである。また、再構成パラメータは、例えば、関心領域、再構成画像サイズ、逆投影位相幅、再構成フィルタ関数などである。
　図2に示すように、X線CT装置1は、大きく分けて、スキャナ2、操作ユニット3、寝台4から構成される。
　スキャナ2は、X線管11(X線発生装置)、検出器12、コリメータ13、駆動装置14、中央制御装置15、X線制御装置16、高電圧発生装置17、スキャナ制御装置18、寝台制御装置19、寝台移動計測装置20、コリメータ制御装置21、プリアンプ22、A/Dコンバータ23などから構成されている。
　中央制御装置15は、操作ユニット3における入力装置6から撮影条件や再構成パラメータを入力し、撮影に必要な制御信号を、コリメータ制御装置21、X線制御装置16、スキャナ制御装置18、寝台制御装置19に送信する。
　コリメータ制御装置21は、制御信号に基づいてコリメータ13の位置を制御する。
　撮影スタート信号を受けて撮影が開始されると、X線制御装置16は、制御信号に基づいて高電圧発生装置17を制御する。高電圧発生装置17は、X線管11(X線発生装置)に管電圧、管電流を印加する。X線管11では、印加された管電圧に応じたエネルギーの電子が陰極から放出され、放出された電子がターゲット(陽極)に衝突することによって電子エネルギーに応じたエネルギーのX線が被検体10に照射される。
　また、スキャナ制御装置18は、制御信号に基づいて駆動装置14を制御する。駆動装置14は、X線管11、検出器12、プリアンプ22等が搭載されているガントリ部を被検体10の周りに周回させる。
　寝台制御装置19は、制御信号に基づいて寝台4を制御する。
　X線管11から照射されるX線は、コリメータ13によって照射領域が制限され、被検体10内の各組織においてX線減弱係数に応じて吸収(減衰)され、被検体10を通過し、X線管11に対向する位置に配置された検出器12によって検出される。検出器12は、2次元方向(チャネル方向およびこれに直交する列方向)に配置された複数の検出素子によって構成される。各検出素子によって受光したX線は、投影データ(以降、「実投影データ」という。)に変換される。すなわち、検出器12によって検出されるX線は、電流に変換され、プリアンプ22によって増幅され、A/Dコンバータ23によってデジタルデータに変換され、LOG変換され、キャリブレーションが行われて実投影データとして演算装置5に入力される。
　このとき、互いに対向するX線管11と検出器12が、被検体10の周囲を回転するので、実投影データは、回転方向の離散的なX線管位置(および対応する検出器位置)において収集される。各々のX線管位置における実投影データの取得単位が、「ビュー」と呼ばれている。
　演算装置5は、再構成演算装置31、画像処理装置32等から構成される。また、入出力装置9は、入力装置6、表示装置7、記憶装置8(記憶部)等から構成される。
　再構成演算装置31は、実投影データを用いて画像再構成処理を行い、再構成画像を生成する。再構成演算装置31は、各ビューの実投影データに再構成フィルタを重畳してフィルタ補正投影データを生成し、フィルタ補正投影データに対して、ビュー方向に重み(以降、「ビュー方向重み」という。)を加重して逆投影処理を行うことによって、被検体10内部のX線減弱係数の分布図として非破壊的に断層像を画像化する。
　再構成演算装置31は、生成される再構成画像を記憶装置8に保存する。また、再構成演算装置31は、表示装置7にCT画像として再構成画像を表示する。あるいは、画像処理装置32が、記憶装置8に保存される再構成画像に対して画像処理を行い、画像処理された再構成画像を表示装置7にCT画像として表示する。
　X線CT装置1は、2次元方向に検出素子が配列された検出器12を用いるマルチスライスCT、検出素子が1列すなわち1次元方向(チャネル方向のみ)に配列された検出器12を用いるシングルスライスCTに大別される。マルチスライスCTでは、検出器12に合わせてX線源であるX線管11から円錐状、もしくは角錐状に広がるX線ビームが照射される。シングルスライスCTでは、X線管11から扇状に広がるX線ビームが照射される。通常、X線CT装置1による撮影では、ガントリ部が、寝台4に載置された被検体10の周りを周回しながら、X線の照射が行われる(但し、スキャノグラム撮影は除く。
)。
　図3(a)に示すように、撮影中に寝台4が固定され、X線管11が被検体10の周りを円軌道状に周回する撮影は、アキシャルスキャンなどと呼ばれる。また、図3(b)に示すように、寝台4が移動し、X線管11が被検体10の周りをらせん軌道状に周回する撮影は、らせんスキャンなどと呼ばれる。
　寝台制御装置20は、アキシャルスキャンの場合、寝台4を静止した状態とする。また、寝台制御装置20は、らせんスキャンの場合、入力装置6を介して入力される撮影条件としての寝台送りの速さに応じて、寝台4を体軸方向に平行移動させる。
　本発明の実施の形態に係るX線CT装置1は、例えば、マルチスライスCTである。また、X線CT装置1のスキャン方式は、例えば、ローテート－ローテート方式(第3世代)である。
　次に、図4、図5を参照しながら、各実施形態の前提となる画像再構成処理について説明する。具体的には、(1)罰則付き加重二乗誤差関数を評価関数として用いる逐次近似法、(2)ビュー方向重み付き逆投影処理を行う逐次近似法について説明する。
　最初に、罰則付き加重二乗誤差関数を評価関数として用いる逐次近似法について説明する。これは、例えば、非特許文献1において提案されており、その更新式は次式によって表される。
Figure JPOXMLDOC01-appb-I000001
　ここで、x(k)は、k回目の逐次更新における画像を表すベクトル(画像ベクトル)であり、ｙは実投影データを表すベクトルである。Aは、画像と投影データを対応付ける行列であり、前述の数学モデルを介して撮影系の特性を表すことから、システムマトリクスと呼ばれる。ATは、Aの転置行列を表す。
　Ax(k)は、画像ベクトルx(k)を投影データベクトルに変換する処理(順投影処理)に相当する。また、AT(・)は、括弧内の投影データベクトルを画像ベクトルに変換する処理(逆投影処理)に相当する。
　Dは、実投影データと順投影データの差分値に加重される重み係数を対角成分に持つ対角行列であり、CTの画像再構成においては検出フォトン数に応じた値を重み係数とする。
　βは、罰則項の強さを調整する任意のパラメータである。罰則項は、逐次近似法の適用による画像の高周波強調効果を抑制する役割を持つ。
　Sは、次式によって表されるベクトルsの各要素の逆数を対角成分に持つ行列である。
Figure JPOXMLDOC01-appb-I000002
　ここで、cは、画像ベクトルと等しい要素数を持ち、全要素に1を値として持つベクトルである。(1)式のRおよび(2)式のR'は、罰則項の1階微分および2階微分の線型作用素であり、ともに画像ベクトルから画像ベクトルへの変換行列である。
　RおよびR'のm行n列の要素は、それぞれ次式のように表される。
Figure JPOXMLDOC01-appb-I000003
　ここで、lmnはm番目とn番目の画素間の距離の逆数である。尚、以上の説明では、罰則項として2次関数を用いた場合を例として示している。
　(1)式に示す更新式は、順投影処理と逆投影処理とで互いに転置の関係にある行列を作用させている。
　次に、ビュー方向重み付き逆投影処理を行う逐次近似法について説明する。これは、例えば、非特許文献2にて提案されている。尚、ビュー方向重み付き逆投影処理自体は、非特許文献3にて提案されている。
　ビュー方向重み付き逆投影処理に対応する逆投影行列をBT(・)とおくと、非特許文献2の提案手法は、次式によって表される。
Figure JPOXMLDOC01-appb-I000004
　αは、逐次近似法の収束の速度と安定性に関わる緩和係数である。順投影処理と逆投影処理とで互いに転置ではない行列を作用させる場合(A≠Bの場合)、逐次近似法が安定に収束するためには、特定の値の範囲において緩和係数を設定する必要がある。
　このことを説明する為に、まず、(5)式を以下のように変形する。
Figure JPOXMLDOC01-appb-I000005
　本発明の実施の形態では、(7)式の行列Qを更新行列と呼ぶことにする。更新行列は、更新式に基づく1回の逐次更新処理当たりの画像の更新量と更新方向を決定する行列である。Iは、単位行列を表す。
　更新行列のスペクトル半径をρ(Q)とおくと、次に示す(8)式の条件の下、(6)式の画像ベクトルx(k+1)は、次に示す(9)式に収束することが知られている。尚、行列のスペクトル半径とは、行列の固有値の絶対値の最小上界のことである。
Figure JPOXMLDOC01-appb-I000006
　ここで、(・)-1は、括弧内の行列の逆行列を表す。
　(7)式の更新行列が、(8)式の条件を満たすように緩和係数αを設定することにより、逐次近似法が安定に収束することになる。
　尚、非特許文献4にて提案されている手法では、ベキ乗法を使用し、更新行列の最大固有値を算出することによって、(8)式に示すスペクトル半径ρ(Q)を算出し、緩和係数αを設定している。しかしながら、べき乗法は、問題を解析的に解く手法ではなく、大規模な行列演算を逐次計算する必要が生じ、計算時間が増大することなる。
　ところで、一般に、高速な寝台移動速度によってらせんスキャンを行う場合、列方向の実投影データが不足し、逐次近似法により画像が得られる領域が制限される問題がある。
この問題については、図4、図5を参照しながら説明する。
　図4(a)では、撮影断面における検出器12とX線管11の配置を示している。また、撮影断面におけるscanFOV41を示している。scanFOV41とは、投影切断面定理に基づき画像が結像する領域である。撮影断面におけるscanFOV41は、ほぼ円形となる。
　図4(b)では、体軸方向における検出器12とX線管11の配置を示している。また、体軸方向におけるscanFOV41を示している。体軸方向におけるscanFOV41は、寝台移動速度によって異なるが、3角形、4角形又は5角形等の多角形となる。
　図5では、寝台移動速度の違いによる体軸方向におけるscanFOV41の違いを示している。図5(a)と図5(b)を比較すると、寝台移動速度が速い程、scanFOV41の体軸方向の幅が狭くなることが分かる。scanFOV41の体軸方向の幅が所定値よりも狭い場合、ある断層像を再構成するために必要な180度分のビューが取得出来ていないことを意味し、ひいては、画像を作成できる領域が制限されることを意味する。
　同様の理由によって、アキシャルスキャン時の体軸方向の周辺視野においても、画像を作成できる領域が制限される。
　更に、撮影断面内においても、検出器12のチャネル方向において被検体10がはみ出す場合、同様の制限を受ける。
　以降では、前述の説明した内容を前提として、本発明の実施の形態を具体的に説明する。
＜第1の実施の形態＞
　第1の実施の形態では、画像再構成処理を実行する前に、逐次近似法を安定に収束させるための緩和係数の値を解析的に算出する。そして、ビュー方向重み付き逆投影処理を行う逐次近似法による画像再構成処理を実行する。
　ビュー方向重み付き逆投影処理は、下記に示すX線CT装置1のスキャンの特徴を利用している。
　(1)図6に示すように、半回転(回転方向に180度)ずつ撮影位相の異なるビューについて、対応するX線の透過経路がほぼ重なるチャネルおよび列が存在し、それらの実投影データは互いに高い相関性を持つ。図6(a)では、チャネル方向における注目ビューの透過経路及び対向位相ビューの透過経路を示している。図6(b)では、列方向における注目ビューの透過経路及び対向位相ビューの透過経路を示している。以降では、このような関係にあるデータを対向位相の実投影データと呼ぶことにする。
　(2)前述のように、X線CT装置1は、X線管11と検出器12を回転させながらデータを収集していくため、実投影データのビュー方向に撮影時間差が生ずる。撮影中の寝台移動速度が一定であれば、撮影時間差はX線管11の体軸方向位置と画像の体軸方向位置の距離に比例する。
　図7は、ビュー方向重みの概念図である。図7に示すように、X線CT装置1の演算装置5は、ビュー方向重みとして、注目する画素を含むスライスの体軸方向位置とX線管11の体軸方向位置の距離に応じて、対向位相の実投影データ間で加算して1になるような重みを加重する。更に、X線CT装置1の演算装置5は、ビュー方向重みとして、注目画素との体軸方向の距離が最も近いX線管11の体軸方向位置のビューを中心として両側の所定ビュー数分のビューを除き0の重みを付ける。以降では、0ではない重みが付けられるビューを逆投影ビューと呼び、逆投影ビューの数を逆投影ビュー数と呼ぶことにする。
　逆投影ビュー数は、らせんピッチ(撮影条件としての寝台送り速度及びガントリ回転速度によって定まる値)、及び検出器12の列方向の長さによって最適な値が定まる。
　このようにビュー方向重みを定義することによって、ビュー方向重みを使用しない場合と比較して時間分解能を向上させることができる。これは、ビュー方向重みによって、画像再構成処理に寄与する複数のビューの時相がほぼ一致するからである。
　ビュー方向重みの定義は、Feldkamp法と同様、全画素において共通の逆投影ビュー数を使用し、スライス毎に中心ビューを決定しても良い。また、ビュー方向重みの定義は、「特開2004－199163」のように、中心ビューを画素毎に決定しても良い。また、ビュー方向重みの定義は、各画素において逆投影ビュー数を決定しても良い。
　実投影データのベクトルのインデックスをi(i＝1、・・・、I)、画像ベクトルのインデックスをj(j＝1、・・・、J)とすると、ビュー方向重みは、実投影データのベクトルと画像ベクトルの要素毎に設定されるため、wijと定義される。i番目の投影データとj番目の画素において、システムマトリクスAの要素をaijとし、ビュー方向重みを乗算したシステムマトリクスBの要素をbijとおくと、bijは次式の通り表される。
Figure JPOXMLDOC01-appb-I000007
　このとき、ビュー方向重みつき逆投影処理を表す逆投影行列はBTとなる。行列BTを逐次近似法の逆投影処理に使用する場合、前述のように(8)式の条件下においてアルゴリズムが収束する。
　ここで、システムマトリクスAの算出例について説明する。
　システムマトリクスAは、前述の通り、画像と投影データを対応付ける行列である。
　X線CT装置1の演算装置5は、X線CT装置1の装置仕様や、入力装置6から入力される撮影条件に基づいて、各要素aijを算出する。システムマトリクスAの計算手法としては、例えば、「Joseph, P. M. (1982).An Improved Algorithm for Reprojecting Rays Through Pixel Images. IEEE Transactions on Medical Imaging MI-1, 192-196.」において提案されているJoseph法が挙げられる。
　システムマトリクスAの計算に用いられるX線CT装置1の装置仕様としては、例えば、以下が挙げられる。
・検出器12の列方向素子数
・検出器12のチャネル方向素子数
・検出器12の列方向1素子のサイズ
・検出器12のチャネル方向1素子のサイズ
・X線管11と回転中心との距離dsod
・X線管11と検出器12の中心(検出器中心)との距離dsid
　また、システムマトリクスAの計算に用いられる撮影条件としては、例えば、以下が挙げられる。尚、意味が分かり難い撮影条件については、設定内容を例示する。
・ビームピッチ(寝台送り速度を回転中心における列方向検出器サイズによって割った値)：「高速に相当する値」又は「低速に相当する値」など。
・撮影回転数
・撮影開始角度
・回転速度
・回転あたりの撮影ビュー数
・撮影方法：「アキシャルスキャン」又は「らせんスキャン」など。
・画素数(断面内)：横方向の画素数×縦方向の画素数
・画像スライス数
・画素サイズ(断面内)
・画像スライス厚
・画像中心位置(断面内)：「回転中心に一致」など。
・画像スライス方向中心位置：「撮影回転数2．5回転目に一致」など。
　前述の実投影データのベクトルの次元数Iは、装置仕様及び撮影条件から定まる。具体的には、I＝検出器12の列方向素子数×検出器12のチャネル方向素子数×回転あたりの撮影ビュー数×撮影回転数である。同様に、画像ベクトルの次元数Jも、撮影条件から定まる。具体的には、J＝横方向の画素数×縦方向の画素数×画像スライス数である。
　例えば、撮影条件の1つであるビームピッチが、「高速に相当する値」の場合と、「低速に相当する値」の場合とでは、システムマトリクスAが異なる。ひいては、逐次近似法が安定に収束する緩和係数αも異なる。
　ところで、前述の非特許文献4では、ベキ乗法により更新行列のスペクトル半径を算出するが、計算時間および計算に要するメモリ使用量が多いという問題がある。そこで、本発明におけるX線CT装置1では、更新行列のスペクトル半径を算出する代わりに、更新行列の作用素ノルムを算出する。
　一般に、行列の作用素ノルムは次のような性質を持つ。
　(1)任意のベクトルxに関して、ベクトルのノルムにp－ノルムを採用したときの任意の正方行列Cの作用素ノルム||C||は、次式の通り表される。
Figure JPOXMLDOC01-appb-I000008
　(2)いかなる作用素ノルム||C||に対しても次式が成り立つ。
Figure JPOXMLDOC01-appb-I000009
　上記の2つの性質を利用して、(8)式は、次式の通り置き換えることができる。
Figure JPOXMLDOC01-appb-I000010
　明らかに、更新行列の作用素ノルムが1より小さい場合((13)式を満たす場合)、更新行列のスペクトル半径は1より小さい((8)式を満たす)。
　以降では、非特許文献3にて提案されている手法に対し、ビュー方向重み付き逆投影処理および緩和係数を導入する例を挙げる。但し、本発明は、この例に限定されず、緩和係数によってスペクトル半径を調整可能な更新行列を更新式に含む逐次近似法であれば適用可能である。
　まず、(1)式にビュー方向重み付き逆投影処理の逆投影行列BTおよび緩和係数αを導入し、次式を得る。
Figure JPOXMLDOC01-appb-I000011
　ここで、Sは、次式にて表されるベクトルsの各要素の逆数を対角成分に持つ行列である。
Figure JPOXMLDOC01-appb-I000012
　さらに、(7)式と同様にして更新行列Qを導出すると、次式のようになる。
Figure JPOXMLDOC01-appb-I000013
　(16)式を(13)式に代入し、次式を得る。
Figure JPOXMLDOC01-appb-I000014
　(17)式を満たすように緩和係数αを決定することにより、(14)式に示す更新式が安定に収束する。
　ここで、(17)式の左辺における行列の作用素ノルムの次数は、いくつでも良い。
　例えば、(17)式において、1－ノルム(1次ノルム)を採用した場合、X線CT装置1の演算装置5は、行列I－α(SBTDA＋βSR)の要素を列方向に加算し、行方向における最大値を、作用素ノルム||I－α(SBTDA＋βSR)||の値とすれば良い。
　以上によって、X線CT装置1の演算装置5は、撮影条件に基づいて緩和係数αを解析的に算出する。これによって、X線CT装置1は、ベキ乗法に代表される固有値問題の反復解法と比較して、高速かつ簡便に緩和係数αを決定できる。
　また、撮影条件に基づいて緩和係数αが決定されることから、緩和係数αを経験的に決定する場合と比較して、アルゴリズムの安定性を確保できる。
　図8を参照しながら、緩和係数の算出処理の詳細について説明する。
　スキャナ2が撮影条件に基づいて被検体10の実投影データを取得し、操作者が入力装置6を介して画像再構成の指示を行うと、X線CT装置1の演算装置5は、入力装置6から入力される撮影条件に基づいて、(17)式における行列A、B、D、R、R'を算出する(ステップ1)。
　次に、演算装置5は、行列I－α(SBTDA＋βSR)の各要素を算出する(ステップ2)。
　次に、演算装置5は、行列I－α(SBTDA＋βSR)の作用素ノルム||I－α(SBTDA＋βSR)||を算出する(ステップ3)。
　そして、演算装置5は、(17)式(所定の条件式)を満たすように、緩和係数αを決定する(ステップ4)。
　このように、X線CT装置1の演算装置5は、撮影条件に基づいて解析的に緩和係数αを算出する。
　X線CT装置1の演算装置5は、緩和係数αを決定すると、撮影条件及び緩和係数αを含む逐次近似法の更新式を用いて、実投影データに対する逐次近似を行うことによって断層像を再構成する。前述の通り、第1の実施の形態における逐次近似法は、順投影処理と逆投影処理とで互いに転置ではない行列を作用させるものである。そして、第1の実施の形態における逐次近似法の更新式は、(14)式である。
　第1の実施の形態では、緩和係数は、(16)式に示す更新行列Qの作用素ノルムに基づいて定義される(17)式を満たすように算出される。これによって、(14)式の更新式に基づく逐次更新処理は、安定に収束する。
　(17)式は、更新行列Qの作用素ノルムが1より小さいことを示す式である。ここで、X線CT装置1の演算装置5は、(14)式に適用される緩和係数αの値を、(17)式の条件式を満たす最大値として算出することが望ましい。すなわち、演算装置5は、(14)式に適用される緩和係数αの値として、(17)式の条件式を満たす範囲内において、出来る限り大きい値を設定することが望ましい。
　これによって、逐次更新処理の反復回数が少なくなり、計算時間が短縮されるとともに、再構成される画像の画質も向上する。
　また、第1の実施の形態における逐次近似法は、ビュー方向重みを用いる逆投影処理を行うものである。演算装置5は、注目画素との体軸方向の距離が最も近いX線管位置のビューを中心として両側に所定ビュー数分のビューを除き0の重みを付けるビュー方向重みを算出し、ビュー方向重みを用いて逆投影処理として作用させる逆投影行列BTを算出し、逆投影行列BTが適用される(14)式の更新式を用いて実投影データに対する逐次近似を行い、被検体10の断層像を再構成する。
　これによって、投影データの冗長性を排除できるとともに、時間分解能を向上できる。ひいては、体動のあるデータに対し逐次近似法を適用した場合の画質劣化を抑制できる。
　前述の説明では、X線CT装置1の演算装置5は、スキャナ2が実投影データを取得すると、撮影条件に基づいて緩和係数αを算出し、自らが算出する緩和係数αを用いて実投影データに対する逐次近似を行うものとした。
　しかしながら、本発明では、緩和係数αの算出処理を実行する装置は、演算装置5に限定されない。
　例えば、X線CT装置1の演算装置5、又は、その他のコンピュータが、様々な撮影条件(例えば、前述したシステムマトリクスAの計算に用いられる撮影条件)ごとに、図8に示す緩和係数αの算出処理を実行し、算出結果をX線CT装置1の記憶装置8に記憶させる。すなわち、記憶装置8は、撮影条件ごとに緩和係数αを記憶する。
　そして、X線CT装置1の演算装置5は、スキャナ2が実投影データを取得すると、入力装置6から入力される撮影条件を検索キーとして、記憶装置8に記憶される緩和係数αの設定値を検索し、検索結果として取得される緩和係数αの設定値を用いて、実投影データに対する逐次近似を行うものとしても良い。
＜第2の実施の形態＞
　第2の実施の形態では、第1の実施の形態によって解析的に算出される緩和係数αを用いて、データ拡張型逆投影処理を行う逐次近似法による画像再構成処理を実行する。
　このような画像再構成処理は、所定時間内に、被検体10の体軸方向に長い範囲を撮影する必要がある場合に有効である。例えば、救急医療などの場面において、画像の画質よりも、撮影時間の短縮を優先する場合に有効である。但し、本発明では、撮影時間の短縮を実現するとともに、画像の画質も一定レベルを維持することが可能となる。
　データ拡張型逆投影処理は、図9に示すように、チャネル方向および列方向に投影データを拡張し、仮想的な投影データ(以降、「拡張投影データ」という。)を用いて、scanFOV41を拡げることが目的である。以降では、データ拡張型逆投影処理の適用により拡張されるscanFOV41を「拡張scanFOV」と呼ぶことにする。
　図9(a)では、チャネル方向における拡張を示している。拡張チャネル42a及び42bは、チャネル方向に拡張される投影データである。また、チャネル方向拡張scanFOV43は、拡張チャネル42a及び42bを用いて、データ拡張型逆投影処理の適用によりscanFOV41が拡張されたものである。
　図9(b)では、列方向における拡張を示している。拡張列44は、列方向に拡張される投影データである。また、列方向拡張scanFOV45は、拡張列44を用いて、データ拡張型逆投影処理の適用によりscanFOV41が拡張されたものである。
　図10を参照しながら、データ拡張型逆投影処理の詳細について説明する。
　図10に示すように、X線CT装置1の演算装置5は、入力装置6から入力される撮影条件に基づいて、拡張列数および拡張チャネル数の算出処理を実行する(ステップ11)。
　図11、図12を参照しながら、演算装置5が拡張列数を算出する例を説明する。但し、本発明はこの例に限定されるわけではなく、操作者が任意に拡張列数を設定してもよい。また、アキシャルスキャンの場合、拡張列数の算出処理は、例えば、特許文献1と同様に行うことができる。
　図11に示すように、X線管11と回転中心との距離をdsodとする。また、X線管11と検出器12の中心(検出器中心52)との距離をdsidとする。また、FOM(Field Of Measurement)51の大きさであるFOMサイズをdφとする。
　ここで、FOMサイズdφは、図12に示すように、再構成FOV53(作成したい画像領域)の頂点のうち、回転中心から最も遠い頂点54と回転中心との距離の2倍である。
　また、体軸方向における検出素子の開口幅をddtc、検出素子の列数をrre、拡張する列数をrimとする。X線管11と検出器中心52を結ぶ直線をq軸とする。検出素子の1列分の長さは一定であるから、rimを算出することによって、拡張列数が求まる。
　ここで、q軸と直交し、かつFOM51に接する2つの平面のうち、X線管11に近い方の平面(以下、「X線管側平面」という。)に注目する。体軸方向において、X線管側平面とX線ビームが交差する長さをdsとすると、dsは、次式の通り表される。
Figure JPOXMLDOC01-appb-I000015
　また、ファン角をΦ、寝台送りの速さをκ(単位は、例えば、1回転当たりの移動距離[mm/rot])とし、観測系がπ＋Φ[rad]だけ回転したとき、体軸方向における寝台4に対する観測系の相対移動量dtは、次式の通り表される。
Figure JPOXMLDOC01-appb-I000016
　ここで、Fは逆投影位相幅であり、本実施形態においては拡張投影データの推定誤差に起因する画質劣化と実投影データの欠損に起因する画質劣化のトレードオフを調整するためのパラメータである。
　(18)式のdsと(19)式のdtが等しい場合、FOM51は、scanFOV41の内部に含まれるため、データ欠損の無い撮影条件となる。そこで、次式に従ってrimを算出する。
Figure JPOXMLDOC01-appb-I000017
　X線CT装置1の演算装置5は、(20)式によって拡張する列数rimを算出する。
　尚、チャネル方向については、演算装置5は、撮影断面内の拡張scanFOV43を内包するようにチャネルを拡張すればよく、拡張チャネル数を容易に算出できる。
　図10の説明に戻る。
　次に、演算装置5は、ステップ11によって算出される拡張列数および拡張チャネル数に基づいて、拡張投影データの作成処理を実行する(ステップ12)。
　ステップ12における拡張投影データの作成処理について説明する。
　(20)式に基づき算出される拡張列について、インデックスをr(r＝1、・・・、rim)とする。また、拡張投影データの投影値をΨ(Ψr：r＝1、・・・、rim)とする。
　本実施の形態では、順投影データと実投影データの差分データから、外挿によって列方向の拡張投影データを推定する。ここでは、最も単純な場合として、実投影データとして得られるデータの中で列方向に最も外側に位置するデータを用いる0次外挿によって列方向の拡張投影データを推定する。
　実投影データから拡張投影データへの変換行列をP(pri：r＝1、・・・、rim、i＝1、・・・、I)とすると、拡張投影データの投影値Ψは、次式の通り表される。
Figure JPOXMLDOC01-appb-I000018
　X線CT装置1の演算装置5は、(21)式によって拡張投影データの投影値Ψを算出する。
　そして、演算装置5は、スキャナ2によって取得される被検体10の実投影データ、及び、ステップ12において作成される拡張投影データを組み合わせて、逆投影処理を実行する(ステップ13)。
　ステップ13における実投影データと拡張投影データを組み合わせた逆投影処理について説明する。
　ステップ12において算出されるΨによって定義される逆投影行列をZ(zjr：j＝1、・・・、J、r＝1、・・・、rim)とする。これによって、拡張投影データの逆投影処理は、ZP(・)と表される。
　また、実投影データの逆投影処理は、AT(・)と表される。
　従って、実投影データと拡張投影データを組み合わせた逆投影処理を示す拡張逆投影行列をBTとすると、次式の通り表される。
Figure JPOXMLDOC01-appb-I000019
　以上から、データ拡張型逆投影処理を行う逐次近似法は、(14)式に、(22)式を代入することによって、定式化することができる。
　以上の説明では、らせんスキャンへの適用について説明したが、本発明は、アキシャルスキャンについても適用できる。また、チャネル方向および列方向について拡張したが、チャネル方向、列方向のどちらか一方のみを拡張してもよい。
　第2の実施の形態における逐次近似法は、実投影データをチャネル方向及び/又は列方向に拡張することによって得られる拡張投影データを用いることによって逆投影処理を行うものである。
　X線CT装置1の演算装置5は、第1の実施の形態と同様に、入力装置6によって入力される撮影条件に基づいて緩和係数αを解析的に算出する。そして、演算装置5は、入力装置6によって入力される撮影条件に基づいて実投影データから拡張投影データへの変換行列Pを算出し、変換行列Pを用いて実投影データ及び拡張投影データを組み合わせた逆投影処理として作用させる拡張逆投影行列BT(＝AT＋ZP)を算出し、拡張逆投影行列BTが適用される(14)式を更新式として、実投影データ及び拡張投影データに対する逐次近似を行うことによって断層像を再構成する。
　これによって、寝台移動速度が高速のらせんスキャン時またはアキシャルスキャン時においてデータ欠損が生ずる撮影条件において、逐次近似法を適用した場合の画質劣化を抑制できる。
　また、第1の実施の形態および第2の実施の形態を組み合わせることも可能である。
　第1の実施の形態および第2の実施の形態を組み合わせることによる逐次近似法は、実投影データをチャネル方向及び/又は列方向に拡張することによって得られる拡張投影データ、及び、ビュー方向重みを用いることによって逆投影処理を行うものである。
　まず、X線CT装置1の演算装置5は、第1の実施の形態と同様に、入力装置6によって入力される撮影条件に基づいて緩和係数αを解析的に算出する。
　次に、演算装置5は、入力装置6によって入力される撮影条件に基づいて実投影データから拡張投影データへの変換行列Pを算出する。また、演算装置5は、注目画素との体軸方向の距離が最も近いX線管位置のビューを中心として両側に所定ビュー数分のビューを除き0の重みを付けるビュー方向重みを算出し、ビュー方向重みを用いて逆投影処理として作用させる逆投影行列BTを算出する。
　そして、演算装置5は、変換行列P及びビュー方向重みを用いて実投影データ及び拡張投影データを組み合わせた逆投影処理として作用させる拡張逆投影行列を算出し、この拡張逆投影行列が適用される(14)式の更新式を用いて実投影データ及び拡張投影データに対する逐次近似を行うことによって断層像を再構成する。
　これによって、投影データの冗長性を排除できるとともに、時間分解能を向上できる。
ひいては、体動のあるデータに対し逐次近似法を適用した場合の画質劣化を抑制できる。
また、寝台移動速度が高速のらせんスキャン時またはアキシャルスキャン時においてデータ欠損が生ずる撮影条件において、逐次近似法を適用した場合の画質劣化を抑制できる。
　以上、本発明の様々な実施形態に関する記述から、本発明の目的が達成されることは明らかである。本発明を詳細にわたって記述すると共に図示しているが、これらは説明及び例示のみを意図したものであって、これらに限定されるものではない。また、本発明の要旨は、特許請求の範囲によってのみ限定されるものとする。
　1　X線CT装置、　2　スキャナ、　3　操作ユニット、　4　寝台、　5　演算装置、　6　入力装置、　7　表示装置、　8　記憶装置、　10　被検体、　11　X線管、　12　検出器、　41　scanFOV、　42a、42b　拡張チャネル、　43　チャネル方向拡張scanFOV、　44　拡張列、　45　列方向拡張scanFOV、　51　FOM、　52　検出器中心、　53　再構成FOV

Claims (10)

  1. 　順投影処理と逆投影処理とで互いに転置ではない行列を作用させる逐次近似法によって被写体の断層像を再構成するX線CT装置であって、
    　撮影条件に基づいて前記被写体の実投影データを取得する撮影部と、
    　前記撮影条件、及び、計算の収束性を定める緩和係数を含む前記逐次近似法の更新式を用いて、前記実投影データに対する逐次近似を行うことによって前記断層像を再構成する演算部と、
    　を具備し、
    　前記緩和係数は、前記撮影条件に基づいて解析的に算出されることを特徴とするX線CT装置。
  2. 　前記緩和係数は、前記更新式の更新量及び更新方向を定める更新行列の作用素ノルムに基づいて定義される所定の条件式を満たすように算出されることを特徴とする請求項1に記載のX線CT装置。
  3. 　前記所定の条件式は、前記更新行列の作用素ノルムが1より小さいことを示す式であることを特徴とする請求項2に記載のX線CT装置。
  4. 　前記緩和係数は、前記所定の条件式を満たす最大値として算出されることを特徴とする請求項3に記載のX線CT装置。
  5. 　前記逐次近似法は、前記実投影データをチャネル方向及び／又は列方向に拡張することによって得られる拡張投影データを用いることによって逆投影処理を行うものであり、
    　前記演算部は、前記撮影条件に基づいて前記実投影データから前記拡張投影データへの変換行列を算出し、前記変換行列を用いて前記実投影データ及び前記拡張投影データを組み合わせた逆投影処理として作用させる拡張逆投影行列を算出し、前記拡張逆投影行列が適用される前記更新式を用いて前記実投影データ及び前記拡張投影データに対する逐次近似を行うことによって前記断層像を再構成することを特徴とする請求項1に記載のX線CT装置。
  6. 　前記逐次近似法は、ビュー方向重みを用いる逆投影処理を行うものであり、
    　前記演算部は、注目画素との体軸方向の距離が最も近いX線管位置のビューを中心として両側に所定ビュー数分のビューを除き0の重みを付ける前記ビュー方向重みを算出し、前記ビュー方向重みを用いて逆投影処理として作用させる逆投影行列を算出し、前記逆投影行列が適用される前記更新式を用いて前記実投影データに対する逐次近似を行うことによって前記断層像を再構成することを特徴とする請求項1に記載のX線CT装置。
  7. 　前記逐次近似法は、前記実投影データをチャネル方向及び／又は列方向に拡張することによって得られる拡張投影データ、及び、ビュー方向重みを用いることによって逆投影処理を行うものであり、
    　前記演算部は、前記撮影条件に基づいて前記実投影データから前記拡張投影データへの変換行列を算出し、注目画素との体軸方向の距離が最も近いX線管位置のビューを中心として両側の所定ビュー数分のビューを除き0の重みを付ける前記ビュー方向重みを算出し、前記変換行列及び前記ビュー方向重みを用いて前記実投影データ及び前記拡張投影データを組み合わせた逆投影処理として作用させる拡張逆投影行列を算出し、前記拡張逆投影行列が適用される前記更新式を用いて前記実投影データ及び前記拡張投影データに対する逐次近似を行うことによって前記断層像を再構成することを特徴とする請求項1に記載のX線CT装置。
  8. 　前記演算部は、前記撮影部が前記実投影データを取得すると、前記撮影条件に基づいて前記緩和係数を算出し、自らが算出する緩和係数を用いて前記実投影データに対する逐次近似を行うことを特徴とする請求項1に記載のX線CT装置。
  9. 　前記撮影条件ごとに前記緩和係数を記憶する記憶部を更に具備し、
    　前記演算部は、前記撮影部が前記実投影データを取得すると、前記撮影条件に基づいて前記記憶部に記憶される前記緩和係数を検索し、検索される前記緩和係数を用いて前記実投影データに対する逐次近似を行うことを特徴とする請求項1に記載のX線CT装置。
  10. 　順投影処理と逆投影処理とで互いに転置ではない行列を作用させる逐次近似法によって被写体の断層像を再構成する画像再構成方法であって、
    　撮影条件に基づいて前記被写体の実投影データを取得するステップと、
    　前記撮影条件、及び、計算の収束性を定める緩和係数を含む前記逐次近似法の更新式を用いて、前記実投影データに対する逐次近似を行うことによって前記断層像を再構成するステップと、
    　を含み、
    　前記緩和係数は、前記撮影条件に基づいて解析的に算出されることを特徴とする画像再構成方法。
PCT/JP2011/073473 2010-10-14 2011-10-13 Ｘ線ｃｔ装置及び画像再構成方法 Ceased WO2012050149A1 (ja)

Priority Applications (3)

Application Number Priority Date Filing Date Title
CN201180049466.5A CN103153192B (zh) 2010-10-14 2011-10-13 X射线ct装置以及图像再构成方法
US13/824,697 US9123156B2 (en) 2010-10-14 2011-10-13 X-ray CT apparatus and image reconstruction method
JP2012538707A JP5828841B2 (ja) 2010-10-14 2011-10-13 Ｘ線ｃｔ装置及び画像再構成方法

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2010231840 2010-10-14
JP2010-231840 2010-10-14

Publications (1)

Publication Number Publication Date
WO2012050149A1 true WO2012050149A1 (ja) 2012-04-19

Family

ID=45938373

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2011/073473 Ceased WO2012050149A1 (ja) 2010-10-14 2011-10-13 Ｘ線ｃｔ装置及び画像再構成方法

Country Status (4)

Country Link
US (1) US9123156B2 (ja)
JP (1) JP5828841B2 (ja)
CN (1) CN103153192B (ja)
WO (1) WO2012050149A1 (ja)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104902818A (zh) * 2013-02-05 2015-09-09 株式会社日立医疗器械 X射线ct装置以及图像重构方法
CN105188543A (zh) * 2013-04-08 2015-12-23 株式会社日立医疗器械 X射线ct装置、重构运算装置以及重构运算方法
JP2015231493A (ja) * 2014-06-10 2015-12-24 株式会社東芝 Ｘ線診断装置
WO2017006764A1 (ja) * 2015-07-08 2017-01-12 株式会社日立製作所 画像演算装置、画像演算方法、および、断層画像撮影装置

Families Citing this family (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP2661736B1 (en) * 2011-01-06 2014-04-30 Koninklijke Philips N.V. Imaging system for imaging an object
CN103501702B (zh) * 2011-04-28 2015-09-02 株式会社日立医疗器械 医用图像处理装置、医用图像处理方法
CN103619259B (zh) * 2011-07-08 2015-12-09 株式会社日立医疗器械 图像重构装置及图像重构方法
DE102012204019B4 (de) * 2012-03-14 2018-02-08 Siemens Healthcare Gmbh Verfahren zur Reduzierung von Bewegungsartefakten
WO2014084291A1 (ja) * 2012-11-30 2014-06-05 株式会社 日立メディコ Ｘ線ｃｔ装置及びその断層画像撮影方法
US10192328B2 (en) * 2013-10-24 2019-01-29 Toshiba Medical Systems Corporation Method for statistical weights design in iterative reconstruction algorithms
JP6256608B2 (ja) * 2014-07-04 2018-01-10 株式会社島津製作所 画像再構成処理方法
EP3428627B1 (en) * 2016-03-11 2021-06-09 Shimadzu Corporation Image reconstruction method, image reconstruction program, and tomography device equipped with same
WO2018115025A1 (en) * 2016-12-21 2018-06-28 Koninklijke Philips N.V. Redundancy weighting for short scan off-center detector x-ray tomography
US10789738B2 (en) * 2017-11-03 2020-09-29 The University Of Chicago Method and apparatus to reduce artifacts in a computed-tomography (CT) image by iterative reconstruction (IR) using a cost function with a de-emphasis operator
JP7093075B2 (ja) * 2018-04-09 2022-06-29 東芝エネルギーシステムズ株式会社 医用画像処理装置、医用画像処理方法、およびプログラム
CN110084759B (zh) * 2019-04-23 2020-06-09 闽南师范大学 一种图像填补方法、终端设备及存储介质
CN114581545A (zh) 2020-12-01 2022-06-03 西安大医集团股份有限公司 医学图像重建方法、装置、计算机设备及存储介质

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2001509400A (ja) * 1997-07-01 2001-07-24 アナロジック コーポレーション 反復円錐ビームｃｔ画像再構成
JP2002065663A (ja) * 2000-05-17 2002-03-05 Ge Medical Systems Global Technology Co Llc イメージング・データを逆畳み込みするための方法及び装置
JP2006105975A (ja) * 2004-09-30 2006-04-20 General Electric Co <Ge> 予備補正を備えたｃｔ再構成方法及びシステム
WO2010016425A1 (ja) * 2008-08-07 2010-02-11 株式会社 日立メディコ Ｘ線ｃｔ画像形成方法及びそれを用いたｘ線ｃｔ装置
JP2010136958A (ja) * 2008-12-13 2010-06-24 Univ Of Tokushima Ｃｔ装置、ｃｔ装置における画像再構成方法、及び電子回路部品

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5414623A (en) * 1992-05-08 1995-05-09 Iowa State University Research Foundation Optoelectronic system for implementation of iterative computer tomography algorithms
US5909476A (en) * 1997-09-22 1999-06-01 University Of Iowa Research Foundation Iterative process for reconstructing cone-beam tomographic images
US6101236A (en) * 1998-10-02 2000-08-08 University Of Iowa Research Foundation Iterative method and apparatus for x-ray computed tomographic fluoroscopy

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2001509400A (ja) * 1997-07-01 2001-07-24 アナロジック コーポレーション 反復円錐ビームｃｔ画像再構成
JP2002065663A (ja) * 2000-05-17 2002-03-05 Ge Medical Systems Global Technology Co Llc イメージング・データを逆畳み込みするための方法及び装置
JP2006105975A (ja) * 2004-09-30 2006-04-20 General Electric Co <Ge> 予備補正を備えたｃｔ再構成方法及びシステム
WO2010016425A1 (ja) * 2008-08-07 2010-02-11 株式会社 日立メディコ Ｘ線ｃｔ画像形成方法及びそれを用いたｘ線ｃｔ装置
JP2010136958A (ja) * 2008-12-13 2010-06-24 Univ Of Tokushima Ｃｔ装置、ｃｔ装置における画像再構成方法、及び電子回路部品

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104902818A (zh) * 2013-02-05 2015-09-09 株式会社日立医疗器械 X射线ct装置以及图像重构方法
CN105188543A (zh) * 2013-04-08 2015-12-23 株式会社日立医疗器械 X射线ct装置、重构运算装置以及重构运算方法
CN105188543B (zh) * 2013-04-08 2017-11-14 株式会社日立制作所 X射线ct装置、重构运算装置以及重构运算方法
JP2015231493A (ja) * 2014-06-10 2015-12-24 株式会社東芝 Ｘ線診断装置
WO2017006764A1 (ja) * 2015-07-08 2017-01-12 株式会社日立製作所 画像演算装置、画像演算方法、および、断層画像撮影装置
JPWO2017006764A1 (ja) * 2015-07-08 2018-04-05 株式会社日立製作所 画像演算装置、画像演算方法、および、断層画像撮影装置
US10410343B2 (en) 2015-07-08 2019-09-10 Hitachi, Ltd. Image computing device, image computing method, and tomograph

Also Published As

Publication number Publication date
JP5828841B2 (ja) 2015-12-09
US20130177132A1 (en) 2013-07-11
CN103153192B (zh) 2015-09-23
CN103153192A (zh) 2013-06-12
US9123156B2 (en) 2015-09-01
JPWO2012050149A1 (ja) 2014-02-24

Similar Documents

Publication Publication Date Title
JP5828841B2 (ja) Ｘ線ｃｔ装置及び画像再構成方法
JP5978429B2 (ja) 医用画像処理装置、医用画像処理方法
US9050003B2 (en) Reconstruction computing device, reconstruction computing method, and X-ray CT apparatus
EP2349008B1 (en) Method for prior image constrained image reconstruction in cardiac cone beam computed tomography
JP5968316B2 (ja) 画像再構成装置及び画像再構成方法
CN101478919B (zh) 改善图像分辨率的系统和方法
JP6312401B2 (ja) 画像処理装置、画像処理方法、及びプログラム
US9662084B2 (en) Method and apparatus for iteratively reconstructing tomographic images from electrocardiographic-gated projection data
US9025838B2 (en) Apparatus and method for hybrid reconstruction of an object from projection data
JP2017189612A (ja) 放射線画像診断装置及び医用画像処理装置
US10565744B2 (en) Method and apparatus for processing a medical image to reduce motion artifacts
US9619902B2 (en) Filter for tomographic reconstructions
JP7062444B2 (ja) 医用画像生成装置及び医用画像生成方法
US20100246918A1 (en) Iterative extra-focal radiation correction in the reconstruction of ct images
JP2016152916A (ja) Ｘ線コンピュータ断層撮像装置及び医用画像処理装置
CN103190927A (zh) 确定运动场和产生运动补偿的ct 图像数据组的方法和系统
US9508164B2 (en) Fast iterative image reconstruction method for 3D computed tomography
US10049468B2 (en) Image reconstruction for computed tomography
US20070253529A1 (en) Systems and methods for digital volumetric laminar tomography
WO2016132880A1 (ja) 演算装置、ｘ線ｃｔ装置、及び画像再構成方法
US11375964B2 (en) Acquisition method, acquisition device, and control program for tomographic image data by means of angular offset
Hou et al. A compressed sensing approach to low-radiation CT reconstruction
Cierniak et al. Ultrafast iterative model-based statistical 3D reconstruction algorithm for x-ray computed tomography

Legal Events

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

Ref document number: 201180049466.5

Country of ref document: CN

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

Ref document number: 11832581

Country of ref document: EP

Kind code of ref document: A1

ENP Entry into the national phase

Ref document number: 2012538707

Country of ref document: JP

Kind code of ref document: A

WWE Wipo information: entry into national phase

Ref document number: 13824697

Country of ref document: US

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 11832581

Country of ref document: EP

Kind code of ref document: A1