WO2013146283A1 - 画像処理装置及び画像処理方法 - Google Patents

画像処理装置及び画像処理方法 Download PDF

Info

Publication number
WO2013146283A1
WO2013146283A1 PCT/JP2013/057118 JP2013057118W WO2013146283A1 WO 2013146283 A1 WO2013146283 A1 WO 2013146283A1 JP 2013057118 W JP2013057118 W JP 2013057118W WO 2013146283 A1 WO2013146283 A1 WO 2013146283A1
Authority
WO
WIPO (PCT)
Prior art keywords
image
region
image processing
feature amount
smoothed
Prior art date
Application number
PCT/JP2013/057118
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/372,330 priority Critical patent/US9454804B2/en
Priority to IN5834DEN2014 priority patent/IN2014DN05834A/en
Priority to CN201380005062.5A priority patent/CN104066378B/zh
Priority to JP2014507658A priority patent/JP6124868B2/ja
Publication of WO2013146283A1 publication Critical patent/WO2013146283A1/ja

Links

Images

Classifications

    • 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/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
    • 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/008—Specific post-processing after tomographic reconstruction, e.g. voxelisation, metal artifact correction
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/50—Image enhancement or restoration using two or more images, e.g. averaging or subtraction
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/77—Retouching; Inpainting; Scratch removal
    • 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/20212—Image combination
    • G06T2207/20224—Image subtraction
    • 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

Definitions

  • the present invention relates to an image processing apparatus that performs image processing on an image obtained by an X-ray CT apparatus or the like, and relates to a technique for removing streaky artifacts while maintaining the edge of a structure.
  • the X-ray CT apparatus reconstructs a tomographic image of a subject by irradiating X-rays from around the subject and backprojecting actual projection data acquired at a plurality of projection angles.
  • Non-Patent Document 1 and Non-Patent Document 2 propose a method of reducing streak-like artifacts by applying the above-described image reconstruction method after smoothing projection data.
  • a method for reducing these streak-like artifacts is referred to as a “projection data smoothing filter”.
  • the image obtained by the projection data smoothing filter is blurred at the edge of a relatively high-contrast structure accompanying the streak-like artifact reduction effect.
  • the effect of reducing streak-like artifacts and the effect of preserving the edge of the structure are in a trade-off relationship, and the latter effect decreases as the former effect increases. Therefore, in the case of reducing strong streak-like artifacts, large blurring occurs at the edges of the structure. Since the edge of the structure is not preferable in clinical diagnosis, an image in which the edge of the structure is maintained and the streak-like artifact is removed (hereinafter referred to as “target image”) is desired.
  • an image obtained from the projection data by the above-described image reconstruction method (hereinafter referred to as “original image”) and the above-described image reconstruction method from the projection data to which the projection data smoothing filter is applied.
  • original image an image obtained from the projection data by the above-described image reconstruction method
  • smoothed images images obtained by the above
  • weighted addition is performed so that the ratio of the pixel value of the smoothed image is increased at the pixel portion where the streak-like artifact is generated on the image and the ratio of the pixel value of the original image is increased at the edge portion of the structure.
  • a low-resolution image having a lower resolution than the original image is created from the original image, an edge is detected from the low-resolution image, and a low-resolution edge image is created.
  • a weighting coefficient is calculated from the edge luminance value by a non-linear function.
  • the difference between the maximum value and the minimum value of the pixels in a predetermined region including the target pixel, the variance value, the density gradient, and the like are used for edge detection.
  • the edge of the structure is emphasized by weighting and adding the image obtained by subtracting the low-resolution image from the original image using the weighting coefficient described above.
  • the image is divided into small areas, the standard deviation value of the small area is calculated, a histogram with the standard deviation value as a class is created for all small areas, and the discrimination class value is determined from the histogram. And each small region is classified into a flat region and a non-flat region. Detecting the edge of the structure on the image is equivalent to identifying flat and non-flat areas. If both can be classified, the pixel value of the original image is used in the non-flat region, and the pixel value of the smoothed image is used in the flat region, so that the weighted addition of the original image and the smoothed image can be performed.
  • the strong streak-like artifact has a characteristic close to a linear structure unlike the granular noise generated in the image. Then, as in the method described in Patent Document 1, when performing edge detection using the difference between the maximum value and the minimum value of the pixel value of the predetermined region including the target pixel, the variance value, the density gradient, A strong streak-like artifact is likely to be detected as an edge.
  • Patent Document 1 and Patent Document 2 are intended to classify noise and structures in a uniform region, and streak-like artifacts. This is because the index for classifying the structure is not used.
  • the present invention has been made in view of the above-described problems, and an object of the present invention is to generate a target image in which the edge of a structure is maintained and streak-like artifacts are removed.
  • An image processing apparatus is provided.
  • an original image generated from original information and corresponding smoothed images generated from the same original information and at least streak-like artifacts are reduced.
  • An image processing apparatus for generating a target image in which edges of a structure are maintained and streak-like artifacts are removed by weight-adding pixels using a weighting coefficient, wherein the feature amounts of the original image and the smoothed image An image processing apparatus comprising: a function shape determining unit that determines a shape of a nonlinear function based on the first step; and a weighting factor calculating unit that calculates the weighting factor based on the nonlinear function.
  • the second invention weights and adds the corresponding pixels of the original image generated from the original information and the smoothed image generated from the same original information and having at least streak-like artifacts reduced by a weighting coefficient.
  • An image processing method comprising: a function shape determining step for determining; and a weighting coefficient calculating step for calculating the weighting coefficient based on the nonlinear function.
  • an image processing apparatus capable of generating a target image in which the edge of a structure is maintained and streak-like artifacts are removed.
  • 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 and the like are included.
  • 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 source), a detector 12 (X-ray detector), a collimator 13, a drive unit 14, a central control unit 15, an X-ray control unit 16, a high voltage generation unit 17, and a scanner control.
  • the device 18 includes a bed control device 19, a bed movement measurement device 20, a collimator control device 21, a preamplifier 22, an A / D converter 23, and the like.
  • 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.
  • 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.
  • 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 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 superimposes the reconstruction filter on the actual projection data of each view to generate filter-corrected projection data, and performs back-projection processing with weighting the filter-corrected projection data in the view direction.
  • 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.
  • each pixel corresponding to the original image generated from the original information acquired by the X-ray CT apparatus 1 and the like and the smoothed image generated from the same original information and at least streak-like artifacts are reduced.
  • An image processing method for generating a target image in which the edge of the structure is maintained and streak-like artifacts are removed by weighted addition using a weighting coefficient will be described.
  • the image processing method of the present invention generates a target image so as to satisfy the following three constraints.
  • the pixel value of the target image is within the range of the pixel value of the original image and the pixel value of the smoothed image related to the corresponding pixel.
  • the target image is a smooth image spatially (that is, in one image space).
  • the constraint condition 1 is the original information acquired by the X-ray CT apparatus 1 (excluding noise). ) Is a condition for not losing.
  • the constraint condition 1 is min (pixel value of the original image, pixel value of the smoothed image) ⁇ pixel value of the target image ⁇ max (pixel value of the original image, pixel value of the smoothed image).
  • Min () is an operator that outputs a minimum value
  • max () is an operator that outputs a maximum value.
  • the constraint condition 2 is a condition for avoiding discontinuity that may occur in the target image.
  • the constraint condition 3 is a condition for mainly removing graininess noise and streak-like artifacts.
  • the image processing method of the present invention is particularly characterized in that the constraint condition 3 is satisfied.
  • the weighting coefficient used for weighted addition is determined by the state coefficient representing the state of the original image and the smoothed image.
  • the state coefficient is a value of a nonlinear function whose variable is the pixel value difference between the target pixel of the original image or the smoothed image and its neighboring pixels.
  • the nonlinear function is a smooth continuous function, has a plurality of arbitrary parameters, and can be arbitrarily adjusted in shape, and is not particularly limited. If the non-linear function is a smooth continuous function, it satisfies the above-mentioned “(constraint condition 2) that the target image is a smooth image spatially (that is, in one image space)”. Another example of the nonlinear function is a logistic function.
  • a state coefficient g s (x) of the target pixel s is calculated by a generalized Gaussian function. It is defined by the following formula.
  • p is an arbitrary parameter for adjusting the slope of the generalized Gaussian function, and is set to the same value for all slices.
  • ⁇ is an arbitrary parameter for adjusting the bending position of the generalized Gaussian function.
  • w sr is a weighting coefficient corresponding to the distance between the target pixel s and the neighboring pixel r, and is defined as, for example, the following equation.
  • l sr is the distance between the target pixel s and the neighboring pixel r.
  • L is the ratio of the pixel size to the detector element size.
  • w sr need not be limited to Equation (2), and may be a function that increases in value as the adjacent pixel r is closer to the target pixel s.
  • the slope p of the set N s, and generalized Gaussian function of neighboring pixels of the pixel of interest s shall be determined empirically.
  • FIG. 3 illustrates changes in the function shape when ⁇ is set to arbitrary constants c, 2c, and 4c with respect to the pixel value difference between the target pixel s and the neighboring pixel r. If the pixel value difference is fixed, the state coefficient is likely to be close to 1 as ⁇ is increased.
  • the state coefficient is a value close to approximately 1.
  • a parameter ⁇ for adjusting the bending position of the function is determined.
  • the computing device 5 (hereinafter referred to as “computing device 5”) of the X-ray CT apparatus 1 determines the shape of the nonlinear function based on the feature amounts of the original image and the smoothed image. The process for determining the shape of the nonlinear function will be described with reference to FIG.
  • the arithmetic device 5 divides the area of the original image and the smoothed image into small areas corresponding to each other. For example, as shown in FIG. 6, the arithmetic device 5 divides the regions of the original image 41 and the smoothed image 42 into a grid and divides it into small regions. The size of the small area is determined empirically.
  • the method for dividing the small area is not limited to this example.
  • the area of the original image 41 and the smoothed image 42 may be divided into small areas by the same division method.
  • the shape of the small area is not limited to a rectangle.
  • a plurality of adjacent pixels may be included in the same small area. Further, one pixel may be included in a plurality of small regions.
  • the arithmetic device 5 calculates the variation value of the original image 41 and the smoothed image 42 from the pixel values of the pixels included in the small region for each small region related to both the original image 41 and the smoothed image 42.
  • the variation value is, for example, a standard deviation value or (maximum value ⁇ minimum value).
  • the variation value is a value indicating variation in each small region, and may be a statistic calculated from the pixel values of the pixels included in each small region. In the following, in order to avoid confusion, the standard deviation value will be described as an example.
  • the arithmetic device 5 calculates the variation value reduction rate of the smoothed image 42 based on the variation value of the original image 41 for each small region.
  • the reduction rate of the standard deviation value is calculated.
  • the reduction rate of the standard deviation value of the smoothed image 42 based on the standard deviation value of the original image 41 is set as ⁇ i and is calculated by the following equation.
  • the arithmetic device 5 calculates the standard deviation value reduction rate ⁇ i in all the small regions based on the equation (3).
  • equation (3) the shape of the nonlinear function is determined using the feature amount (for example, the reduction rate ⁇ i ) calculated from the pixel values of both the original image 41 and the smoothed image 42, and is shown in FIG.
  • Patent Documents 1 and 2 related to a method of mixing two images having the same original information
  • Patent Documents 1 and 2 related to a method of mixing two images having the same original information
  • only feature values calculated from pixel values of pixels included in only one image are used.
  • a standard deviation value calculated from pixel values of pixels included in only one image it is impossible to distinguish whether noise has been reduced or a structure has been blurred by an arbitrary smoothing process. Therefore, the conventional technique cannot generate a target image from which streak-like artifacts are removed as in the present invention.
  • the arithmetic device 5 extracts a feature amount calculation region from a set of small regions based on the reduction rate ⁇ i .
  • the arithmetic device 5 extracts a small region having the maximum reduction rate as a feature amount calculation region from small regions included in a predetermined range (for example, in the same slice).
  • the predetermined range may be within a plurality of slices or all slices.
  • the arithmetic device 5 determines the shape of the nonlinear function from the feature amount calculated from the pixel value of the feature amount calculation region (kth small region).
  • the feature amount is calculated using the standard deviation value ⁇ k (org) of the original image 41 calculated in step S202 and the standard deviation value ⁇ k (smt) of the smoothed image 42 calculated in step S203.
  • ⁇ k (org) the standard deviation value of the original image 41 calculated in step S202
  • ⁇ k (smt) of the smoothed image 42 calculated in step S203.
  • Arbitrary constant is set to ⁇ (0 ⁇ ⁇ 1), and when the pixel value difference is ⁇ k (generic symbol of ⁇ k (org) and ⁇ k (smt)) , the state coefficient takes the value of ⁇ .
  • ⁇ is a real value such as 0.99, 0.98, 0.97, for example.
  • the parameter ⁇ for adjusting the bending position of the function is defined as the following expression.
  • Equation (4) parameters for adjusting the bending position of the nonlinear function calculated in each of the original image 41 and the smoothed image 42 are set as ⁇ k (org) and ⁇ k (smt) .
  • the arithmetic device 5 dynamically determines the shape of the nonlinear function for each predetermined range (for example, for each slice) by the processing shown in FIG.
  • One of the features in the first embodiment is that a region where the reduction rate ⁇ i of the standard deviation value of the smoothed image 42 with respect to the original image 41 is maximized is extracted as a feature amount calculation region, and the standard deviation in the feature amount calculation region is extracted. The point is to determine the shape of the nonlinear function using the value.
  • the projection data smoothing filter described above has a characteristic of greatly reducing streak-like artifacts. From this characteristic, the standard deviation value reduction rate ⁇ i becomes large in a small region including streak-like artifacts. Further, the projection data smoothing filter has a characteristic that it is accompanied by blurring of the structure. From this characteristic, the standard deviation value of the original image increases in a small region including a structure, and the denominator of Equation (3) increases, resulting in a reduction in standard deviation value reduction rate ⁇ i . That is, there is a high possibility that streak-like artifacts are included in the region where the standard deviation value reduction rate ⁇ i is maximum.
  • Step S102> The arithmetic device 5 calculates the state coefficients of the original image 41 and the smoothed image 42 using a nonlinear function whose shape is determined in step S101.
  • the processing for calculating the state coefficient follows equation (1).
  • the arithmetic device 5 calculates a weighting coefficient for each pixel of the original image 41 and the smoothed image 42 using the state coefficient calculated in step S102.
  • the weighting coefficient related to the target pixel s is set to ⁇ s .
  • the calculation process of ⁇ s follows any of the following formula (5), formula (6), or formula (7).
  • weighting coefficient values of the corresponding pixels are calculated using the state coefficients of both the original image 41 and the smoothed image 42 as shown in the equations (5) to (7). It is a point to do.
  • Equation (5) and Equation (7) are obtained when the state coefficient is a value close to 0 in at least one of the original image 41 and the smoothed image 42, and it is determined that the target pixel s corresponds to the edge portion of the structure.
  • the weighting coefficient takes a value close to zero.
  • both the original image 41 and the smoothed image 42 are weighted when it is determined that the state coefficient has a value close to 1 and the target pixel s corresponds to a flat region (for example, a region considered to be the same organ). The coefficient is close to 1.
  • Equation (6) indicates that the weighting coefficient is close to 0 when the state coefficients of both the original image 41 and the smoothed image 42 are close to 0. Further, in the equation (6), when the state coefficients of both the original image 41 and the smoothed image 42 are close to 1, the weighting coefficient is close to 1. On the other hand, in the equation (6), when the difference between the values of the state coefficients of the original image 41 and the smoothed image 42 is large, the average value of both is the weighting coefficient.
  • Equation (5) or Equation (7) when it is desired to sharply restore the edge of the structure in the target image. Further, when it is desired to smoothly restore the edge of the structure in the target image, it is desirable to use Equation (6). Both are used properly depending on the application.
  • Equation (8) By performing weighted addition using Equation (8), the above-mentioned ⁇ (Constraint 1)
  • the pixel value of the target image is within the range of the pixel value of the original image and the pixel value of the smoothed image related to the corresponding pixel. To meet.
  • FIG. 7 shows a schematic diagram of an original image 41, a smoothed image 42, and a target image 43.
  • streak-like artifacts 51 are generated in the original image 41.
  • the smoothed image 42 the edge blur 52 of the structure occurs in the smoothed image 42.
  • a target image 43 in which the edge of the structure is maintained and the streak-like artifact 51 is removed is generated from the original image 41 and the smoothed image 42. can do.
  • the edge of the structure is maintained and the streak-like artifact 51 is removed.
  • a small region in which the reduction rate of the variation value of the smoothed image 42 with respect to the original image 41 is maximized is extracted as the feature amount calculation region, and the nonlinear function is calculated based on the feature amount of the feature amount calculation region.
  • a parameter ⁇ for adjusting the bending position is determined.
  • the relationship between the variation value reduction rate and the variation value of the original image 41 or the smoothed image 42 is analyzed in more detail, and a feature amount calculation region is extracted. This improves the accuracy of extracting a flat region including the streak-like artifact 51 as a feature amount calculation region.
  • the computing device 5 is the upper m of the standard deviation value reduction rate ⁇ i in a small area included in a predetermined range (for example, in the same slice). Extract regions. That is, the arithmetic unit 5 extracts from the small area where the standard deviation value reduction rate ⁇ i is the largest to the m-th order small area.
  • m is an arbitrary constant and is determined empirically.
  • the arithmetic unit 5 extracts, for each of the original image 41 and the smoothed image 42, a small region having the maximum standard deviation value among the upper m regions as a feature amount calculation region. Subsequent processing, such as feature amount calculation processing, nonlinear function shape determination processing, state coefficient calculation processing, weighting coefficient calculation processing, weighted addition processing, and the like are the same as those in the first embodiment.
  • the accuracy of extracting a flat region including streak-like artifacts 51 as a feature amount calculation region is improved.

Landscapes

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

Abstract

　構造物のエッジが保たれ、かつストリーク状のアーチファクトが除去される目的画像を生成することが可能な画像処理装置などを提供するために、演算装置は、原画像と平滑化画像の特徴量に基づき、非線形関数の形状を決定する(S101)。次に、演算装置は、S101において形状が決定される非線形関数によって、原画像及び平滑化画像の状態係数を算出する。(S102)。次に、演算装置は、S102において算出される状態係数を用いて、原画像及び平滑化画像の各画素について加重係数を算出する(S103)。次に、演算装置は、原画像と平滑化画像の加重加算を行い、目的画像を生成する(S104)。

Description

画像処理装置及び画像処理方法
　本発明は、X線CT装置などによって得られる画像に画像処理を施す画像処理装置に係り、構造物のエッジを保って、ストリーク状のアーチファクトを除去する技術に関する。
　X線CT装置は、被検体の周囲からX線を照射し、複数の投影角度において取得される実投影データを逆投影することにより、被検体の断層像を再構成する。
　前述のようなX線CT装置などによって得られる原情報に基づいて、肩等の扁平な部位を再構成する際、X線の照射線量が不十分な場合、X線減弱の大きい方向にストリーク状のアーチファクトが生じる。非特許文献1や非特許文献2では、投影データを平滑化した上で前述の画像再構成法を適用することによって、ストリーク状のアーチファクトを低減する手法が提案されている。以降では、これらのストリーク状のアーチファクトを低減する方法を「投影データ平滑化フィルタ」と呼ぶ。
　投影データ平滑化フィルタによって得られる画像は、ストリーク状のアーチファクトの低減効果に付随して比較的高コントラストな構造物のエッジにおいて暈けを伴う。ストリーク状のアーチファクトの低減効果と構造物のエッジ保存効果は、トレードオフの関係にあり、前者の効果を大きくするほど後者の効果が小さくなる。従って、強いストリーク状のアーチファクトを低減させる場合には、構造物のエッジに大きな暈けを生じてしまう。構造物のエッジの暈けは臨床診断上好ましくないため、構造物のエッジが保たれ、かつストリーク状のアーチファクトが除去される画像(以降「目的画像」と呼ぶ。)が望まれる。
　この問題の解決策として、投影データから前述の画像再構成法によって得られる画像(以降、「原画像」と呼ぶ。)と、投影データ平滑化フィルタを適用した投影データから前述の画像再構成法によって得られる画像(以降、「平滑化画像」と呼ぶ。)を適切に加重加算する方法が挙げられる。すなわち、画像上でストリーク状のアーチファクトを生ずる部分の画素では平滑化画像の画素値、構造物のエッジ部分では原画像の画素値の比率がそれぞれ大きくなるように加重加算することで、前述の目的画像が得られる。このような原情報を同じくする二枚の画像を混合する手法は、これまでにいくつか提案されている。
　特許文献1に記載の手法では、原画像から原画像よりも解像度が低い低解像度画像を作成し、低解像度画像からエッジを検出して低解像度なエッジ画像を作成し、低解像度なエッジ画像のエッジ輝度値から非線形関数によって加重係数を算出する。このとき、エッジの検出には注目画素を含む所定の領域の画素の値の最大値と最小値の差や分散値、濃度勾配などを使用する。原画像に対し、原画像から低解像度画像を差分した画像を前述の加重係数で加重加算することで、構造物のエッジを強調する。このような手法によって、原画像と平滑化画像の加重加算を行うことができる。
　特許文献2に記載の手法では、画像を小領域に分割し、小領域の標準偏差値を算出し、全小領域を対象として標準偏差値を階級とするヒストグラムを作成し、ヒストグラムから判別階級値を算出して各小領域を平坦領域と非平坦領域に分類する。画像上の構造物のエッジを検出することは、平坦領域と非平坦領域を識別することに等しい。両者が分類できれば、非平坦領域では原画像の画素値を使用し、平坦領域では平滑化画像の画素値を使用することで、原画像と平滑化画像の加重加算を行うことができる。
特許第3700798号公報 特開平07-93543号公報
J. Hsieh., "Adaptive streak artifact reduction in computed tomography resulting from excessive x-ray photon noise," Med. Phys., Vol.25,No.11, pp.2139-2147, 1998 T. Li et. al., "Nonlinear Sinogram Smoothing for Low-Dose X-Ray CT," IEEE. Trans. Nucl. Sci., Vol.51, No.5,pp.2505-2513, 2004
　ところで、強いストリーク状のアーチファクトは、画像に生ずる粒状性のノイズと異なり、直線状の構造物に近い特性を持つ。そうすると、特許文献1に記載の手法のように、注目画素を含む所定の領域の画素の値の最大値と最小値の差や分散値、濃度勾配などを使用してエッジの検出を行う場合、強いストリーク状のアーチファクトがエッジとして検出される可能性が高い。
　また、前述のように、強いストリーク状のアーチファクトは、構造物に近い特性を持つ。そうすると、特許文献2に記載の手法のように、標準偏差値を基に構造物の判定を行う場合、ストリーク状のアーチファクトを含む領域は非平坦領域として認識されるおそれがある。
　ストリーク状のアーチファクトを構造物や構造物のエッジとして誤認識した場合、ストリーク状のアーチファクトの一部ないし全部が、原画像と平滑化画像の加重加算後の画像に復元してしまう。つまり、特許文献1や特許文献2に記載の手法では、目的画像を得ることができない。
　目的画像を得ることができない本質的な理由は、特許文献1や特許文献2に記載の手法では、一様領域における粒状性のノイズと構造物を分類することを目的としており、ストリーク状のアーチファクトと構造物を分類するための指標を用いていないためである。
　本発明は、前述した問題点に鑑みてなされたものであり、その目的とすることは、構造物のエッジが保たれ、かつストリーク状のアーチファクトが除去される目的画像を生成することが可能な画像処理装置などを提供することである。
　前述した目的を達成するための第1の発明は、原情報から生成される原画像、及び同一の前記原情報から生成され、少なくともストリーク状のアーチファクトが低減されている平滑化画像の対応する各画素を加重係数によって加重加算して、構造物のエッジが保たれ、かつストリーク状のアーチファクトが除去される目的画像を生成する画像処理装置であって、前記原画像及び前記平滑化画像の特徴量に基づき、非線形関数の形状を決定する関数形状決定部と、前記非線形関数に基づいて前記加重係数を算出する加重係数算出部と、を具備することを特徴とする画像処理装置である。
　第2の発明は、原情報から生成される原画像、及び同一の前記原情報から生成され、少なくともストリーク状のアーチファクトが低減されている平滑化画像の対応する各画素を加重係数によって加重加算して、構造物のエッジが保たれ、かつストリーク状のアーチファクトが除去される目的画像を生成する画像処理方法であって、前記原画像及び前記平滑化画像の特徴量に基づき、非線形関数の形状を決定する関数形状決定ステップと、前記非線形関数に基づいて前記加重係数を算出する加重係数算出ステップと、を含むことを特徴とする画像処理方法である。
　本発明により、構造物のエッジが保たれ、かつストリーク状のアーチファクトが除去される目的画像を生成することが可能な画像処理装置などを提供することができる。
X線CT装置の全体概観図 X線CT装置の全体構成図 非線形関数の関数形状の例を示す図 全体の処理フローを示す図 関数形状決定ステップに係る処理フローを示す図 原画像及び平滑化画像の領域分割の例を示す図 原画像、平滑化画像、目的画像の例を示す図
　以下図面に基づいて、本発明の実施形態を詳細に説明する。最初に、図1、図2を参照しながら、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(X線検出器)、コリメータ13、駆動装置14、中央制御装置15、X線制御装置16、高電圧発生装置17、スキャナ制御装置18、寝台制御装置19、寝台移動計測装置20、コリメータ制御装置21、プリアンプ22、A/Dコンバータ23などから構成されている。
　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線の照射が行われる。
　中央制御装置15は、操作ユニット3における入力装置6から撮影条件や再構成パラメータを入力し、撮影に必要な制御信号を、コリメータ制御装置21、X線制御装置16、スキャナ制御装置18、寝台制御装置19に送信する。コリメータ制御装置21は、制御信号に基づいてコリメータ13の位置を制御する。
　撮影スタート信号を受けて撮影が開始されると、X線制御装置16は、制御信号に基づいて高電圧発生装置17を制御する。高電圧発生装置17は、X線管11に管電圧、管電流を印加する。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などによって取得される原情報から生成される原画像、及び同一の原情報から生成され、少なくともストリーク状のアーチファクトが低減されている平滑化画像の対応する各画素を加重係数によって加重加算して、構造物のエッジが保たれ、かつストリーク状のアーチファクトが除去される目的画像を生成する画像処理方法について説明する。
　本発明の画像処理方法は、以下の3つの制約条件を満たすように、目的画像を生成する。
　(制約条件1)目的画像の画素値は、対応する画素に係る原画像の画素値及び平滑化画像の画素値の範囲内とする。
　(制約条件2)目的画像は、空間的に(つまり1枚の画像空間において)、滑らかな画像とする。
　(制約条件3)目的画像は、平滑化画像と比較して、同一の臓器と思われる領域(＝平坦領域の1つ)において、ストリーク状のアーチファクトなどのノイズが増えないようにする。
　制約条件1は、X線CT装置1などによって取得される原情報(但し、ノイズを除く。
)を失わないための条件である。制約条件1を言い換えると、min(原画像の画素値、平滑化画像の画素値)≦目的画像の画素値≦max(原画像の画素値、平滑化画像の画素値)、となる。なお、min()は最小値を出力する演算子であり、max()は最大値を出力する演算子である。
　制約条件2は、目的画像で生じうる不連続性を回避するための条件である。制約条件3は、主に粒状性のノイズとストリーク状のアーチファクトを除去するための条件である。
本発明の画像処理方法は、特に、制約条件3を満たすことが特徴である。
　加重加算に用いる加重係数は、原画像と平滑化画像の状態を表す状態係数によって決定する。状態係数は、原画像もしくは平滑化画像の注目画素とその近接画素の画素値差を変数とする非線形関数の値とする。
　本発明の実施の形態では、非線形関数として一般化ガウス関数を使用する。但し、非線形関数は、滑らかな連続関数であって、複数の任意のパラメータを持ち、その形状が任意に調整可能であれば良く、特に限定する必要はない。非線形関数が、滑らかな連続関数であれば、前述の「(制約条件2)目的画像は、空間的に(つまり1枚の画像空間において)、滑らかな画像とする。」を満たす。非線形関数の他の例としては、ロジスティック関数などが挙げられる。
　画像ベクトルをx＝｛x1，・・・，xJ｝、注目画素sに対する近接画素rの集合をNsとおいて、一般化ガウス関数により、注目画素sの状態係数gs(x)を次式で定義する。
Figure JPOXMLDOC01-appb-M000001
　ここで、pは、一般化ガウス関数の傾きを調整する任意のパラメータであり、全てのスライスに対して同じ値とする。νは、一般化ガウス関数の折れ曲がり位置を調整する任意のパラメータである。wsrは、注目画素sと近接画素rの間の距離に応じた重み係数であり、例えば次式のように定義する。
Figure JPOXMLDOC01-appb-M000002
　ここで、lsrは、注目画素sと近接画素rの間の距離である。Lは、検出器の素子の大きさに対する画素の大きさの比である。但し、wsrは、式(2)に限定する必要はなく、近接画素rが注目画素sに近いほど、値が大きくなる関数であれば良い。また、注目画素sの近傍画素の集合Ns、及び一般化ガウス関数の傾きpは、経験的に決定するものとする。
　式(1)において、注目画素sと近接画素rの画素値が等しければ、状態係数は1となり、注目画素sと近接画素rとの画素値の差が大きいほど、状態係数は0に近づく。注目画素sと近接画素rとの画素値差に対し、νを任意の定数cおよび2c、4cと設定したときの関数形状の変化を図3に例示する。画素値差を固定して考えると、νを大きくする程、状態係数が1に近い値になり易くなる。本発明では、1枚の断層像内において、ストリーク状のアーチファクトによって生じると思われる注目画素sと近接画素rとの画素値差に対して、状態係数がほぼ1に近い値になるように、関数の折れ曲がり位置を調整するためのパラメータνを決定する。
　［第1の実施形態］
　以降では、図4に示す処理フローの流れに沿って、適宜他の図を参照しながら、本発明の第1の実施形態について説明する。
　＜ステップS101＞
　X線CT装置1の演算装置5(以降、「演算装置5」と表記する。)は、原画像と平滑化画像の特徴量に基づき、非線形関数の形状を決定する。非線形関数の形状を決定する処理については、図5を参照しながら説明する。
　＜ステップS201＞
　演算装置5は、原画像及び平滑化画像の領域を、互いに対応する小領域に分割する。例えば、図6に示すように、演算装置5は、原画像41及び平滑化画像42の領域を格子状に区切り、小領域に分割する。小領域の大きさは、経験的に決める。
　但し、小領域の分割方法は、この例に限定されない。本発明では、原画像41及び平滑化画像42の領域が、同一の分割方法によって小領域に分割されれば良い。また、小領域の形状も、矩形に限定されない。本発明では、近接する複数の画素が同一の小領域に含まれていれば良い。また、一の画素が複数の小領域に重複して含まれていても良い。
　＜ステップS202、S203＞
　演算装置5は、原画像41及び平滑化画像42の両方に係る各小領域に対して、小領域に含まれる画素の画素値から、原画像41及び平滑化画像42のばらつき値を算出する。
ばらつき値は、例えば、標準偏差値や、(最大値－最小値)などである。ばらつき値は、各小領域のばらつきを示す値であって、各小領域に含まれる画素の画素値から算出される統計量であれば良い。以降では、混乱を避ける為に、標準偏差値を例にして説明する。
　ここで、説明のため、小領域に通し番号i＝1，・・・，Ｉを振り、i番目の小領域に注目する。i番目の小領域における原画像41の標準偏差値をσi (org)、平滑化画像42の標準偏差値をσi (smt)とおく。
　＜ステップS204＞
　演算装置5は、小領域ごとに、原画像41のばらつき値を基準とする平滑化画像42のばらつき値の低減率を算出する。本実施形態では、標準偏差値の低減率を算出する。i番目の小領域について、原画像41の標準偏差値を基準とする平滑化画像42の標準偏差値の低減率をρiとおき、次式で算出する。
Figure JPOXMLDOC01-appb-M000003
　演算装置5は、式(3)に基づき、全ての小領域において、標準偏差値の低減率ρiを算出する。式(3)のように、原画像41及び平滑化画像42の両方の画素値から算出される特徴量(例えば、低減率ρi)を用いて、非線形関数の形状を決定し、図4に示す処理を行うことによって、前述の「(制約条件3)目的画像は、平滑化画像と比較して、同一の臓器と思われる領域において、ストリーク状のアーチファクトなどのノイズが増えないようにする。」を満たす。
　尚、原情報を同じくする二枚の画像を混合する手法に関する従来技術(例えば、特許文献1や2)では、一方の画像のみに含まれる画素の画素値から算出される特徴量しか用いていない。例えば、一方の画像のみに含まれる画素の画素値から算出される標準偏差値では、任意の平滑化処理によって、ノイズが減ったのか、それとも構造物がぼけたのか、の区別ができない。従って、従来技術では、本発明のようにストリーク状のアーチファクトが除去される目的画像を生成することができない。
　＜ステップS205＞
　演算装置5は、低減率ρiに基づいて、小領域の集合から特徴量算出領域を抽出する。
第1の実施の形態では、演算装置5は、所定の範囲内(例えば、同一スライス内)に含まれる小領域の中で、低減率が最大の小領域を、特徴量算出領域として抽出する。尚、所定の範囲は、複数のスライス内としても良いし、全てのスライス内としても良い。
　以降では、k番目の小領域を、特徴量算出領域として抽出したものとして説明する。
　＜ステップS206＞
　演算装置5は、特徴量算出領域(k番目の小領域)の画素値から算出される特徴量から、非線形関数の形状を決定する。以下では、ステップS202において算出される原画像41の標準偏差値σk (org)及びステップS203において算出される平滑化画像42の標準偏差値σk (smt)を用いて、特徴量を算出する例を説明する。但し、本発明は、この例に限定されるわけではなく、標準偏差値に代えて、他のばらつき値(例えば、小領域内の画素の最大値および最小値の差)を用いても良い。
　任意の定数をα(0＜α＜1)とおき、画素値差がσk(σk (org)とσk (smt)の総称記号)のとき、状態係数がαの値をとるように、非線形関数の形状を設定する(図3参照)。αは、例えば、0.99、0.98、0.97などの実数値である。例えば、式(1)の一般化ガウス関数であれば、関数の折れ曲がり位置を調整するためのパラメータνは、次式のように定義する。
Figure JPOXMLDOC01-appb-M000004
　ここで、tは1ではない正の係数であり、既に撮影済の画像から得られる統計的な情報から、ストリーク状のアーチファクトなどのノイズの影響を考慮して、経験的に決定する。例えば、ガウシアンノイズ(粒状性のノイズ)しか発生しない画像であれば、t＝1とし、式(4)にtを導入する必要はない。しかし、本発明は、ガウシアンノイズよりも強いノイズであるストリーク状のアーチファクトが発生する画像を対象とするので、式(4)に1より大きい値のtを導入する。
　式(4)において、原画像41及び平滑化画像42のそれぞれにおいて算出される非線形関数の折れ曲がり位置を調整するためのパラメータを、νk (org)およびνk (smt)とおく。
　演算装置5は、図5に示す処理によって、非線形関数の形状を、所定の範囲毎(例えば、スライス毎)に動的に決定する。第1の実施形態における特徴の一つは、原画像41に対する平滑化画像42の標準偏差値の低減率ρiが最大となる領域を特徴量算出領域として抽出し、特徴量算出領域における標準偏差値を用いて非線形関数の形状を決定する点にある。
　前述の投影データ平滑化フィルタは、ストリーク状のアーチファクトを大幅に低減するという特性を持つ。この特性から、ストリーク状のアーチファクトを含む小領域において、標準偏差値の低減率ρiが大きくなる。また、投影データ平滑化フィルタは、構造物のぼけを伴うという特性を持つ。この特性から、構造物を含む小領域では原画像の標準偏差値が大きくなり、式(3)の分母が大きくなるので、結果として、標準偏差値の低減率ρiが小さくなる。つまり、標準偏差値の低減率ρiが最大となる領域では、ストリーク状のアーチファクトを含む可能性が高い。
　従って、標準偏差値の低減率ρiが最大となる領域を特徴量算出領域として抽出すると、ストリーク状のアーチファクトを含む小領域を特徴量算出領域として抽出することになる。そして、ストリーク状のアーチファクトを含む小領域における標準偏差値を式(4)に代入し、非線形関数の折れ曲がり位置を調整するためのパラメータνの値を決定することによって、ひいては、ストリーク状のアーチファクトが除去される目的画像を得ることができる。
　図4の説明に戻る。以降では、非線形関数によって、原画像41と平滑化画像42のそれぞれの状態を表す状態係数を算出し、原画像41と平滑化画像42の状態係数から加重加算の加重係数を算出する処理について説明する。
　＜ステップS102＞
　演算装置5は、ステップS101において形状が決定される非線形関数によって、原画像41及び平滑化画像42の状態係数を算出する。状態係数の算出処理は、式(1)に従う。
　＜ステップS103＞
　演算装置5は、ステップS102において算出される状態係数を用いて、原画像41及び平滑化画像42の各画素について加重係数を算出する。原画像41の画像ベクトルをx(org)＝｛x1 (org)，・・・，xJ (org)｝とおき、平滑化画像42の画像ベクトルをx(smt)＝｛x1 (smt)，・・・，xJ (smt)｝とおく。画像の加重加算において、注目画素sに係る加重係数をλsとおく。λsの算出処理は、以下に示す式(5)、式(6)又は式(7)のいずれかに従う。
Figure JPOXMLDOC01-appb-M000005
Figure JPOXMLDOC01-appb-M000006
Figure JPOXMLDOC01-appb-M000007
　本発明の特徴の一つは、式(5)～式(7)に示すように、原画像41と平滑化画像42の両方の状態係数を用いて、対応する画素の加重係数の値を算出する点である。
　式(5)及び式(7)は、原画像41と平滑化画像42の少なくともいずれかにおいて、状態係数が0に近い値をとり、注目画素sが構造物のエッジ部分にあたると判断される場合に、加重係数が0に近い値をとる。一方、原画像41と平滑化画像42のいずれも、状態係数が1に近い値をとり、注目画素sが平坦領域(例えば、同一の臓器と思われる領域)にあたると判断される場合に、加重係数が1に近い値をとる。
　式(6)は、原画像41と平滑化画像42の両方の状態係数が0に近い値をとる場合に、加重係数が0に近い値をとる。また、式(6)は、原画像41と平滑化画像42の両方の状態係数が1に近い値をとる場合に、加重係数が1に近い値をとる。一方、式(6)は、原画像41と平滑化画像42の状態係数の値の差が大きい場合には、両者の平均値が加重係数となる。
　以上より、目的画像において、構造物のエッジを鋭く回復させたい場合には、式(5)又は式(7)を用いることが望ましい。また、目的画像において、構造物のエッジを滑らかに回復させたい場合には、式(6)を用いることが望ましい。両者は、用途によって使い分けるものとする。
　＜ステップS104＞
　演算装置5は、ステップS103において算出される加重係数を用いて、原画像41と平滑化画像42の加重加算を行い、目的画像を生成する。加重加算後の目的画像の画像ベクトルをx(mrg)＝｛x1 (mrg)，・・・，xJ (mrg)｝とする。演算装置5は、目的画像の注目画素sについて、次式を用いて加重加算を行う。
Figure JPOXMLDOC01-appb-M000008
　式(8)を用いて加重加算を行うことによって、前述の「(制約条件1)目的画像の画素値は、対応する画素に係る原画像の画素値及び平滑化画像の画素値の範囲内とする。」を満たす。
　図7には、原画像41、平滑化画像42、目的画像43の模式図が図示されている。図7に示すように、原画像41では、ストリーク状のアーチファクト51が発生している。
また、平滑化画像42では、構造物のエッジの暈け52が発生している。図4に示す本発明の画像処理方法によれば、このような原画像41及び平滑化画像42から、構造物のエッジが保たれ、かつストリーク状のアーチファクト51が除去される目的画像43を生成することができる。図7に示すように、目的画像43では、構造物のエッジが保たれており、かつストリーク状のアーチファクト51が除去されている。
　［第2の実施形態］
　以降では、本発明における第2の実施形態について説明する。尚、第1の実施形態と共通する内容については説明を省略する。
　第1の実施形態では、原画像41に対する平滑化画像42のばらつき値の低減率が最大となる小領域を特徴量算出領域として抽出し、特徴量算出領域の特徴量に基づいて、非線形関数の折れ曲がり位置を調整するためのパラメータνを決定している。このとき、ばらつき値の低減率と原画像41又は平滑化画像42のばらつき値の関係をより詳細に解析し、特徴量算出領域を抽出する。これによって、ストリーク状のアーチファクト51を含む平坦領域を特徴量算出領域として抽出する精度が向上する。
　第2の実施形態では、演算装置5は、図5のステップS205において、所定の範囲内(例えば、同一スライス内)に含まれる小領域の中で、標準偏差値の低減率ρiの上位m領域を抽出する。つまり、演算装置5は、標準偏差値の低減率ρiが最大の小領域からm番目の順位の小領域までを抽出する。ここで、mは、任意の定数であり、経験的に決定する。
　次に、演算装置5は、原画像41および平滑化画像42のそれぞれについて、上位m領域の中で、標準偏差値が最大となる小領域を、特徴量算出領域として抽出する。後続の処理である特徴量の算出処理、非線形関数の形状決定処理、状態係数の算出処理、加重係数の算出処理、加重加算処理などは、第1の実施形態と同様である。
　第2の実施の形態によれば、ストリーク状のアーチファクト51を含む平坦領域を特徴量算出領域として抽出する精度が向上する。
　本発明の様々な実施形態に関する以上の記述から、本発明の目的が達成されることは明らかである。本発明を詳細にわたって記述すると共に図示したが、これらは説明及び例示のみを意図したものであって、これらに限定されるものではない。また、本発明の要旨は、特許請求の範囲によってのみ限定されるものとする。
　1　X線CT装置、　2　スキャナ、　3　操作ユニット、　4　寝台、　5　演算装置、　6　入力装置、　7　表示装置、　8　記憶装置、　10　被検体、　11　X線管、　41　原画像、　42　平滑化画像、　43　目的画像、　51　ストリーク状のアーチファクト、　52　構造物のエッジの暈け

Claims (14)

  1. 　原情報から生成される原画像、及び同一の前記原情報から生成され、少なくともストリーク状のアーチファクトが低減されている平滑化画像の対応する各画素を加重係数によって加重加算して、構造物のエッジが保たれ、かつストリーク状のアーチファクトが除去される目的画像を生成する画像処理装置であって、
    　前記原画像及び前記平滑化画像の特徴量に基づき、非線形関数の形状を決定する関数形状決定部と、
    　前記非線形関数に基づいて前記加重係数を算出する加重係数算出部と、
    を具備することを特徴とする画像処理装置。
  2. 　前記関数形状決定部は、
    　前記原画像及び前記平滑化画像の領域を、互いに対応する小領域に分割する分割部と、
    　前記原画像及び前記平滑化画像の両方に係る前記小領域に対してばらつき値を算出し、前記小領域ごとに前記原画像のばらつき値を基準とする前記平滑化画像のばらつき値の低減率を算出する低減率算出部と、
    　前記低減率に基づいて、前記小領域の集合から特徴量算出領域を抽出する領域抽出部と、
    　前記特徴量算出領域の画素値から、前記特徴量を算出する特徴量算出部と、
    を含むことを特徴とする請求項1に記載の画像処理装置。
  3. 　前記領域抽出部は、所定の範囲において前記低減率が最大の前記小領域を、前記特徴量算出領域として抽出する
    ことを特徴とする請求項2に記載の画像処理装置。
  4. 　前記領域抽出部は、所定の範囲において前記低減率が最大の前記小領域から一定順位の前記小領域までの中で、前記ばらつき値が最大の前記小領域を、前記特徴量算出領域として抽出する
    　ことを特徴とする請求項2に記載の画像処理装置。
  5. 　前記非線形関数は、前記原画像又は前記平滑化画像の注目画素とその近接画素との画素値差を変数とする滑らかな連続関数であることを特徴とする請求項1に記載の画像処理装置。
  6. 　前記非線形関数は、前記画素値差が小さいほど1に近づき、前記画素値が大きいほど0に近づく関数であることを特徴とする請求項5に記載の画像処理装置。
  7. 　前記加重係数算出部は、前記非線形関数を用いて前記原画像及び前記平滑化画像のそれぞれの状態を表す状態係数を算出し、前記状態係数に基づいて前記加重係数を算出することを特徴とする請求項1に記載の画像処理装置。
  8. 　原情報から生成される原画像、及び同一の前記原情報から生成され、少なくともストリーク状のアーチファクトが低減されている平滑化画像の対応する各画素を加重係数によって加重加算して、構造物のエッジが保たれ、かつストリーク状のアーチファクトが除去される目的画像を生成する画像処理方法であって、
    　前記原画像及び前記平滑化画像の特徴量に基づき、非線形関数の形状を決定する関数形状決定ステップと、
    　前記非線形関数に基づいて前記加重係数を算出する加重係数算出ステップと、
    を含むことを特徴とする画像処理方法。
  9. 　前記関数形状決定ステップは、
    　前記原画像及び前記平滑化画像の領域を、互いに対応する小領域に分割する分割ステップと、
    　前記原画像及び前記平滑化画像の両方に係る前記小領域に対してばらつき値を算出し、前記小領域ごとに前記原画像のばらつき値を基準とする前記平滑化画像のばらつき値の低減率を算出する低減率算出ステップと、
    　前記低減率に基づいて、前記小領域の集合から特徴量算出領域を抽出する領域抽出ステップと、
    　前記特徴量算出領域の画素値から、前記特徴量を算出する特徴量算出ステップと、
    を含むことを特徴とする請求項8に記載の画像処理方法。
  10. 　前記領域抽出ステップは、所定の範囲において前記低減率が最大の前記小領域を、前記特徴量算出領域として抽出する
    ことを特徴とする請求項9に記載の画像処理方法。
  11. 　前記領域抽出ステップは、所定の範囲において前記低減率が最大の前記小領域から一定順位の前記小領域までの中で、前記ばらつき値が最大の前記小領域を、前記特徴量算出領域として抽出する
    　ことを特徴とする請求項9に記載の画像処理方法。
  12. 　前記非線形関数は、前記原画像又は前記平滑化画像の注目画素とその近接画素との画素値差を変数とする滑らかな連続関数であることを特徴とする請求項8に記載の画像処理方法。
  13. 　前記非線形関数は、前記画素値差が小さいほど1に近づき、前記画素値が大きいほど0に近づく関数であることを特徴とする請求項12に記載の画像処理方法。
  14. 　前記加重係数算出ステップは、前記非線形関数を用いて前記原画像及び前記平滑化画像のそれぞれの状態を表す状態係数を算出し、前記状態係数に基づいて前記加重係数を算出することを特徴とする請求項8に記載の画像処理方法。
PCT/JP2013/057118 2012-03-27 2013-03-14 画像処理装置及び画像処理方法 WO2013146283A1 (ja)

Priority Applications (4)

Application Number Priority Date Filing Date Title
US14/372,330 US9454804B2 (en) 2012-03-27 2013-03-14 Image processing device and image processing method
IN5834DEN2014 IN2014DN05834A (ja) 2012-03-27 2013-03-14
CN201380005062.5A CN104066378B (zh) 2012-03-27 2013-03-14 图像处理装置以及图像处理方法
JP2014507658A JP6124868B2 (ja) 2012-03-27 2013-03-14 画像処理装置及び画像処理方法

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2012-070823 2012-03-27
JP2012070823 2012-03-27

Publications (1)

Publication Number Publication Date
WO2013146283A1 true WO2013146283A1 (ja) 2013-10-03

Family

ID=49259547

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2013/057118 WO2013146283A1 (ja) 2012-03-27 2013-03-14 画像処理装置及び画像処理方法

Country Status (5)

Country Link
US (1) US9454804B2 (ja)
JP (1) JP6124868B2 (ja)
CN (1) CN104066378B (ja)
IN (1) IN2014DN05834A (ja)
WO (1) WO2013146283A1 (ja)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2015100539A (ja) * 2013-11-26 2015-06-04 日立アロカメディカル株式会社 超音波診断装置
US10217248B2 (en) 2014-08-27 2019-02-26 General Electric Company Method for removing streak from detector cell with performance difference
JP2020078377A (ja) * 2018-11-12 2020-05-28 株式会社日立製作所 画像再構成装置および画像再構成方法
JP2020156582A (ja) * 2019-03-25 2020-10-01 キヤノンメディカルシステムズ株式会社 超音波診断装置、医用画像処理装置および医用画像処理プログラム

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
DE102013221603A1 (de) * 2013-10-24 2015-04-30 Siemens Aktiengesellschaft CT-System mit Recheneinheit und Verfahren zur Rekonstruktion und Befundung von CT-Bilddarstellungen
JP6778158B2 (ja) 2017-07-27 2020-10-28 株式会社日立製作所 Ｘ線ｃｔ装置、画像生成方法および画像生成プログラム
JP6987352B2 (ja) * 2017-11-17 2021-12-22 富士フイルムヘルスケア株式会社 医用画像処理装置および医用画像処理方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2001198120A (ja) * 1999-12-08 2001-07-24 Koninkl Philips Electronics Nv 再構成画像を組合わせる方法
JP3700798B2 (ja) * 1996-08-19 2005-09-28 富士写真フイルム株式会社 画像処理方法および装置
JP2007044275A (ja) * 2005-08-10 2007-02-22 Hitachi Medical Corp マルチエナジーｘ線ｃｔ装置
JP2010226694A (ja) * 2009-02-24 2010-10-07 Ricoh Co Ltd 画像処理装置及び画像処理方法

Family Cites Families (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP1156451B1 (en) * 1995-09-29 2004-06-02 Fuji Photo Film Co., Ltd. Image processing method and apparatus
US6667815B1 (en) * 1998-09-30 2003-12-23 Fuji Photo Film Co., Ltd. Method and apparatus for processing images
US7409092B2 (en) * 2002-06-20 2008-08-05 Hrl Laboratories, Llc Method and apparatus for the surveillance of objects in images
JP4413504B2 (ja) 2003-02-13 2010-02-10 株式会社東芝 医用画像処理装置、医用画像処理方法および医用画像処理プログラム
JP4069943B2 (ja) * 2003-12-03 2008-04-02 株式会社ニコン ノイズ除去の強弱を画面内でコントロールする画像処理装置、画像処理プログラム、画像処理方法、および電子カメラ
JP4780374B2 (ja) * 2005-04-21 2011-09-28 Ｎｋワークス株式会社 粒状ノイズ抑制のための画像処理方法及びプログラム及びこの方法を実施する粒状抑制処理モジュール
EP1882449A4 (en) 2005-05-18 2010-05-05 Hitachi Medical Corp RADIOGRAPHY DEVICE AND IMAGE PROCESSING PROGRAM
JP5042465B2 (ja) * 2005-05-18 2012-10-03 株式会社日立メディコ 放射線撮影装置、画像処理方法
JP5253835B2 (ja) * 2008-02-19 2013-07-31 株式会社キーエンス 画像生成装置、画像生成方法及びコンピュータプログラム
KR101481551B1 (ko) * 2008-06-03 2015-01-13 엘지전자 주식회사 영상 노이즈 제거 장치 및 방법
JP5143038B2 (ja) * 2009-02-02 2013-02-13 オリンパス株式会社 画像処理装置及び画像処理方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP3700798B2 (ja) * 1996-08-19 2005-09-28 富士写真フイルム株式会社 画像処理方法および装置
JP2001198120A (ja) * 1999-12-08 2001-07-24 Koninkl Philips Electronics Nv 再構成画像を組合わせる方法
JP2007044275A (ja) * 2005-08-10 2007-02-22 Hitachi Medical Corp マルチエナジーｘ線ｃｔ装置
JP2010226694A (ja) * 2009-02-24 2010-10-07 Ricoh Co Ltd 画像処理装置及び画像処理方法

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2015100539A (ja) * 2013-11-26 2015-06-04 日立アロカメディカル株式会社 超音波診断装置
US10217248B2 (en) 2014-08-27 2019-02-26 General Electric Company Method for removing streak from detector cell with performance difference
JP2020078377A (ja) * 2018-11-12 2020-05-28 株式会社日立製作所 画像再構成装置および画像再構成方法
JP7077208B2 (ja) 2018-11-12 2022-05-30 富士フイルムヘルスケア株式会社 画像再構成装置および画像再構成方法
JP2020156582A (ja) * 2019-03-25 2020-10-01 キヤノンメディカルシステムズ株式会社 超音波診断装置、医用画像処理装置および医用画像処理プログラム
JP7297485B2 (ja) 2019-03-25 2023-06-26 キヤノンメディカルシステムズ株式会社 超音波診断装置、医用画像処理装置および医用画像処理プログラム

Also Published As

Publication number Publication date
JPWO2013146283A1 (ja) 2015-12-10
JP6124868B2 (ja) 2017-05-10
US9454804B2 (en) 2016-09-27
CN104066378B (zh) 2016-10-12
US20150010224A1 (en) 2015-01-08
IN2014DN05834A (ja) 2015-05-15
CN104066378A (zh) 2014-09-24

Similar Documents

Publication Publication Date Title
JP6124868B2 (ja) 画像処理装置及び画像処理方法
Liu et al. Total variation-stokes strategy for sparse-view X-ray CT image reconstruction
JP7139119B2 (ja) 医用画像生成装置及び医用画像生成方法
JP5978429B2 (ja) 医用画像処理装置、医用画像処理方法
US9265475B2 (en) Methods and apparatus for scatter correction for CBCT system and cone-beam image reconstruction
US7379575B2 (en) Method for post- reconstructive correction of images of a computer tomograph
EP3195265B1 (en) Iterative image reconstruction with a sharpness driven regularization parameter
US10395397B2 (en) Metal artifacts reduction for cone beam CT
US8965078B2 (en) Projection-space denoising with bilateral filtering in computed tomography
US8885903B2 (en) Method and apparatus for statistical iterative reconstruction
US10789738B2 (en) 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
JP2017196404A (ja) 画像処理装置、ｘ線ｃｔ装置及び画像処理方法
JP2019516460A (ja) 空間とスペクトル情報に基づく複数エネルギーのｃｔ画像におけるノイズ制御のためのシステムと方法
Zeng et al. Penalized weighted least-squares approach for multienergy computed tomography image reconstruction via structure tensor total variation regularization
JP2019130302A (ja) 医用画像処理装置及びｘ線ｃｔシステム
WO2014167935A1 (ja) Ｘ線ｃｔ装置、再構成演算装置、及び再構成演算方法
Van Aarle et al. Super-resolution for computed tomography based on discrete tomography
US10813616B2 (en) Variance reduction for monte carlo-based scatter estimation
JP2011203160A (ja) Ｘ線ｃｔ画像再構成方法及びｘ線ｃｔ画像再構成プログラム
CN114387359A (zh) 一种三维x射线低剂量成像方法及装置
CN112842370A (zh) 用于x射线成像中的参数噪声调制的方法和系统
CN107106114B (zh) 运算装置、x射线ct装置及图像重构方法
Xie et al. Scatter correction for cone-beam computed tomography using self-adaptive scatter kernel superposition
JP2013119021A (ja) Ｘ線ｃｔ装置及び画像処理方法
US20050018889A1 (en) Systems and methods for filtering images

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

Country of ref document: EP

Kind code of ref document: A1

ENP Entry into the national phase

Ref document number: 2014507658

Country of ref document: JP

Kind code of ref document: A

WWE Wipo information: entry into national phase

Ref document number: 14372330

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

Country of ref document: EP

Kind code of ref document: A1