WO2012147471A1 - 医用画像処理装置、医用画像処理方法 - Google Patents

医用画像処理装置、医用画像処理方法 Download PDF

Info

Publication number
WO2012147471A1
WO2012147471A1 PCT/JP2012/059136 JP2012059136W WO2012147471A1 WO 2012147471 A1 WO2012147471 A1 WO 2012147471A1 JP 2012059136 W JP2012059136 W JP 2012059136W WO 2012147471 A1 WO2012147471 A1 WO 2012147471A1
Authority
WO
WIPO (PCT)
Prior art keywords
penalty term
weight
true subset
projection data
term weight
Prior art date
Application number
PCT/JP2012/059136
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 US14/112,276 priority Critical patent/US9123098B2/en
Priority to CN201280020678.5A priority patent/CN103501702B/zh
Priority to JP2013511984A priority patent/JP5978429B2/ja
Publication of WO2012147471A1 publication Critical patent/WO2012147471A1/ja

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/0002—Inspection of images, e.g. flaw detection
    • G06T7/0012—Biomedical image inspection
    • 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
    • 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
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/70—Denoising; Smoothing
    • 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/583—Calibration using calibration phantoms
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10072—Tomographic images
    • G06T2207/10081—Computed x-ray tomography [CT]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20172—Image enhancement details
    • G06T2207/20182—Noise reduction or smoothing in the temporal domain; Spatio-temporal filtering
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/30—Subject of image; Context of image processing
    • G06T2207/30004—Biomedical image processing
    • 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

Definitions

  • the present invention relates to a medical image processing apparatus or the like that performs correction processing of projection data and / or reconstruction processing of a reconstructed image by a successive approximation method including weights and penalty terms according to the output of a detector in an evaluation function. is there.
  • the X-ray CT apparatus has an X-ray detector in which detection elements are arranged in the channel direction and the column direction.
  • the projection data becomes a position of the discrete X-ray tube (opposing X It can also be said that the position of the line detector.)
  • a unit for acquiring projection data at each X-ray tube position is referred to as a “view”.
  • the arithmetic unit obtains the filter corrected projection data by superimposing the reconstruction filter in the channel direction of the projection data for each view and each column, and in the view direction with respect to the filter corrected projection data.
  • a tomographic image is imaged non-destructively as a distribution map of the X-ray attenuation coefficient inside the subject.
  • an image reconstruction method a method for imaging a tomographic image from projection data
  • an image obtained by the image reconstruction method is referred to as a CT image.
  • Image reconstruction methods are roughly divided into analysis methods and successive approximation methods.
  • 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 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 projection data can be considered as mathematical models and statistical models, respectively. Artifacts) and quantum noise on the image can be reduced.
  • the Feldkamp method which is an analysis method, or an improved method of the Feldkamp method has been mainly used because of the small amount of calculation.
  • the practical use of the successive approximation method has begun to be studied.
  • a technique for improving image quality by a successive approximation method has been proposed.
  • the successive approximation method is roughly divided into the following two types depending on the estimated variable.
  • (1) A method to construct an evaluation function using the projection value of projection data as a variable
  • (2) Method for Constructing Evaluation Function Using Pixel Values of Image Data as Variables
  • the method (1) is applied to projection data correction processing, which is preprocessing of image reconstruction processing.
  • projection data correction processing to which the successive approximation method (1) is applied is referred to as “sequential approximation projection data correction processing”.
  • Non-Patent Document 1 proposes a successive approximation projection data correction process using a weighted square error function with penalties expressed by the following equation as an evaluation function.
  • 1,..., I,..., I are serial numbers that uniquely identify combinations of detector channels, detector columns, and views.
  • p i and y i are provisional projection data and actual projection data specified by the serial number i, respectively.
  • p ⁇ p 1 ,..., p i ,..., p I ⁇ is a vector of provisional projection data.
  • ⁇ i is a set of detection elements adjacent to the detection element specified by the serial number i.
  • ⁇ ik is a constant indicating the correlation between the detection element specified by the serial number i and the detection element specified by the serial number k.
  • ⁇ ik is determined empirically.
  • ⁇ (p i ⁇ p k ) is a potential function with the contrast between the provisional projection data p i and the provisional projection data p k as a variable.
  • a quadratic function is used as the potential function.
  • d i is a weighting factor weighted to the difference value between the actual projection data y i and the forward projection data.
  • d i since is a value reflecting the detector output of the detection element identified by serial number i, and later, the d i is referred to as "detector output weight".
  • the detector output is a value that depends on the number of detected photons.
  • ⁇ is an arbitrary constant.
  • Non-Patent Document 2 proposes a successive approximation reconstruction process using a weighted square error function with penalties expressed by the following equation as an evaluation function.
  • Equation (2) 1,..., J,..., J are serial numbers that uniquely identify the pixels (two-dimensional coordinates) of the image data.
  • X j is a pixel value in the pixel specified by the serial number j.
  • X ⁇ X 1 ,... X j ,..., X J ⁇ is a vector representing image data.
  • a ij is a matrix element that associates image data with projection data. This matrix is called a “system matrix” because it represents the characteristics of the imaging system of the X-ray CT apparatus through a mathematical model.
  • X j is a set of pixels adjacent to the pixel specified by the serial number j.
  • w jk is a constant indicating the correlation between the pixel specified by the serial number i and the pixel specified by the serial number k.
  • w jk is the reciprocal of the distance between the center position of the pixel specified by the serial number j and the center position of the pixel specified by the serial number k.
  • Non-Patent Document 2 proposes a quadratic function as the simplest potential function.
  • Non-Patent Document 3 proposes a method of changing the shape of the potential function according to an arbitrary parameter.
  • the provisional projection data p i in the update process is expressed as the actual projection data y i and the detector output weight d i according to the first term of the equation (1). While being restricted by the above, the penalty term of the second term is smoothed according to the contrast with the provisional projection data p k that is close. Therefore, it can be said that the arbitrary constant ⁇ in the equation (1) is a parameter for adjusting the degree of smoothing of the provisional projection data p i .
  • Equation (2) can be said to be a parameter for adjusting the degree of smoothing of the image data Xj .
  • Non-Patent Document 2 the arbitrary constant ⁇ is a constant value throughout the successive approximation reconstruction process.
  • Non-Patent Document 4 proposes a method for changing ⁇ . If ⁇ remains constant, the CT image has non-uniform spatial resolution and noise reduction characteristics depending on the pixel position. Therefore, in Non-Patent Document 4, by changing ⁇ in each pixel, in addition to the contrast of adjacent pixels, the nonuniformity of spatial resolution and noise reduction characteristics is alleviated.
  • the observer of the CT image has a quantum noise reduction effect and an artifact reduction effect (hereinafter both If they are not distinguished, they are referred to as “noise reduction effect”).
  • the size of the detector output weight d i is responsible for the difference in noise reduction effect of the CT images between sites.
  • the successive approximation method based on the equations (1) and (2) is applied from the chest to the abdomen. Lung small attenuation of X-rays to (organ that is a hollow) is present, the value is large element i of the detector output weight d i is compared with the number present chest, the value of the detector output weight d i In the abdomen where the element i having a large value is small, the restrictions on the detector output weights d i in the expressions (1) and (2) are weak. That is, if the same image processing is performed on the chest and abdomen, the noise reduction effect is different. Eventually, a CT image having a variation in noise reduction effect for each part, which is undesirable for an observer, is generated.
  • Non-Patent Document 4 when ⁇ is changed in each pixel, it is possible to suppress variation in the noise reduction effect in the part. However, will be large and small equalization of the detector output weight d i, the noise reduction effect of the CT image is insufficient, i.e. the reduction of quantum noise and artifacts becomes insufficient. Furthermore, the method according to Non-Patent Document 4 has a problem that it takes a huge amount of time to determine ⁇ in each pixel.
  • the present invention has been made in view of the above-described problems, and the object of the present invention is included in projection data without greatly increasing the amount of calculation for projection data over a plurality of parts. It is an object of the present invention to provide a medical image processing apparatus or the like that can create a medical image in which uniform quantum noise and artifact reduction effects are achieved in all parts.
  • the first invention provides a correction process for projection data and / or by a successive approximation method including a detector output weight, which is a weight according to the detector output, and a penalty term in the evaluation function.
  • a medical image processing apparatus for performing reconstruction processing of a reconstructed image, wherein information that uniquely identifies a combination of a channel of the detector, a column of the detector, and a view that is an acquisition unit of the projection data is a collective element
  • a true subset determining unit for determining one or a plurality of true subsets based on imaging conditions and reconstruction conditions from the entire set, and the set elements included in the true subset for each true subset
  • a penalty term weight calculation unit for calculating a penalty term weight that is a weight related to the penalty term based on the detector output weight corresponding to the penetrating term, and the successive approximation method using the penalty term weight for each true subset Run
  • a medical image processing apparatus comprising a.
  • the second invention performs the projection data correction process and / or the reconstructed image reconstruction process by the successive approximation method including the detector output weight, which is a weight according to the detector output, and the penalty term in the evaluation function.
  • a medical image processing method to be performed wherein an imaging condition is obtained from an entire set having information that uniquely identifies a combination of a view of a channel of the detector, a row of the detector, and a view that is an acquisition unit of the projection data.
  • a medical image processing method comprising: calculating a penalty term weight that is a weight related to the penalty term; and executing the successive approximation method using the penalty term weight for each true subset.
  • a medical image in which uniform quantum noise and artifact reduction effects are achieved in all parts included in the projection data without significantly increasing the amount of calculation for the projection data over a plurality of parts. Can be created.
  • the X-ray CT apparatus 1 performs processing of data obtained from the scanner 2 on which the X-ray tube 11 and the detector 12 are mounted, the bed 4 on which the subject 10 is placed, and the detector 12.
  • An arithmetic device 5 an input device 6 such as a mouse, a trackball, a keyboard, and a touch panel, and a display device 7 for displaying a reconstructed image (CT image) and the like are included.
  • the operator inputs shooting conditions and reconstruction conditions via the input device 6.
  • the imaging conditions include, for example, a bed feeding speed, tube current, tube voltage, an imaging range (slice position range), and the number of imaging views per revolution.
  • the reconstruction conditions are, for example, a region of interest, a reconstructed image size (reconstructed image size), a reconstruction filter function, and 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 control device 15 inputs imaging conditions and reconstruction conditions from the input device 6 in the operation unit 3, and sends control signals necessary for imaging to the collimator control device 21, the X-ray control device 16, the scanner control device 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 couch controller 19 controls the couch 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 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 actual projection data can also be said to be a discrete X-ray tube position in the rotation direction (an opposing detector position). ).
  • the unit of acquisition of actual projection data at each X-ray tube position is a “view”.
  • 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 (input unit), a display device 7 (display unit), 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 superimposes a reconstruction filter on the actual projection data of each view to generate filter-corrected projection data, and weights the view-corrected weight on the filter-corrected projection data to perform backprojection processing.
  • a tomographic image is imaged non-destructively 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. Further, the reconstruction calculation device 31 displays the 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 it 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).
  • the imaging mode 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.
  • a photographing mode in which photographing is performed with the bed 4 fixed and the bed 4 is moved to the next photographing position is called step-and-shoot scanning. Since the axial scan is considered as a step-and-shoot scan in which the bed 4 is moved to the photographing position only once, both are hereinafter collectively referred to as a step-and-shoot scan.
  • the imaging mode in which the bed 4 continuously moves and the X-ray tube 11 circulates around the subject 10 in a spiral trajectory is called a spiral scan.
  • the bed control device 20 keeps the bed 4 stationary while taking a picture.
  • the bed control device 20 translates the bed 4 in the direction of the body axis during shooting according to the speed of the bed feeding as the shooting 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).
  • FIG. 3 with reference to FIG. 4, the formula (1) will be described detector output weights d i in (2).
  • FIG. 3 a cross-sectional image 41 (CT image) imaged by the X-ray CT apparatus 1 by a reconstruction process, a vertical view (a) and a horizontal view of the cross-sectional image 41 of the shoulder portion of the phantom simulating a human body
  • CT image CT image
  • a direction view (b) is shown schematically.
  • the dimension of the column of the detector 12 is not considered here, but is limited to the two dimensions of the view and the channel 42.
  • the horizontal view (b) there are many channels 42 in which the X-ray transmission length in the human body is long and the detector output shows a small value compared to the vertical view (a). That is, in the horizontal view (b), the detector output includes a large amount of noise. Note that this causes streak-like artifacts mainly in the horizontal direction in the CT image reconstructed by the image reconstruction method of
  • the detector output weights d i in Equations (1) and (2) are values corresponding to the detector outputs, the detector output weights d i corresponding to the vertical view (a) are more in the horizontal direction.
  • the detector output weights d i corresponding to the view (b) are smaller.
  • the constraint by the first term in equations (1) and (2) is relatively weaker, and the second term (penalty term) ) Will be strongly affected by the smoothing effect. Therefore, in the successive approximation method using the evaluation functions such as Expressions (1) and (2), noise included in the detector output can be averaged, and streak-like artifacts generated in the CT image can be effectively reduced. Further, with the same principle, the successive approximation method can effectively reduce the quantum noise generated in the CT image.
  • the magnitude of the detector output weight d i causes a difference in the noise reduction effect of the CT image between the parts.
  • FIG. 4 for the chest and abdomen of a phantom simulating a human body, a coronal plane image 43 (CT image) imaged by the X-ray CT apparatus 1 by reconstruction processing and a view corresponding to the chest of the coronal plane image 43 are displayed.
  • CT image coronal plane image
  • the range (a) and the view range (b) corresponding to the abdomen are schematically shown.
  • the channel dimensions of the detector 12 are not considered here, but are limited to the two dimensions of the view and column 44. 45 indicates the scanning direction, that is, the direction opposite to the direction in which the bed 4 is moved.
  • the average detector output value in the chest is larger than the average detector output value in the abdomen. This is because the amount of attenuation of X-rays is relatively small compared to the abdomen and the like because the chest includes the lung (hollow part).
  • the variation in noise reduction effect due to this part is made uniform by medical image processing described later.
  • the medical image processing apparatus of the present invention may be the arithmetic device 5 included in the X-ray CT apparatus 1 or a general-purpose computer not included in the X-ray CT apparatus 1. Furthermore, the X-ray CT apparatus 1 and the medical image processing apparatus may not be connected via a network. In order to avoid confusion, the arithmetic device 5 will be described below as the medical image processing device of the present invention.
  • the input unit, the display unit, and the storage unit included in the medical image processing apparatus of the present invention may be the input device 6, the display device 7, and the storage device 8 included in the X-ray CT apparatus 1, or the X-ray CT.
  • a device provided in a general-purpose computer not included in the device 1 or an external device may be used.
  • the input device 6, the display device 7, and the storage device 8 will be described below as an input unit, a display unit, and a storage unit included in the medical image processing apparatus of the present invention.
  • the arithmetic device 5 performs projection data correction processing by a successive approximation method including a detector output weight (weight according to the output of the detector 12) and a penalty term in the evaluation function.
  • the operator inputs imaging conditions and reconstruction conditions, and a desired calculation time and desired noise reduction rate of the successive approximation process via the input device 6 (step 1).
  • the arithmetic unit 5 may automatically set the values stored in the storage unit 8 together with the projection data. For example, the arithmetic device 5 displays the value of the photographing condition stored in the storage device 8 on the display device 7, and the operator only needs to perform the confirmation work.
  • the arithmetic device 5 displays default values, options, and the like on the display device 7 so as to support the operator's input work.
  • the arithmetic device 5 may display options such as “high speed”, “medium speed”, and “low speed” on the display device 7. In this case, the arithmetic device 5 stores the calculation time for each option in the storage device 8. Then, the arithmetic device 5 sets a desired calculation time according to the option selected by the operator.
  • the arithmetic device 5 may display options such as “high image quality”, “medium image quality”, and “low image quality” on the display device 7. In this case, the arithmetic device 5 stores the noise reduction rate for each option in the storage device 8. Then, the arithmetic device 5 sets a desired noise reduction rate in accordance with the option selected by the operator.
  • the computing device 5 calculates one or more true subsets ⁇ S 1 ,..., S m ,..., S M ⁇ based on the imaging conditions and reconstruction conditions input in step 1. Determine (step 2).
  • step-and-shoot scanning In the case of the step-and-shoot scan, it is assumed that the same part is imaged while the position of the bed 4 is fixed. Then, the arithmetic device 5 determines a set element of each true subset so that projection data collected while the position of the bed 4 is fixed belongs to the same true subset.
  • the number of divisions of the entire set ⁇ that is, the number M of the true subset S m is equal to the number of times the position of the bed 4 is moved.
  • number of views V 1 to be photographed at each position of the bed 4 ⁇ ⁇ ⁇ , V m, ⁇ ⁇ ⁇ , and V M.
  • the number of views overall shooting V V 1 + ⁇ + V m + ⁇ + V M.
  • the number of set elements of the subset S m is the view number V m to be photographed at each position of the bed 4, the number of channels C of the detector 12, the number obtained by multiplying the number of columns R of the detector 12 V m ⁇ C ⁇ R.
  • the arithmetic unit 5 determines the number of set elements included in the true subset S m based on the number of views determined as the calculation unit of the back projection process in the successive approximation method, that is, the number of back projection views. To do. This is because the arithmetic unit 5 needs to determine each true subset S m so as to have projection data necessary for creating one cross-sectional image as a set element.
  • the number of backprojection views corresponds to the phase width tw of the backprojection process (hereinafter referred to as “backprojection phase width tw”) as an imaging condition.
  • the value is obtained by multiplying the set number of shooting views V d per round.
  • the value indicating the phase is 2, etc.
  • the backprojection phase width tw is an arbitrary parameter, and in Non-Patent Document 5, it is determined empirically according to the bed feeding speed and the reconstructed image size.
  • the computing device 5 stores the back projection phase width tw in the storage device 8 in advance for each bed feeding speed and reconstructed image size. Then, the arithmetic device 5 acquires a single back projection phase width tw corresponding to the bed feeding speed set as the imaging condition and the reconstructed image size set as the reconstruction condition from the storage device 8. Next, the arithmetic unit 5 sets a value obtained by multiplying the number of photographic views V d per revolution set as the photographic condition by the back projection phase width tw acquired from the storage device 8 as the number of back projection views.
  • the true subset may be determined after replacing the number of backprojection views with 1/2 round (minimum unit of backprojection). Note that the value of the reconstruction condition is applied to the back projection process regardless of this replacement.
  • the number of divisions of the entire set ⁇ that is, the number M of the true subset S m is a quotient obtained by dividing the number of views V of the entire imaging by the number of backprojection views tw ⁇ V d .
  • Last View entire photographing V is, if not divisible by back projection view number tw ⁇ V d, remainder view, included at the end of the subset S M.
  • the number of set elements in each true subset S m is the same, and the number tw ⁇ V d ⁇ C multiplied by the backprojection phase width tw, the number of captured views V d per round, the number of channels C, and the number of columns R ⁇ R.
  • the calculation unit 5 for each subset S m determined in step 2, based on the detector output weights d i corresponding to the set elements as fall within the true subset S m, according to the penalty term Weights (penal term weights) ⁇ 1 ,..., ⁇ m ,..., ⁇ M ⁇ are calculated (step 3). Details of the penalty term weight calculation processing will be described later with reference to FIGS.
  • the arithmetic unit 5 executes a successive approximation method for i ⁇ S m using the penalty term weight ⁇ m for each true subset S m (step 4). That is, the arithmetic unit 5 executes the successive approximation projection data correction process while changing the penalty term weight ⁇ m for each part.
  • the following expression is an example of an evaluation function.
  • the technical idea of the present invention is to change the penalty term weight in the evaluation function for each part, if the evaluation function includes the detector output weight and the penalty term, the numerical analysis used for the optimization It can be applied regardless of the law.
  • Equation (3) Deriving the update formula in the successive approximation method from Equation (3) can be realized by a known method.
  • the update equation can be derived by using the GauSS-Seidel method.
  • Arithmetic unit 5 derives the update equation in the successive approximation from Equation (3), a penalty term weight beta m of the calculated detector output weights d i and the true part each set S m, the derived updated formula Substituting and correcting the projection data by the successive approximation method.
  • the operator can arbitrarily select one of the two in consideration of the desired calculation time and the accuracy of the penalty term weight estimation.
  • the arithmetic unit 5 calculates a penalty term weight value (hereinafter referred to as ⁇ penalty '') that achieves a noise reduction effect equivalent to the noise reduction rate with respect to the reference phantom for each noise reduction rate.
  • Term term reference value a penalty term weight value (hereinafter referred to as ⁇ penalty '') that achieves a noise reduction effect equivalent to the noise reduction rate with respect to the reference phantom for each noise reduction rate.
  • Term term reference value is stored in the storage device 8 in advance.
  • the arithmetic unit 5 registers penalties weight reference values ⁇ * 10 , ⁇ * 20 , ⁇ * 30 ,... For noise reduction rates of 10%, 20%, 30%,.
  • the stored noise reduction rate table is stored in the storage device 8.
  • the penalties weight reference values ⁇ * 10 , ⁇ * 20 , ⁇ * 30 , ... are taken as a circular phantom or a phantom simulating a human body as a reference phantom in advance. The value obtained by executing the configuration process.
  • the arithmetic device 5 acquires a penalty term weight reference value corresponding to the input noise reduction rate from the storage device 8.
  • the calculation unit 5, for each subset S m calculates a representative value based on the detector output weights d i corresponding to the set elements contained within each subset S m.
  • the calculation unit 5 calculates a representative value based on the detector output weights d i corresponding to the set elements contained within each subset S m.
  • the calculation device 5 each subset S m, the average and median of the detector output weight d i corresponding to the set elements as fall within the true subset S m, as well as the detector output weights d i It is desirable that any one of the three values of the class that divides the whole by a predetermined ratio in the histogram in which is a class is a representative value for each true subset S m .
  • the arithmetic device 5 selects a table value corresponding to the desired noise reduction rate input by the operator from the noise reduction rate table stored in the storage device 8 (step 11).
  • the computing device 5 creates a subset T m of the true subset S m from the set elements included in the true subset S m according to the desired calculation time input by the operator (step 12).
  • the arithmetic device 5 creates the lower set T m by thinning out the set elements included in S m based on a predetermined thinning rule.
  • the view is thinned out repeatedly. This is because when the view is skipped, the influence on the image quality of the finally created cross-sectional image is small.
  • the number of captured views Vd per rotation is 360.
  • the number of set elements of the sub-set T m is 4/5 of the number of set elements of the true subset S m. It becomes.
  • the thinning amount can be adjusted according to the desired calculation time input by the operator. In this way, the calculation time can be adjusted in accordance with the intention of the operator by creating the subset T m of the true subset S m and executing the processing described later on the subset T m .
  • the arithmetic unit 5 constructs a histogram using the detector output weights d i corresponding to the set elements of the lower set T m created in step 12 as a class (step 13).
  • calculation method of the detector output weight d i may use any known technique.
  • the arithmetic device 5 may perform calculation by multiplying the air data after performing inverse logarithmic transformation on the projection data.
  • FIG. 7 shows an example of a histogram.
  • the horizontal axis (class) is the detector output weight d i
  • the vertical axis (frequency) is the number of set elements corresponding to the detector output weight d i equal to the value of each class.
  • Histograms have different distributions depending on the size and part of the subject and can be said to be an index that well expresses their characteristics.
  • the maximum peak exists where the class is a value close to zero. Further, the histogram of the shoulder portion 51, overall, a lower value of the detector output weights d i. This is because the shoulder 51 includes many bones with a large amount of X-ray attenuation.
  • the maximum peak exists where the class is close to 1. Comparing the histograms of the chest 52 and the abdomen 53, the frequency of the class 1 to 3 is the histogram of the chest 52> the histogram of the abdomen 53, and the frequency of the class 5 or more is the histogram of the chest 52 ⁇ the histogram of the abdomen 53. ing. This is because the chest 52 includes a lung with a small amount of X-ray attenuation.
  • the arithmetic unit 5 calculates the area of each histogram shown in FIG. 7, and specifies a class that divides the area of the histogram into a predetermined ratio as a representative value (step 14).
  • the representative value may be a median value or an average value of the detector output weights d i included in i ⁇ S m . The operator can arbitrarily select these.
  • the computing device 5 calculates the penalty term weight ⁇ m of the true subset S m by multiplying the table value selected in Step 11 by the class (representative value) specified in Step 14 ( Step 15).
  • the penalty term weight ⁇ m of the true subset S m can be calculated without significantly increasing the amount of calculation.
  • the second penalty term weight calculation process is based on the following two findings. (1) It is empirically clear that when the detector output weight is a constant and the successive approximation method is executed with an arbitrary penalty term weight, almost the same noise reduction rate can be obtained regardless of the subject and the part. (2) There is a high correlation between the correction amount of the projection value (pixel value) by the successive approximation method and the noise reduction rate.
  • the arithmetic unit 5 calculates the true subset S for the correction amount calculation function for calculating the correction amount of the projection data when the successive approximation method is executed. Substitute the detector output weights d i and projection data corresponding to the set elements included in m .
  • the process for determining the predetermined correction amount reference value is as follows. Similar to the first penalty term weight calculation process, the arithmetic device 5 stores the noise reduction rate table in the storage device 8, and from the storage device 8, the penalty term weight reference value corresponding to the input noise reduction rate is stored. get. Then, the calculation device 5 calculates the correction amount of the projection data when the successive approximation method is executed using the penalty term weight reference value acquired from the storage device 8 with the detector output weight as a constant, and the predetermined correction Use quantity reference value.
  • FIG. 8 shows processing for calculating the penalty term weight ⁇ m for the true subset S m , but the same applies to other true subsets.
  • the arithmetic device 5 selects a table value corresponding to the desired noise reduction rate input by the operator from the noise reduction rate table stored in the storage device 8 (step 21). Note that the stored table values are different between the noise reduction rate table used for the first penalty term weight calculation process and the noise reduction rate table used for the second penalty term weight calculation process due to the difference in processing contents.
  • the computing device 5 creates a subset T m of the true subset S m from the set elements included in the true subset S m according to the desired calculation time input by the operator (step 22).
  • the process of creating the subset T m is the same as in Step 12.
  • the calculation unit 5 the selected table value at step 21, and, using the projection data corresponding to the set elements of subset T m created in step 22, the detector output weight d i any A projection value after application of the successive approximation projection data correction process is estimated with a constant.
  • the approximate value is substituted. Assuming that the projected value of interest is calculated from only the adjacent channels, columns, and views, the approximate value can be easily calculated. Further, the approximate value may be calculated using a known method after replacing the inverse matrix necessary for calculating the projection value after the iteration with an approximate inverse matrix.
  • the arithmetic unit 5 calculates an error between the approximate value and the projection value before application of the successive approximate projection data correction process, and sets the sum of the errors as a correction amount reference value (step 23).
  • a known index such as an absolute error or a square error can be used for the error calculation process.
  • the arithmetic device 5 substitutes the detector output weight di and the projection data corresponding to the set element of the sub-set T m created in Step 23 into the correction amount calculation function using the penalty term weight ⁇ m as a variable. . Then, the arithmetic unit 5 determines the penalty term weight so that the value of the correction amount calculation function is equal to the correction amount reference value calculated in step 23 (step 24).
  • a known numerical analysis method for example, a bisection method
  • the penalty term weight ⁇ m of the true subset S m can be calculated with high accuracy.
  • the penalty term weight ⁇ m used for the successive approximation projection data correction process is set to a constant value in the true subset S m .
  • the noise reduction effect of the CT image can be achieved.
  • a second embodiment will be described.
  • the arithmetic device 5 performs reconstruction processing of a reconstructed image by a successive approximation method that includes the detector output weight and the penalty term in the evaluation function.
  • equation (4) as in equation (3), an update equation in the successive approximation method can be derived from equation (4).
  • Arithmetic device 5 a penalty term weight beta m of the detector for each output weight d i and subset S m is calculated as in the first embodiment, it substitutes the derived updated formula, successive approximation Thus, the reconstructed image is reconstructed.
  • the penalty term weight ⁇ m used for the successive approximation reconstruction process is set to a constant value in the true subset S m as in the first embodiment.
  • the noise reduction effect of the CT image can be achieved as in the conventional case.
  • ⁇ Third embodiment> A third embodiment will be described.
  • the arithmetic device 5 performs the projection data correction process and the reconstructed image reconstruction process by the successive approximation method including the detector output weight and the penalty term in the evaluation function.
  • the third embodiment is a combination of the first embodiment and the second embodiment. That is, the arithmetic device 5, as in the first embodiment, the formula (3) derives the update equation in the successive approximation from the detector output weights d i and calculated in the same manner as in the first embodiment The penalty term weight ⁇ m for each true subset S m is substituted into the derived update formula, and the projection data is corrected by the successive approximation method.
  • the arithmetic device 5 as in the second embodiment, the formula (4) to derive the update equation in the successive approximation from the detector output weights d i and calculated in the same manner as in the first embodiment
  • the penalty term weight ⁇ m for each true subset S m is substituted into the derived update formula, and the reconstructed image is reconstructed by the successive approximation method.
  • the penalty term weight ⁇ m used for the successive approximation projection data correction process and the successive approximation reconstruction process is set to a constant value within the true subset S m , so that In the approximate projection data correction process and the successive approximation reconstruction process, the noise reduction effect of the CT image can be achieved as in the conventional case.
  • the true subset S m is calculated between the regions by calculating the penalty term weight ⁇ m for each true subset S m based on the detector output weights d i reflecting the subject information and the projection values. Variations in the noise reduction effect can be suppressed.

Landscapes

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

Abstract

　複数の部位に亘る投影データに対して、投影データに含まれる全ての部位において一様なノイズ低減効果が達成された医用画像を作成するために、演算装置５は、撮影条件及び再構成条件に基づいて、1又は複数の真部分集合を決定する(ステップ2)。次に、演算装置5は、真部分集合毎に、真部分集合内に含まれる集合要素に対応する検出器出力重みに基づいて、罰則項重みを算出する(ステップ3)。次に、演算装置5は、真部分集合毎の罰則項重みを用いて逐次近似法を実行する(ステップ4)。

Description

医用画像処理装置、医用画像処理方法
　本発明は、検出器の出力に応じた重み及び罰則項を評価関数に含む逐次近似法によって、投影データの補正処理及び/又は再構成画像の再構成処理を行う医用画像処理装置等に関するものである。
　X線CT装置は、チャネル方向及び列方向に検出素子が配列されるX線検出器を有する。そして、X線CT装置による撮影では、互いに対向するX線管とX線検出器を被検体の周囲に回転させることにより、投影データが周回方向の離散的なX線管の位置(対向するX線検出器の位置とも言える。)において収集される。以下では、各々のX線管位置における投影データの取得単位を「ビュー」と呼ぶ。投影データを収集した後、演算装置が、各ビュー及び各列について、投影データのチャネル方向に再構成フィルタを重畳してフィルタ補正投影データを得て、フィルタ補正投影データに対して、ビュー方向に重みを加重しながら逆投影することによって、被検体内部のX線減弱係数の分布図として非破壊的に断層像を画像化する。以下では、投影データから断層像を画像化する方法を「画像再構成法」と呼び、画像再構成法によって得られた画像をCT画像と呼ぶ。
　画像再構成法は、解析法と逐次近似法に大別される。解析法は、投影切断面定理に基づき解析的に問題を解く方法である。逐次近似法は、投影データの取得に至る観測系を数学的にモデル化し、数学モデルに基づいて反復法により最良の画像を推定する方法である。
　両者を比較すると、解析法の利点としては、投影データから直接的に再構成画像が得られるため、計算量が圧倒的に少ないことが挙げられる。一方、逐次近似法の利点としては、投影データの取得に至る物理的過程及び投影データに含まれる統計的な揺らぎをそれぞれ数学モデルや統計モデルとして考慮できるため、解析法において発生するアーチファクト(コーンビームアーチファクトなど)や画像上の量子ノイズを低減できることが挙げられる。
　従来、マルチスライスCTにおける画像再構成法としては、計算量が少ないことを採用理由として、解析法であるＦｅｌdkamp法、またはＦｅｌdkamp法を改良した手法が主に用いられている。しかし、近年のコンピュータの高性能化に伴い、逐次近似法も実用化が検討され始めている。特に、逐次近似法によって画質の向上を図る技術が提案されている。
　ここでは、推定される変数の違いによって、逐次近似法を次の2種類に大別する。
(1)投影データの投影値を変数として評価関数を構築する手法
(2)画像データの画素値を変数として評価関数を構築する手法
　(1)の手法は、画像再構成処理の前処理である投影データの補正処理などに適用される。以下では、(1)の逐次近似法が適用される投影データの補正処理を「逐次近似投影データ補正処理」と呼ぶ。
　(2)の手法は、画像再構成処理に適用される。以下では、(2)の逐次近似法が適用される画像再構成処理を「逐次近似再構成処理」と呼ぶ。
　逐次近似投影データ補正処理では、投影データの評価指標を事前に設定しておき、評価指標を数値化した評価値が最大値もしくは最小値をとるように投影データを逐次更新する。評価指標は、更新の対象である暫定投影データと実投影データとの矛盾や確率的な尤もらしさなどが用いられる。一例として、非特許文献1には、次式によって表される罰則付き加重二乗誤差関数を評価関数として用いる逐次近似投影データ補正処理が提案されている。
Figure JPOXMLDOC01-appb-M000001
 　(1)式において、1，・・・i，・・・，Iは、検出器のチャネル、検出器の列、及びビューの組み合わせを一意に識別する通し番号である。pi及びyiは、それぞれ通し番号iによって特定される暫定投影データ及び実投影データである。p＝｛p1，・・・pi，・・・，pI｝は、暫定投影データのベクトルである。
κiは、通し番号iによって特定される検出素子と近接する検出素子の集合である。
　νikは、通し番号iによって特定される検出素子と、通し番号kによって特定される検出素子との相関性を示す定数である。非特許文献1では、νikを経験的に決定している。
　Ψ(pi－pk)は、暫定投影データpiと暫定投影データpkとのコントラストを変数とするポテンシャル関数である。非特許文献1では、ポテンシャル関数として、2次関数を使用している。
　diは、実投影データyiと順投影データの差分値に加重される重み係数である。CT画像の画像処理では、diは通し番号iによって特定される検出素子の検出器出力を反映した値である為、以降では、diを「検出器出力重み」と呼ぶ。尚、検出器出力は、検出フォトン数に依存する値である。
　βは、任意定数である。
　逐次近似再構成処理では、画像データの評価指標を事前に設定しておき、評価指標を数値化した評価値が最大値もしくは最小値をとるように画像データを逐次更新する。評価指標は、更新の過程において暫定画像データが投影データに変換されたデータである順投影データと、実投影データとの矛盾や確率的な尤もらしさなどが用いられる。一例として、非特許文献2には、次式によって表される罰則付き加重二乗誤差関数を評価関数として用いる逐次近似再構成処理が提案されている。
Figure JPOXMLDOC01-appb-M000002
 　(2)式において、1，・・・j，・・・，Ｊは、画像データの画素(2次元座標)を一意に識別する通し番号である。Xjは、通し番号jによって特定される画素における画素値である。X＝｛X1，・・・Xj，・・・，XＪ｝は、画像データを表すベクトルである。
　aijは、画像データと投影データを対応付ける行列の要素である。この行列は、数学モデルを介してX線CT装置の撮影系の特性を表すことから、「システムマトリクス」と呼ばれている。
　Xjは、通し番号jによって特定される画素と近接する画素の集合である。wjkは、通し番号iによって特定される画素と、通し番号kによって特定される画素との相関性を示す定数である。一般的には、wjkは、通し番号jによって特定される画素の中心位置と、通し番号kによって特定される画素の中心位置との距離の逆数である。
　Ψ(Xj－Xk)は、通し番号jによって特定される画素における画素値Xjと、通し番号kによって特定される画素における画素値Xkとのコントラストを変数とするポテンシャル関数である。非特許文献2では、最も簡単なポテンシャル関数として、2次関数が提案されている。また、非特許文献3には、ポテンシャル関数の形状を任意のパラメータによって変化させる手法が提案されている。
　その他の定数は、(1)式と同様である。
　(1)式から種々の数値解析法により導出される更新式において、更新の過程における暫定投影データpiは、(1)式の第1項により実投影データyi及び検出器出力重みdiによる制約を受けながら、第2項の罰則項により近接する暫定投影データpkとのコントラストに応じて平滑化される。従って、(1)式の任意定数βは、暫定投影データpiの平滑化の度合いを調整する為のパラメータであると言える。
　同様に、(2)式の任意定数βは、画像データXjの平滑化の度合いを調整する為のパラメータと言える。
　非特許文献2及び非特許文献3では、任意定数βは、逐次近似再構成処理全体を通して一定値である。
　一方、非特許文献4では、βを変更する手法が提案されている。βが一定のままでは、CT画像は、画素の位置によって空間分解能及びノイズ低減特性が不均一になる。そこで、非特許文献4では、各画素においてβを変更することによって、近接する画素のコントラストに加え、空間分解能及びノイズ低減特性の不均一性を緩和している。
T. Li et. al., "Nonlinear SinogramSmoothing for Low-DoSe X-Ray CT,"IEEE. TranS. Nucl. Sci., Vol.51, No.5,pp.2505-2513, 2004 K. Saueret. al., "A Local Update Strategy for Iterative ReconStruCTion form ProjeCTionS,"IEEE.TranS. Signal. Proc., Vol.41, No.2, pp.534-548, 1993 J. B. Thibault et. al., "Athree-dimenSional StatiStical approach to improved image quality for multiSlicehelical CT,"Med. PhyS., Vol.34, No.11, pp.4526-4544, 2007 H. R. Shi et. al., "Quadratic Regularization DeSign for 2-D CT,"IEEE. TranS. Med. Imag., Vol.28, No.5, pp645-656, 2009 T. Goto et. al., "Weihgted-Feldkamp algorithm with SeleCTive narroweSt cone angle data for conebeam CT,"Proc. Of. The Eigth Meeting on Fully Three-dimenSional Image ReconStruCTion in Radiology and Nuclear Medocine, pp.189-192, 2005 J. Wang et. al., "Penalized weighted leaSt-SquareS approach to Sonogram noiSe reduCTion and image reconStruCTion for low-doSe X-ray computed tomography,"IEEE. TranS. Med. imag.,Vol.25, No.10, pp.1272-1283, 2006
　ところで、逐次近似投影データ補正処理や逐次近似再構成処理を複数の部位(胸部及び腹部等)に適用する際、CT画像の観察者にとっては、量子ノイズ低減効果及びアーチファクト低減効果(以下、両者を区別しない場合には「ノイズ低減効果」と呼ぶ。)が、部位に因らず同等であることが望ましい。
　しかしながら、複数部位に逐次近似投影データ補正処理や逐次近似再構成処理を適用する場合、検出器出力重みdiの大小は、部位間におけるCT画像のノイズ低減効果の違いの原因となる。例えば、胸部から腹部にかけて、式(1)や式(2)に基づく逐次近似法を適用した場合を考える。肺(中空になっている臓器)が存在する為にX線の減衰が少なく、検出器出力重みdiの値が大きい要素iが数多く存在する胸部と比較すると、検出器出力重みdiの値が大きい要素iが少ない腹部では、(1)式及び(2)式における検出器出力重みdiの制約が弱い。すなわち、胸部と腹部に対して同じ画像処理を施すと、ノイズ低減効果が異なることになる。ひいては、観察者にとって好ましくない、部位ごとにノイズ低減効果のばらつきがあるCT画像が生成されてしまう。
　また、非特許文献4のように、各画素においてβを変更する場合、部位におけるノイズ低減効果のばらつきを抑制することはできる。しかしながら、検出器出力重みdiの大小が均等化されてしまい、CT画像のノイズ低減効果が不十分、すなわち量子ノイズ及びアーチファクトの低減が不十分になってしまう。更に、非特許文献4による手法では、各画素においてβを決定する処理に膨大な時間を要するという問題もある。
　本発明は、前述した問題点に鑑みてなされたものであり、その目的とすることは、複数の部位に亘る投影データに対して、計算量を大幅に増加させることなく、投影データに含まれる全ての部位において一様な量子ノイズ及びアーチファクトの低減効果が達成された医用画像を作成することができる医用画像処理装置等を提供することである。
　前述した目的を達成するために第1の発明は、検出器の出力に応じた重みである検出器出力重み、及び罰則項を評価関数に含む逐次近似法によって、投影データの補正処理及び/又は再構成画像の再構成処理を行う医用画像処理装置であって、前記検出器のチャネル、前記検出器の列、及び前記投影データの取得単位であるビューの組み合わせを一意に識別する情報を集合要素とする全体集合から、撮影条件及び再構成条件に基づいて1又は複数の真部分集合を決定する真部分集合決定部と、前記真部分集合毎に、前記真部分集合内に含まれる前記集合要素に対応する前記検出器出力重みに基づいて、前記罰則項に係る重みである罰則項重みを算出する罰則項重み算出部と、前記真部分集合毎の前記罰則項重みを用いて前記逐次近似法を実行する逐次近似法実行部と、を備える医用画像処理装置である。
　第2の発明は、検出器の出力に応じた重みである検出器出力重み、及び罰則項を評価関数に含む逐次近似法によって、投影データの補正処理及び/又は再構成画像の再構成処理を行う医用画像処理方法であって、前記検出器のチャネル、前記検出器の列、及び前記投影データの取得単位であるビューの組み合わせを一意に識別する情報を集合要素とする全体集合から、撮影条件及び再構成条件に基づいて1又は複数の真部分集合を決定するステップと、前記真部分集合毎に、前記真部分集合内に含まれる前記集合要素に対応する前記検出器出力重みに基づいて、前記罰則項に係る重みである罰則項重みを算出するステップと、前記真部分集合毎の前記罰則項重みを用いて前記逐次近似法を実行するステップと、を含む医用画像処理方法である。
　本発明により、複数の部位に亘る投影データに対して、計算量を大幅に増加させることなく、投影データに含まれる全ての部位において一様な量子ノイズ及びアーチファクトの低減効果が達成された医用画像を作成することができる。
X線CT装置1の全体外観図 X線CT装置1の構成図 人体を模擬したファントムの肩部の断面像、並びに断面像に対する鉛直方向のビュー及び水平方向のビューを示す模式図 人体を模擬したファントムの冠状面像、並びに冠状面像の胸部に対応するビューの範囲及び腹部に対応するビューの範囲を示す模式図 医用画像処理全体の流れを示すフローチャート 第1の罰則項重み算出処理を示すフローチャート ヒストグラムの一例を示す図 第2の罰則項重み算出処理を示すフローチャート
　以下図面に基づいて、本発明の実施形態を詳細に説明する。最初に、図1、図2を参照しながら、X線CT装置1の構成を説明する。
　図1に示すように、X線CT装置1は、X線管11や検出器12が搭載されるスキャナ2、被検体10を載置する寝台4、検出器12から得られるデータの処理を行う演算装置5、マウス、トラックボール、キーボード、タッチパネルなどの入力装置6、及び再構成画像(CT画像)などを表示する表示装置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線の照射が行われる(但し、スキャノグラム撮影は除く。
　撮影中に寝台4が固定され、X線管11が被検体10の周りを円軌道状に周回する撮影態様は、アキシャルスキャンなどと呼ばれる。特に、寝台4を固定して撮影し、寝台4を次の撮影位置に移動させることを繰り返す撮影態様は、ステップ・アンド・シュートスキャンなどと呼ばれる。アキシャルスキャンは、1回だけ寝台4を撮影位置に移動させるステップ・アンド・シュートスキャンと考えられることから、以降では、両者を纏めてステップ・アンド・シュートスキャンと表記する。
　また、寝台4が連続的に移動し、X線管11が被検体10の周りをらせん軌道状に周回する撮影態様は、らせんスキャンなどと呼ばれる。
　寝台制御装置20は、ステップ・アンド・シュートスキャンの場合、撮影している間、寝台4を静止した状態とする。また、寝台制御装置20は、らせんスキャンの場合、入力装置6を介して入力される撮影条件としての寝台送りの速さに応じて、撮影している間、寝台4を体軸方向に平行移動させる。
　X線CT装置1は、例えば、マルチスライスCTである。また、X線CT装置1のスキャン方式は、例えば、ローテート－ローテート方式(第3世代)である。
　次に、図3、図4を参照しながら、式(1)、(2)における検出器出力重みdiについて説明する。
　図3では、人体を模擬したファントムの肩部について、X線CT装置1が再構成処理によって画像化した断面像41(CT画像)と、この断面像41に対する鉛直方向のビュー(a)及び水平方向のビュー(b)を模式的に示している。説明を分かり易くする為、ここでは、検出器12の列の次元を考えず、ビュー及びチャネル42の2つの次元に限定して考える。水平方向のビュー(b)では、鉛直方向のビュー(a)と比較して人体におけるX線の透過長が長く、検出器出力が小さな値を示すチャネル42が多い。つまり、水平方向のビュー(b)では、検出器出力に大きなノイズが含まれる。尚、このことにより、解析法の画像再構成法によって再構成されたCT画像では、主に水平方向にストリーク状のアーチファクトが生ずる。
　式(1)及び(2)の検出器出力重みdiは、検出器出力に応じた値である為、鉛直方向のビュー(a)に対応する検出器出力重みdiよりも、水平方向のビュー(b)に対応する検出器出力重みdiの方が小さな値となる。つまり、水平方向のビュー(b)では、鉛直方向のビュー(a)と比較して、式(1)及び(2)の第1項による制約が相対的に弱くなり、第2項(罰則項)による平滑化効果を強く受けることになる。従って、式(1)及び(2)などの評価関数を用いる逐次近似法では、検出器出力に含まれるノイズを平均化し、CT画像において生ずるストリーク状のアーチファクトを効果的に低減することができる。また、同様の原理によって、逐次近似法では、CT画像において生ずる量子ノイズを効果的に低減することができる。
　CT画像の量子ノイズ低減効果及びストリークアーチファクト低減効果(＝ノイズ低減効果)は、部位に因らず同程度であることが望ましい。しかし、複数部位に亘り両者を適用する場合、検出器出力重みdiの大小は、部位間におけるCT画像のノイズ低減効果の違いの原因となる。
　図4では、人体を模擬したファントムの胸部及び腹部について、X線CT装置1が再構成処理によって画像化した冠状面像43(CT画像)と、この冠状面像43の胸部に対応するビューの範囲(a)及び腹部に対応するビューの範囲(b)を模式的に示している。
説明を分かり易くする為、ここでは、検出器12のチャネルの次元を考えず、ビュー及び列44の2つの次元に限定して考える。45は、スキャン方向、つまり寝台4を移動させる方向と正反対の方向を示している。
　図4のように、胸部から腹部をX線CT装置1によって画像化する場合、胸部における平均的な検出器出力の値は、腹部における平均的な検出器出力の値よりも大きい。これは、胸部は肺(中空部)を含む為、腹部等と比較して、X線の減弱量が相対的に小さいからである。
　つまり、式(1)及び(2)などの評価関数を用いて、胸部から腹部にかけて、部位によらず同様の逐次近似法を適用する場合、大きい値の検出器出力重みdiが多い胸部と比較して、小さい値の検出器出力重みdiが多い腹部では、式(1)及び(2)の制約が相対的に弱くなる。従って、式(1)及び(2)などの評価関数を用いる逐次近似法では、CT画像において生ずるストリーク状のアーチファクト及び量子ノイズの低減効果(＝ノイズ低減効果)が、部位によってばらついてしまう。
　本発明の医用画像処理装置では、後述する医用画像処理によって、この部位によるノイズ低減効果のばらつきを一様にする。
　本発明の医用画像処理装置が実行する各ステップは、実投影データが収集された後に実行される。従って、本発明の医用画像処理装置は、X線CT装置1に含まれる演算装置5でも良いし、X線CT装置1に含まれない汎用のコンピュータでも良い。更に言えば、X線CT装置1と医用画像処理装置は、ネットワークを介して接続されていなくても良い。
混乱を避ける為に、以下では、演算装置5を本発明の医用画像処理装置として説明する。
　同様に、本発明の医用画像処理装置が備える入力部、表示部、及び記憶部は、X線CT装置1に含まれる入力装置6、表示装置7、及び記憶装置8でも良いし、X線CT装置1に含まれない汎用のコンピュータが備える装置又は外部装置でも良い。混乱を避ける為に、以下では、入力装置6、表示装置7、及び記憶装置8を本発明の医用画像処理装置が備える入力部、表示部、及び記憶部として説明する。
　また、実投影データは、ファンビーム系で収集されるが、以降ではパラレルビーム系に変換したデータを前提として説明する。
　＜第1の実施の形態＞
　図5～図8を参照しながら、第1の実施の形態について説明する。第1の実施の形態では、本発明を逐次近似投影データ補正処理に適用する場合について説明する。第1の実施の形態では、演算装置5が、検出器出力重み(検出器12の出力に応じた重み)及び罰則項を評価関数に含む逐次近似法によって、投影データの補正処理を行う。
　図5に示すように、操作者は、入力装置6を介して、撮影条件及び再構成条件、並びに、逐次近似処理の所望計算時間及び所望ノイズ低減率を入力する(ステップ1)。
　撮影条件については、演算装置5は、投影データとともに記憶装置8に記憶されている値を自動的に設定すれば良い。例えば、演算装置5が、記憶装置8に記憶されている撮影条件の値を表示装置7に表示し、操作者は確認作業のみ行えば良い。
　また、再構成条件、並びに、逐次近似処理の所望計算時間及び所望ノイズ低減率については、演算装置5は、表示装置7にデフォルト値や選択肢などを表示し、操作者の入力作業を支援するようにしても良い。
　例えば、所望計算時間の入力支援として、演算装置5は、表示装置7に「高速」、「中速」、「低速」などの選択肢を表示するようにしても良い。この場合、演算装置5は、記憶装置8に選択肢ごとの計算時間を記憶しておく。そして、演算装置5は、操作者に選択された選択肢に応じて、所望計算時間を設定する。
　また、例えば、所望ノイズ低減率の入力支援として、演算装置5は、表示装置7に「高画質」、「中画質」、「低画質」などの選択肢を表示するようにしても良い。この場合、演算装置5は、記憶装置8に選択肢ごとのノイズ低減率を記憶しておく。そして、演算装置5は、操作者に選択された選択肢に応じて、所望ノイズ低減率を設定する。
　次に、演算装置5は、ステップ1において入力された撮影条件及び再構成条件に基づいて、1又は複数の真部分集合｛S1、・・・、Sm、・・・、SM｝を決定する(ステップ2)。
　全体集合Ωは、検出器12のチャネル、検出器12の列、及びビューの組み合わせを一意に識別する通し番号｛1，・・・i，・・・，I｝を集合要素とする。より詳細には、検出器12のチャネル数がC、検出器12の列数がR、撮影全体のビュー数がVの場合、Ω＝｛1、2、3、・・・、・・・、C×R×V－2、C×R×V－1、C×R×V｝となる。
　ここで、SmがΩの真部分集合であるとは、Sm⊆Ω(SmはΩの部分集合)、かつ、Sm≠Ωが成り立つことである。更に、以下では、任意の真部分集合の組(Si、Sj)について、Si∩Sj＝φが成り立つものとして説明する。
　まず、ステップ・アンド・シュートスキャンの場合について説明する。ステップ・アンド・シュートスキャンの場合、寝台4の位置が固定されている間は、同一の部位が撮影されていると仮定する。そして、演算装置5は、寝台4の位置が固定されている間に収集された投影データが、同一の真部分集合に属するように、各真部分集合の集合要素を決定する。
　ステップ・アンド・シュートスキャンの場合、全体集合Ωの分割数、すなわち真部分集合Smの数Mは、寝台4の位置を移動させる回数に等しい。
　ここで、寝台4の各位置において撮影されるビュー数V1、・・・、Vm、・・・、VMとする。この場合、撮影全体のビュー数V＝V1＋・・・＋Vm＋・・・＋VMである。また、各真部分集合Smの集合要素の数は、寝台4の各位置において撮影されるビュー数Vm、検出器12のチャネル数C、検出器12の列数Rを乗算した数Vm×C×Rとなる。
　従って、各真部分集合｛S1、・・・、Sm、・・・、SM｝は、S1＝｛1、2、・・・、V1×C×R｝、S2＝｛V1×C×R＋1、・・・、V1×C×R＋V2×C×R｝、・・・、Sm＝｛(V1＋V2＋・・・＋Vm－1)×C×R＋1、・・・、(V1＋V2＋・・・＋Vm－1＋Vm)×C×R｝、・・・、　SM＝｛(V1＋V2＋・・・＋VM－1)×C×R＋1、・・・、(V1＋V2＋・・・＋VM－1＋VM)×C×R｝となる。
　例えば、寝台4の各位置において、常に1周回(360°)分のビューが撮影され、1周回(＝360°)当たりのビュー数をV360とする。この場合、V1＝・・・＝Vm＝・・・＝VM＝V360、各真部分集合Smの集合要素の数＝V360×C×Rとなり、全ての真部分集合Smの集合要素の数が同一となる。
　次に、らせんスキャンの場合について説明する。らせんスキャンの場合、演算装置5は、逐次近似法における逆投影処理の計算単位として定められるビューの数、すなわち逆投影ビュー数に基づいて、真部分集合Smに含まれる集合要素の数を決定する。これは、演算装置5は、1枚の断面像を作成するために必要な投影データを集合要素として持つように、各真部分集合Smを決定する必要があるからである。
　例えば、非特許文献5に記載の方法によって逐次近似処理を行う場合、逆投影ビュー数は、逆投影処理の位相幅tw(以下、「逆投影位相幅tw」と呼ぶ。)に、撮影条件として設定される周回当たりの撮影ビュー数Vdを乗じた値となる。
　逆投影位相幅twは、1周回(＝360°)分の位相を示す値が1、半周回(＝180°)分の位相を示す値が1/2、2周回(＝720°)分の位相を示す値が2、などになる。
　ここで、逆投影位相幅twは任意のパラメータであり、非特許文献5では、寝台送り速度及び再構成画像サイズに応じて経験的に決定している。
　そこで、演算装置5は、寝台送り速度及び再構成画像サイズごとに、逆投影位相幅twを予め記憶装置8に記憶しておくことが望ましい。そして、演算装置5は、記憶装置8から、撮影条件として設定される寝台送り速度、及び再構成条件として設定される再構成画像サイズに対応する単一の逆投影位相幅twを取得する。次に、演算装置5は、撮影条件として設定される周回当たりの撮影ビュー数Vdと、記憶装置8から取得される逆投影位相幅twとを乗算した値を逆投影ビュー数とする。
　もし操作者が計算速度を優先する場合には、逆投影ビュー数を1/2周回(逆投影の最小単位)と置き換えた上で、真部分集合を決定しても良い。尚、逆投影処理は、この置き換えによらず、再構成条件の値を適用する。
　らせんスキャンの場合、全体集合Ωの分割数、すなわち真部分集合Smの数Mは、撮影全体のビュー数Vを、逆投影ビュー数tw×Vdによって除したときの商となる。撮影全体のビュー数Vが、逆投影ビュー数tw×Vdによって割りきれない場合、余りのビューは、最後の真部分集合SMに含める。
　各真部分集合Smの集合要素の数は、全て同一であり、逆投影位相幅tw、周回当たりの撮影ビュー数Vd、チャネル数C、列数Rを乗算した数tw×Vd×C×Rとなる。
　従って、各真部分集合｛S1、・・・、Sm、・・・、SM｝は、S1＝｛1、2、・・・、tw×Vd×C×R｝、S2＝｛tw×Vd×C×R＋1、・・・、2×tw×Vd×C×R｝、・・・、Sm＝｛(m－1)×tw×Vd×C×R＋1、・・・、m×tw×Vd×C×R｝、・・・、SM＝｛(M－1)×tw×Vd×C×R＋1、・・・、M×tw×Vd×C×R｝となる。
　図5の説明に戻る。次に、演算装置5は、ステップ2において決定された真部分集合Sm毎に、真部分集合Sm内に含まれる集合要素に対応する検出器出力重みdiに基づいて、罰則項に係る重み(罰則項重み)｛β1、・・・、βm、・・・、βM｝を算出する(ステップ3)。罰則項重み算出処理の詳細は、図6～図8を参照しながら後述する。
　次に、演算装置5は、i∈Smについて、真部分集合Sm毎の罰則項重みβmを用いて逐次近似法を実行する(ステップ4)。すなわち、演算装置5は、部位毎に罰則項重みβmを変更しながら、逐次近似投影データ補正処理を実行する。次式は、評価関数の一例である。
Figure JPOXMLDOC01-appb-M000003
 　ここで、本発明の技術的思想は、評価関数における罰則項重みを部位毎に変更することであるため、検出器出力重み及び罰則項を含む評価関数であれば、最適化に使用する数値解析法に因らずに適用することができる。
　式(3)から逐次近似法における更新式を導出することは、公知の手法によって実現することができる。例えば、非特許文献6のように、GauSS-Seidel法を用いることによって更新式を導出することができる。
　演算装置5は、式(3)から逐次近似法における更新式を導出し、算出された検出器出力重みdi及び真部分集合Sm毎の罰則項重みβmを、導出された更新式に代入し、逐次近似法によって、投影データの補正処理を行う。
　以降では、図6～図8を参照しながら、2種類の罰則項重み算出処理の詳細について説明する。第1の罰則項重み算出処理については、図6、図7を参照しながら説明する。第2の罰則項重み算出処理については、図8を参照しながら説明する。
　操作者は、所望計算時間や罰則項重みの推定精度を考慮して、両者のうちいずれか一方を任意に選択することができる。
　最初に、第1の罰則項重み算出処理について説明する。
　第1の罰則項重み算出処理では、演算装置5は、ノイズ低減率ごとに、基準ファントムに対してノイズ低減率と同程度のノイズ低減効果が達成される罰則項重みの値(以下、「罰則項重み基準値」と呼ぶ。)を予め記憶装置8に記憶する。
　例えば、演算装置5は、ノイズ低減率10％、20％、30％、・・・に対して、罰則項重み基準値β＊ 10、β＊ 20、β＊ 30、・・・の値が登録されているノイズ低減率テーブルを記憶装置8に記憶しておく。罰則項重み基準値β＊ 10、β＊ 20、β＊ 30、・・・の値は、予め円形の水ファントムや人体を模擬したファントムを基準ファントムとし、基準ファントムを実際に撮影し、画像再構成処理を実行して得られた値とする。
　そして、演算装置5は、記憶装置8から、入力されるノイズ低減率に対応する罰則項重み基準値を取得する。次に、演算装置5は、真部分集合Sm毎に、真部分集合Sm毎内に含まれる集合要素に対応する検出器出力重みdiに基づいて代表値を算出する。次に、演算装置5は、真部分集合Sm毎の代表値と、記憶装置8から取得される罰則項重み基準値とを乗算した値を、真部分集合Sm毎の罰則項重みβmとする。
　特に、演算装置5は、真部分集合Sm毎に、真部分集合Sm内に含まれる集合要素に対応する検出器出力重みdiの平均値及び中央値、並びに、検出器出力重みdiを階級とするヒストグラムにおいて全体を所定の割合によって分割する階級の3つの値の中でいずれか1つを、真部分集合Sm毎の代表値とすることが望ましい。
　図6に示すフローチャートでは、代表値として、検出器出力重みdiを階級とするヒストグラムにおいて全体を所定の割合によって分割する階級を算出する例を示している。
　また、図6に示すフローチャートでは、真部分集合Smに対する罰則項重みβmの算出処理が示されているが、他の真部分集合についても同様である。
　図6に示すように、演算装置5は、記憶装置8に記憶されているノイズ低減率テーブルから、操作者が入力した所望ノイズ低減率に対応するテーブル値を選択する(ステップ11)。
　次に、演算装置5は、操作者が入力した所望計算時間に応じて、真部分集合Smに含まれる集合要素から、真部分集合Smの下位集合Tmを作成する(ステップ12)。
　演算装置5は、予め定められた間引きルールに基づいて、Smに含まれる集合要素を間引くことによって、下位集合Tmを作成する。
　間引きルールとしては、例えば、ビューをとびとびに間引くことが考えられる。ビューをとびとびに間引く場合には、最終的に作成される断面像の画質に与える影響が少ないからである。例えば、周回方向に1°ずつビューを取得している場合、周回当たりの撮影ビュー数Vdは、360個となる。360個をとびとびに間引くことによって、180個のビューだけが残る為、下位集合Tmの集合要素の数は、真部分集合Smの集合要素の数の1/2となる。
　また、間引きルールとしては、例えば、各ビューの左端及び右端のチャネルを間引くことが考えられる。各ビューの左端及び右端のチャネルは、被検体を通過しないX線を受光する場合が多い為、データを間引いたとしても、最終的に作成される断面像の画質に与える影響が少ないからである。例えば、チャネル数Cに対して、左端から10％及び右端から10％のデータを間引くことによって、下位集合Tmの集合要素の数は、真部分集合Smの集合要素の数の4/5となる。
　これらの間引きルールは、組み合わせて適用することも可能である。
　いずれの間引きルール(又は複数の間引きルールの組み合わせ)であっても、操作者が入力した所望計算時間に応じて間引き量を加減することが可能である。このように、真部分集合Smの下位集合Tmを作成し、下位集合Tmに対して後述する処理を実行することによって、操作者の意図に合わせて計算時間を調整することができる。
　次に、演算装置5は、ステップ12において作成された下位集合Tmの集合要素に対応する検出器出力重みdiを階級として、ヒストグラムを構築する(ステップ13)。ここで、検出器出力重みdiの算出手法は、公知のあらゆる手法を使用することができる。演算装置5は、例えば、非特許文献6のように、投影データを逆対数変換した後、エアデータを乗算することによって算出するようにしても良い。
　図7には、ヒストグラムの一例が示されている。図7に示すヒストグラムは、横軸(階級)が検出器出力重みdi、縦軸(頻度)が各階級の値に等しい検出器出力重みdiに対応する集合要素の数である。
　ヒストグラムは、被写体の大きさや部位によって異なる分布となり、それらの特徴をよく表す指標と言える。
　肩部51のヒストグラムは、階級が0に近い値のところに最大のピークが存在する。また、肩部51のヒストグラムは、全体的に、検出器出力重みdiの値が低い。これは、肩部51が、X線の減弱量が大きい骨部を多く含むからである。
　また、胸部52及び腹部53のヒストグラムは、階級が1に近い値のところに最大のピークが存在する。胸部52及び腹部53のヒストグラムを比較すると、階級が1～3の頻度は、胸部52のヒストグラム＞腹部53のヒストグラムとなり、階級が5以上の頻度は、胸部52のヒストグラム＜腹部53のヒストグラムとなっている。これは、胸部52が、X線の減弱量が小さい肺を含むからである。
　尚、被写体の大きさが異なれば、最大のピークが存在する階級が異なるので、図7に示すヒストグラムは、被写体の大きさの特徴もよく表していると言える。
　図6の説明に戻る。次に、演算装置5は、図7に示す各ヒストグラムの面積を算出し、ヒストグラムの面積を所定の割合に分割する階級を代表値として特定する(ステップ14)。この他、代表値は、前述した通り、i∈Smに含まれる検出器出力重みdiの中央値や平均値としても良い。操作者は、これらを任意に選択することができる。
　次に、演算装置5は、ステップ11において選択されたテーブル値と、ステップ14において特定された階級(代表値)を乗算することによって、真部分集合Smの罰則項重みβmを算出する(ステップ15)。
　第1の罰則項重み算出処理では複雑な計算を行わないことから、計算量を大幅に増やすことなく、真部分集合Smの罰則項重みβmを算出することができる。
　次に、第2の罰則項重み算出処理について説明する。第2の罰則項重み算出処理は、以下の2つの知見に基づく。
(1)検出器出力重みを定数とし、逐次近似法を任意の罰則項重みによって実行した場合、被写体及び部位によらず、ほぼ同一のノイズ低減率を得られることが経験的に明らかである。
(2)逐次近似法による投影値(画素値)の修正量とノイズ低減率との間には、高い相関がある。
　尚、これらの知見は、逐次近似投影データ補正処理及び逐次近似再構成処理の両方について言える。
　前述の知見を生かして、第2の罰則項重み算出処理では、演算装置5は、逐次近似法を実行した場合の投影データの修正量を算出する修正量算出関数に対して、真部分集合Sm内に含まれる集合要素に対応する検出器出力重みdi及び投影データを代入する。ここで、修正量算出関数には、真部分集合Sm毎の罰則項重みβmが変数として含まれる。そして、演算装置5は、修正量算出関数の値(＝投影データの修正量の総和)と所定の修正量基準値(＝投影データの修正量の総和の基準値)の誤差が最小となるように、真部分集合毎Smの罰則項重みβmを決定する。
　所定の修正量基準値の決定処理は、次の通りである。演算装置5は、第1の罰則項重み算出処理と同様、ノイズ低減率テーブルを記憶装置8に記憶しておき、記憶装置8から、入力されるノイズ低減率に対応する罰則項重み基準値を取得する。そして、演算装置5は、検出器出力重みを定数とし、記憶装置8から取得される罰則項重み基準値を用いて逐次近似法を実行した場合の投影データの修正量を算出し、所定の修正量基準値とする。
　以下、図8を参照しながら、第2の罰則項重み算出処理の詳細を説明する。図8に示すフローチャートでは、真部分集合Smに対する罰則項重みβmの算出処理が示されているが、他の真部分集合についても同様である。
　図8に示すように、演算装置5は、記憶装置8に記憶されているノイズ低減率テーブルから、操作者が入力した所望ノイズ低減率に対応するテーブル値を選択する(ステップ21)。尚、処理内容の違いから、第1の罰則項重み算出処理に用いるノイズ低減率テーブルと、第2の罰則項重み算出処理に用いるノイズ低減率テーブルとは、記憶されているテーブル値が異なる。
　次に、演算装置5は、操作者が入力した所望計算時間に応じて、真部分集合Smに含まれる集合要素から、真部分集合Smの下位集合Tmを作成する(ステップ22)。下位集合Tmの作成処理は、ステップ12と同様である。
　次に、演算装置5は、ステップ21において選択されたテーブル値、及び、ステップ22において作成された下位集合Tmの集合要素に対応する投影データを用いて、検出器出力重みdiを任意の定数とした上で、逐次近似投影データ補正処理の適用後の投影値を推定する。ここで、計算量及び計算に要するメモリの観点から、適用後の投影値を解析的に算出することが難しいため、近似値を代用する。注目する投影値が、近接するチャネル、列、ビューのみから算出されると仮定すると、近似値は容易に算出できる。また、公知の方法を用いて、反復後の投影値の算出に必要となる逆行列を近似逆行列に置き換えた上で、近似値を算出しても良い。そして、演算装置5は、近似値と、逐次近似投影データ補正処理の適用前の投影値との誤差を算出し、誤差の総和を修正量基準値とする(ステップ23)。ここで、誤差の算出処理には、絶対誤差や二乗誤差等公知の指標を用いることができる。
　次に、演算装置5は、罰則項重みβmを変数とした修正量算出関数に、ステップ23において作成された下位集合Tmの集合要素に対応する検出器出力重みdi及び投影データを代入する。そして、演算装置5は、修正量算出関数の値が、ステップ23において算出された修正量基準値と等しくなるように、罰則項重みを決定する(ステップ24)。ここで、修正量算出関数の値と修正量基準値との誤差の最小化には、公知の数値解析法(例えば二分法)を適用することができる。
　第2の罰則項重み算出処理は前述の知見に基づくものであることから、精度良く真部分集合Smの罰則項重みβmを算出することができる。
　以上のように、第1の実施の形態では、逐次近似投影データ補正処理に用いられる罰則項重みβmを、真部分集合Sm内では一定値とすることによって、逐次近似投影データ補正処理において、従来と同様に、CT画像のノイズ低減効果を達成することができる。
　また、被写体の情報が反映された検出器出力重みdi及び投影値に基づいて真部分集合Sm毎に罰則項重みβmを算出することによって、部位間でのノイズ低減効果のばらつきを抑制することができる。
　＜第2の実施の形態＞
　第2の実施の形態について説明する。第2の実施の形態では、本発明を逐次近似再構成処理に適用する場合について説明する。つまり、第2の実施の形態では、演算装置5が、検出器出力重み及び罰則項を評価関数に含む逐次近似法によって、再構成画像の再構成処理を行う。
　第1の実施の形態と第2の実施の形態の違いは、第1の実施の形態における式(3)が、次式になることのみである。
Figure JPOXMLDOC01-appb-M000004
 　式(4)についても、式(3)と同様、式(4)から逐次近似法における更新式を導出することができる。演算装置5は、第1の実施の形態と同様に算出された検出器出力重みdi及び真部分集合Sm毎の罰則項重みβmを、導出された更新式に代入し、逐次近似法によって、再構成画像の再構成処理を行う。
　以上のように、第2の実施の形態では、逐次近似再構成処理に用いられる罰則項重みβmを、第1の実施の形態と同様、真部分集合Sm内では一定値とすることによって、逐次近似再構成処理において、従来と同様に、CT画像のノイズ低減効果を達成することができる。
　また、第1の実施の形態と同様、被写体の情報が反映された検出器出力重みdi及び投影値に基づいて真部分集合Sm毎に罰則項重みβmを算出することによって、部位間でのノイズ低減効果のばらつきを抑制することができる。
　＜第3の実施の形態＞
　第3の実施の形態について説明する。第3の実施の形態では、本発明を逐次近似投影データ補正処理及び逐次近似再構成処理の両方に適用する場合について説明する。つまり、第3の実施の形態では、演算装置5が、検出器出力重み及び罰則項を評価関数に含む逐次近似法によって、投影データの補正処理及び再構成画像の再構成処理を行う。
　第3の実施の形態は、第1の実施の形態と第2の実施の形態の組み合わせとなる。つまり、演算装置5は、第1の実施の形態と同様、式(3)から逐次近似法における更新式を導出し、第1の実施の形態と同様に算出された検出器出力重みdi及び真部分集合Sm毎の罰則項重みβmを、導出された更新式に代入し、逐次近似法によって、投影データの補正処理を行う。また、演算装置5は、第2の実施の形態と同様、式(4)から逐次近似法における更新式を導出し、第1の実施の形態と同様に算出された検出器出力重みdi及び真部分集合Sm毎の罰則項重みβmを、導出された更新式に代入し、逐次近似法によって、再構成画像の再構成処理を行う。
　以上のように、第3の実施の形態では、逐次近似投影データ補正処理及び逐次近似再構成処理に用いられる罰則項重みβmを、真部分集合Sm内では一定値とすることによって、逐次近似投影データ補正処理及び逐次近似再構成処理において、従来と同様に、CT画像のノイズ低減効果を達成することができる。
　また、真部分集合Smは、被写体の情報が反映された検出器出力重みdi及び投影値に基づいて真部分集合Sm毎に罰則項重みβmを算出することによって、部位間でのノイズ低減効果のばらつきを抑制することができる。
　以上、添付図面を参照しながら、本発明に係る医用画像処理装置等の好適な実施形態について説明したが、本発明はかかる例に限定されない。当業者であれば、本願で開示した技術的思想の範疇内において、各種の変更例又は修正例に想到し得ることは明らかであり、それらについても当然に本発明の技術的範囲に属するものと了解される。
　1　X線CT装置、　4　寝台、　5　演算装置、　6　入力装置、　7　表示装置、　8　記憶装置、　11　X線管、　12　検出器、　41　断面像、　42　チャネル、　43冠状面像、　44　列、　45スキャン方向

Claims (8)

  1. 　検出器の出力に応じた重みである検出器出力重み、及び罰則項を評価関数に含む逐次近似法によって、投影データの補正処理及び／又は再構成画像の再構成処理を行う医用画像処理装置であって、
    　前記検出器のチャネル、前記検出器の列、及び前記投影データの取得単位であるビューの組み合わせを一意に識別する情報を集合要素とする全体集合から、撮影条件及び再構成条件に基づいて1又は複数の真部分集合を決定する真部分集合決定部と、
    　前記真部分集合毎に、前記真部分集合内に含まれる前記集合要素に対応する前記検出器出力重みに基づいて、前記罰則項に係る重みである罰則項重みを算出する罰則項重み算出部と、
    　前記真部分集合毎の前記罰則項重みを用いて前記逐次近似法を実行する逐次近似法実行部と、
    　を備える医用画像処理装置。
  2. 　前記真部分集合決定部は、らせんスキャンの場合、前記逐次近似法における逆投影処理の計算単位として定められる前記ビューの数である逆投影ビュー数に基づいて、前記真部分集合に含まれる前記集合要素の数を決定する
    請求項1に記載の医用画像処理装置。
  3. 　前記真部分集合決定部は、
    　寝台送り速度及び前記再構成画像の大きさごとに、前記逐次近似法における逆投影処理の位相幅である逆投影位相幅を予め記憶部に記憶し、
    　前記記憶部から、前記撮影条件として設定される寝台送り速度、及び前記再構成条件として設定される前記再構成画像の大きさに対応する前記逆投影位相幅を取得し、
    　前記撮影条件として設定される周回当たりの撮影ビュー数と、前記記憶部から取得される前記逆投影位相幅とを乗算した値を前記逆投影ビュー数とする
    　請求項2に記載の医用画像処理装置。
  4. 　前記罰則項重み算出部は、
    　ノイズ低減率ごとに、基準ファントムに対して前記ノイズ低減率と同程度のノイズ低減が達成される前記罰則項重みである罰則項重み基準値を予め記憶部に記憶し、
    　前記記憶部から、入力される前記ノイズ低減率に対応する前記罰則項重み基準値を取得し、
    　前記真部分集合毎に、前記真部分集合内に含まれる前記集合要素に対応する前記検出器出力重みに基づいて代表値を算出し、
    　前記真部分集合毎の前記代表値と、前記記憶部から取得される前記罰則項重み基準値とを乗算した値を、前記真部分集合毎の前記罰則項重みとする
    　請求項1に記載の医用画像処理装置。
  5. 　前記罰則項重み算出部は、前記真部分集合毎に、前記真部分集合内に含まれる前記集合要素に対応する前記検出器出力重みの平均値及び中央値、並びに、前記検出器出力重みを階級とするヒストグラムにおいて全体を所定の割合によって分割する階級の3つの値の中でいずれか1つを、前記真部分集合毎の前記代表値とする
    　請求項4に記載の医用画像処理装置。
  6. 　前記罰則項重み算出部は、前記真部分集合毎の前記罰則項重みが変数であって、前記逐次近似法を実行した場合の前記投影データの修正量を算出する修正量算出関数に対して、前記真部分集合内に含まれる前記集合要素に対応する前記検出器出力重み及び前記投影データを代入して算出される前記修正量算出関数の値と修正量基準値との誤差が最小となるように、前記真部分集合毎の前記罰則項重みを決定する
    　請求項1に記載の医用画像処理装置。
  7. 　前記罰則項重み算出部は、
    　ノイズ低減率ごとに、基準ファントムに対して前記ノイズ低減率と同程度のノイズ低減効果が達成される前記罰則項重みである罰則項重み基準値を予め記憶部に記憶し、
    　前記記憶部から、入力される前記ノイズ低減率に対応する前記罰則項重み基準値を取得し、
    　前記検出器出力重みを定数とし、前記記憶部から取得される前記罰則項重み基準値を用いて前記逐次近似法を実行した場合の前記投影データの修正量を算出し、前記修正量基準値とする
    　請求項6に記載の医用画像処理装置。
  8. 　検出器の出力に応じた重みである検出器出力重み、及び罰則項を評価関数に含む逐次近似法によって、投影データの補正処理及び／又は再構成画像の再構成処理を行う医用画像処理方法であって、
    　前記検出器のチャネル、前記検出器の列、及び前記投影データの取得単位であるビューの組み合わせを一意に識別する情報を集合要素とする全体集合から、撮影条件及び再構成条件に基づいて1又は複数の真部分集合を決定するステップと、
    　前記真部分集合毎に、前記真部分集合内に含まれる前記集合要素に対応する前記検出器出力重みに基づいて、前記罰則項に係る重みである罰則項重みを算出するステップと、
    　前記真部分集合毎の前記罰則項重みを用いて前記逐次近似法を実行するステップと、
    　を含む医用画像処理方法。
PCT/JP2012/059136 2011-04-28 2012-04-04 医用画像処理装置、医用画像処理方法 WO2012147471A1 (ja)

Priority Applications (3)

Application Number Priority Date Filing Date Title
US14/112,276 US9123098B2 (en) 2011-04-28 2012-04-04 Medical image processing device and medical image processing method, applying weighted penalty term in iterative approximation
CN201280020678.5A CN103501702B (zh) 2011-04-28 2012-04-04 医用图像处理装置、医用图像处理方法
JP2013511984A JP5978429B2 (ja) 2011-04-28 2012-04-04 医用画像処理装置、医用画像処理方法

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2011-100571 2011-04-28
JP2011100571 2011-04-28

Publications (1)

Publication Number Publication Date
WO2012147471A1 true WO2012147471A1 (ja) 2012-11-01

Family

ID=47071996

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2012/059136 WO2012147471A1 (ja) 2011-04-28 2012-04-04 医用画像処理装置、医用画像処理方法

Country Status (4)

Country Link
US (1) US9123098B2 (ja)
JP (1) JP5978429B2 (ja)
CN (1) CN103501702B (ja)
WO (1) WO2012147471A1 (ja)

Cited By (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2014128576A (ja) * 2012-11-30 2014-07-10 Canon Inc 画像処理装置、画像処理方法、及びプログラム
WO2014123041A1 (ja) * 2013-02-05 2014-08-14 株式会社 日立メディコ Ｘ線ｃｔ装置及び画像再構成方法
WO2015108097A1 (ja) * 2014-01-20 2015-07-23 株式会社 日立メディコ Ｘ線ｃｔ装置、画像処理装置、及び画像再構成方法
WO2015137011A1 (ja) * 2014-03-14 2015-09-17 株式会社日立メディコ Ｘ線ｃｔ装置、及び処理装置
CN105188543A (zh) * 2013-04-08 2015-12-23 株式会社日立医疗器械 X射线ct装置、重构运算装置以及重构运算方法
CN105451658A (zh) * 2013-08-08 2016-03-30 株式会社日立医疗器械 X射线ct装置以及校正处理装置
JP2016188786A (ja) * 2015-03-30 2016-11-04 三菱電機株式会社 窓関数決定装置、パルス圧縮装置、レーダ信号解析装置、レーダ装置、窓関数決定方法およびプログラム
WO2016199716A1 (ja) * 2015-06-12 2016-12-15 株式会社日立製作所 Ｘ線ｃｔ装置および逐次修正パラメータ決定方法
CN110070588A (zh) * 2019-04-24 2019-07-30 上海联影医疗科技有限公司 Pet图像重建方法、系统、可读存储介质和设备

Families Citing this family (14)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9123098B2 (en) * 2011-04-28 2015-09-01 Hitachi Medical Corporation Medical image processing device and medical image processing method, applying weighted penalty term in iterative approximation
WO2013008702A1 (ja) * 2011-07-08 2013-01-17 株式会社 日立メディコ 画像再構成装置及び画像再構成方法
IN2014DN05824A (ja) * 2012-04-24 2015-05-15 Hitachi Medical Corp
JP6218334B2 (ja) * 2012-11-30 2017-10-25 株式会社日立製作所 Ｘ線ｃｔ装置及びその断層画像撮影方法
CN104318536B (zh) 2014-10-21 2018-03-20 沈阳东软医疗系统有限公司 Ct图像的校正方法及装置
JPWO2016132880A1 (ja) * 2015-02-16 2017-11-30 株式会社日立製作所 演算装置、ｘ線ｃｔ装置、及び画像再構成方法
JP2017122705A (ja) * 2016-01-06 2017-07-13 三菱電機株式会社 算出方法、判定方法、選別方法および選別装置
CN106994021B (zh) * 2016-01-22 2022-10-14 通用电气公司 一种计算ct影像上的噪声的方法及装置
CN108780052B (zh) * 2016-03-11 2020-11-17 株式会社岛津制作所 图像重构处理方法、图像重构处理程序以及安装有该程序的断层摄影装置
US10650621B1 (en) 2016-09-13 2020-05-12 Iocurrents, Inc. Interfacing with a vehicular controller area network
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
JP7282487B2 (ja) * 2018-06-07 2023-05-29 キヤノンメディカルシステムズ株式会社 医用画像診断装置
JP7077208B2 (ja) 2018-11-12 2022-05-30 富士フイルムヘルスケア株式会社 画像再構成装置および画像再構成方法
US10679385B1 (en) 2018-12-17 2020-06-09 General Electric Company System and method for statistical iterative reconstruction and material decomposition

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPS6017568A (ja) * 1983-07-11 1985-01-29 Hitachi Ltd 画像処理方法および装置
JP2009540966A (ja) * 2006-06-22 2009-11-26 ゼネラル・エレクトリック・カンパニイ 画像の分解能を高めるシステム及び方法
WO2010062956A1 (en) * 2008-11-26 2010-06-03 Wisconsin Alumni Research Foundation Method for prior image constrained image reconstruction in cardiac cone beam computed tomography
JP2010136958A (ja) * 2008-12-13 2010-06-24 Univ Of Tokushima Ｃｔ装置、ｃｔ装置における画像再構成方法、及び電子回路部品

Family Cites Families (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5909476A (en) * 1997-09-22 1999-06-01 University Of Iowa Research Foundation Iterative process for reconstructing cone-beam tomographic images
US8204173B2 (en) * 2003-04-25 2012-06-19 Rapiscan Systems, Inc. System and method for image reconstruction by using multi-sheet surface rebinning
WO2005077278A1 (ja) * 2004-02-16 2005-08-25 Hitachi Medical Corporation 断層撮影像の再構成方法及び断層撮影装置
EP1861825B1 (en) * 2005-03-17 2008-10-08 Philips Intellectual Property & Standards GmbH Method and device for the iterative reconstruction of cardiac images
CN102105106B (zh) * 2008-08-07 2013-12-25 株式会社日立医疗器械 X射线ct图像形成方法和应用了该方法的x射线ct装置
US20110052023A1 (en) * 2009-08-28 2011-03-03 International Business Machines Corporation Reconstruction of Images Using Sparse Representation
CN101901472B (zh) * 2010-07-07 2012-12-19 清华大学 一种基于矩阵秩最小化的非刚性鲁棒批量图像对齐方法
CN103153192B (zh) * 2010-10-14 2015-09-23 株式会社日立医疗器械 X射线ct装置以及图像再构成方法
US9123098B2 (en) * 2011-04-28 2015-09-01 Hitachi Medical Corporation Medical image processing device and medical image processing method, applying weighted penalty term in iterative approximation

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPS6017568A (ja) * 1983-07-11 1985-01-29 Hitachi Ltd 画像処理方法および装置
JP2009540966A (ja) * 2006-06-22 2009-11-26 ゼネラル・エレクトリック・カンパニイ 画像の分解能を高めるシステム及び方法
WO2010062956A1 (en) * 2008-11-26 2010-06-03 Wisconsin Alumni Research Foundation Method for prior image constrained image reconstruction in cardiac cone beam computed tomography
JP2010136958A (ja) * 2008-12-13 2010-06-24 Univ Of Tokushima Ｃｔ装置、ｃｔ装置における画像再構成方法、及び電子回路部品

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
H.R.SHI ET AL.: "Quadratic Regularization Design for 2-D CT", IEEE TRANS.MED.IMG., vol. 28, no. 5, May 2009 (2009-05-01), pages 645 - 656 *
J.WANG ET AL.: "Multiscale Penalized Weighted Least-Squares Sinogram Restoration for Low-Dose X-Ray Computed Tomography", CONF.PROC.IEEE ENG.MED.BIOL.SOC., vol. 55, no. 3, March 2008 (2008-03-01), pages 1022 - 1031 *
TAIGA GOTO ET AL.: "Chikuji Kinji Saikosei Algorithm no Kaihatsu", JAPANESE SOCIETY OF RADIOLOGICAL TECHNOLOGY SOKAI GAKUJUTSU TAIKAI YOKOSHU, vol. 67, 25 February 2011 (2011-02-25), pages 184 *
Z.TIAN ET AL.: "Low-dose CT reconstruction via edge-preserving total variation regularization.", PHYS.MED.BIOL., vol. 56, no. 18, September 2011 (2011-09-01), pages 5949 - 5967 *

Cited By (20)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2014128576A (ja) * 2012-11-30 2014-07-10 Canon Inc 画像処理装置、画像処理方法、及びプログラム
JPWO2014123041A1 (ja) * 2013-02-05 2017-02-02 株式会社日立製作所 Ｘ線ｃｔ装置及び画像再構成方法
WO2014123041A1 (ja) * 2013-02-05 2014-08-14 株式会社 日立メディコ Ｘ線ｃｔ装置及び画像再構成方法
US9715745B2 (en) 2013-02-05 2017-07-25 Hitachi, Ltd. X-ray CT apparatus and image reconstruction method
CN104902818A (zh) * 2013-02-05 2015-09-09 株式会社日立医疗器械 X射线ct装置以及图像重构方法
CN105188543A (zh) * 2013-04-08 2015-12-23 株式会社日立医疗器械 X射线ct装置、重构运算装置以及重构运算方法
CN105451658A (zh) * 2013-08-08 2016-03-30 株式会社日立医疗器械 X射线ct装置以及校正处理装置
JPWO2015108097A1 (ja) * 2014-01-20 2017-03-23 株式会社日立製作所 Ｘ線ｃｔ装置、画像処理装置、及び画像再構成方法
WO2015108097A1 (ja) * 2014-01-20 2015-07-23 株式会社 日立メディコ Ｘ線ｃｔ装置、画像処理装置、及び画像再構成方法
US9974495B2 (en) 2014-01-20 2018-05-22 Hitachi, Ltd. X-ray CT apparatus, image processing device, and image reconstruction method
WO2015137011A1 (ja) * 2014-03-14 2015-09-17 株式会社日立メディコ Ｘ線ｃｔ装置、及び処理装置
JPWO2015137011A1 (ja) * 2014-03-14 2017-04-06 株式会社日立製作所 Ｘ線ｃｔ装置、及び処理装置
US20170119335A1 (en) * 2014-03-14 2017-05-04 Hitachi, Ltd. X-ray ct device and processing device
US10368824B2 (en) 2014-03-14 2019-08-06 Hitachi, Ltd. X-ray CT device and processing device
JP2016188786A (ja) * 2015-03-30 2016-11-04 三菱電機株式会社 窓関数決定装置、パルス圧縮装置、レーダ信号解析装置、レーダ装置、窓関数決定方法およびプログラム
WO2016199716A1 (ja) * 2015-06-12 2016-12-15 株式会社日立製作所 Ｘ線ｃｔ装置および逐次修正パラメータ決定方法
JPWO2016199716A1 (ja) * 2015-06-12 2018-03-01 株式会社日立製作所 Ｘ線ｃｔ装置および逐次修正パラメータ決定方法
US10210633B2 (en) 2015-06-12 2019-02-19 Hitachi, Ltd. X-ray CT device and sequential correction parameter determination method
CN110070588A (zh) * 2019-04-24 2019-07-30 上海联影医疗科技有限公司 Pet图像重建方法、系统、可读存储介质和设备
CN110070588B (zh) * 2019-04-24 2023-01-31 上海联影医疗科技股份有限公司 Pet图像重建方法、系统、可读存储介质和设备

Also Published As

Publication number Publication date
US9123098B2 (en) 2015-09-01
US20140226887A1 (en) 2014-08-14
CN103501702A (zh) 2014-01-08
JPWO2012147471A1 (ja) 2014-07-28
JP5978429B2 (ja) 2016-08-24
CN103501702B (zh) 2015-09-02

Similar Documents

Publication Publication Date Title
JP5978429B2 (ja) 医用画像処理装置、医用画像処理方法
JP6937157B2 (ja) 放射線画像診断装置及び医用画像処理装置
Liu et al. Total variation-stokes strategy for sparse-view X-ray CT image reconstruction
JP5530637B2 (ja) 画像再構成の方法及びシステム
JP5828841B2 (ja) Ｘ線ｃｔ装置及び画像再構成方法
JP6492005B2 (ja) Ｘ線ｃｔ装置、再構成演算装置、及び再構成演算方法
JP6470837B2 (ja) Ｘ線ｃｔ装置および逐次修正パラメータ決定方法
JP5968316B2 (ja) 画像再構成装置及び画像再構成方法
US8885903B2 (en) Method and apparatus for statistical iterative reconstruction
JP6215449B2 (ja) Ｘ線ｃｔ装置、及び処理装置
WO2014041889A1 (ja) Ｘ線ｃｔ装置およびｘ線ｃｔ画像の処理方法
US10813616B2 (en) Variance reduction for monte carlo-based scatter estimation
JP2016152916A (ja) Ｘ線コンピュータ断層撮像装置及び医用画像処理装置
JPWO2014123041A1 (ja) Ｘ線ｃｔ装置及び画像再構成方法
JP6124868B2 (ja) 画像処理装置及び画像処理方法
JP5637768B2 (ja) コンピュータ断層撮影画像の生成方法およびコンピュータ断層撮影装置
JP2018110867A (ja) 医用画像生成装置及び医用画像生成方法
JP2015231528A (ja) Ｘ線コンピュータ断層撮像装置及び医用画像処理装置
WO2016132880A1 (ja) 演算装置、ｘ線ｃｔ装置、及び画像再構成方法
US20220375038A1 (en) Systems and methods for computed tomography image denoising with a bias-reducing loss function
US9508164B2 (en) Fast iterative image reconstruction method for 3D computed tomography
US11353411B2 (en) Methods and systems for multi-material decomposition

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

Country of ref document: EP

Kind code of ref document: A1

ENP Entry into the national phase

Ref document number: 2013511984

Country of ref document: JP

Kind code of ref document: A

WWE Wipo information: entry into national phase

Ref document number: 14112276

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

Country of ref document: EP

Kind code of ref document: A1