WO2012143979A1 - Ｘ線診断装置およびｘ線診断用プログラム - Google Patents

Ｘ線診断装置およびｘ線診断用プログラム Download PDF

Info

Publication number
WO2012143979A1
WO2012143979A1 PCT/JP2011/002377 JP2011002377W WO2012143979A1 WO 2012143979 A1 WO2012143979 A1 WO 2012143979A1 JP 2011002377 W JP2011002377 W JP 2011002377W WO 2012143979 A1 WO2012143979 A1 WO 2012143979A1
Authority
WO
WIPO (PCT)
Prior art keywords
harmonic
image
ray
diagnostic apparatus
unit
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Ceased
Application number
PCT/JP2011/002377
Other languages
English (en)
French (fr)
Inventor
直紀 長谷川
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Shimadzu Corp
Original Assignee
Shimadzu Corp
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Shimadzu Corp filed Critical Shimadzu Corp
Priority to EP11864008.5A priority Critical patent/EP2700359B1/en
Priority to PCT/JP2011/002377 priority patent/WO2012143979A1/ja
Priority to JP2013510742A priority patent/JP5527481B2/ja
Priority to US14/113,176 priority patent/US9232925B2/en
Priority to CN201180070344.4A priority patent/CN103491875B/zh
Publication of WO2012143979A1 publication Critical patent/WO2012143979A1/ja
Anticipated expiration legal-status Critical
Ceased legal-status Critical Current

Links

Images

Classifications

    • 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/42—Arrangements for detecting radiation specially adapted for radiation diagnosis
    • A61B6/4291—Arrangements for detecting radiation specially adapted for radiation diagnosis the detector being combined with a grid or grating
    • 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 OR CALCULATING; COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/10—Image enhancement or restoration using non-spatial domain filtering
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/20—Image enhancement or restoration using local operators
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/70—Denoising; Smoothing
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; 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/52—Devices using data or image processing specially adapted for radiation diagnosis
    • A61B6/5205—Devices using data or image processing specially adapted for radiation diagnosis involving processing of raw data to produce diagnostic data
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/52—Devices using data or image processing specially adapted for radiation diagnosis
    • A61B6/5258—Devices using data or image processing specially adapted for radiation diagnosis involving detection or reduction of artifacts or noise
    • A61B6/5282—Devices using data or image processing specially adapted for radiation diagnosis involving detection or reduction of artifacts or noise due to scatter
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; 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/10116—X-ray image
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; 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/20004—Adaptive image processing
    • G06T2207/20012—Locally adaptive
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; 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/20048—Transform domain processing
    • G06T2207/20056—Discrete and fast Fourier transform, [DFT, FFT]
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; 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

Definitions

  • the present invention relates to an X-ray diagnosis apparatus and an X-ray diagnosis program for removing a moire pattern of a grid reflected in an X-ray image.
  • a conventional X-ray diagnostic apparatus includes an X-ray tube 103 that irradiates a subject M with X-rays, and an X-ray detector (for example, a flat panel) that detects X-rays transmitted through the subject M.
  • FPD Type X-ray detector
  • X-ray detector 104 scattered ray removal grid
  • the grid 105 is configured by alternately arranging an absorber (for example, lead) that absorbs X-rays and a transmissive body (for example, aluminum or air) that transmits X-rays. X-rays scattered and incident from an oblique direction are absorbed by the absorber by the grid 105, and only the X-rays that have passed through the transmissive body are detected by the X-ray detector 104. Thereby, a clear image can be acquired.
  • an absorber for example, lead
  • a transmissive body for example, aluminum or air
  • the conventional X-ray diagnostic apparatus 101 has a problem in that the moire pattern of the grid 105 is reflected in the acquired image due to the difference between the resolution of the X-ray detector 104 and the density of the grid 105.
  • Various proposals have been made for countermeasures against the moire pattern of the grid 105 (see, for example, Patent Document 1).
  • the conventional X-ray diagnostic apparatus 101 shown in FIG. 11 includes a moire pattern removal processing unit 121 and removes the first and second harmonics of the moire pattern of the grid 105 reflected in the image. .
  • a moire pattern removal processing unit 121 removes the first and second harmonics of the moire pattern of the grid 105 reflected in the image.
  • one-dimensional Fourier transform is performed on an image in which a moire pattern is reflected, and a peak frequency indicating the first harmonic of the moire pattern is detected from the Fourier transform result. Since the second harmonic is twice (integer multiple) the first harmonic, it is calculated from the peak frequency of the first harmonic.
  • a filter for extracting the first harmonic and the second harmonic of the moire pattern is created, and using this filter, the first harmonic and Extract the second harmonic. Then, the first harmonic and the second harmonic of the moire pattern are removed by subtracting the first harmonic and the second harmonic from the image in which the moire pattern is reflected.
  • the Nyquist frequency Ny indicates the limit of the spatial frequency that can be sampled at the pixel pitch ⁇ (mm).
  • the resolution of the X-ray detector 104 is hereinafter abbreviated as “3.33 lp / mm” as appropriate.
  • the density of the grid 105 is, for example, 50 lines / cm. This 50 line / cm is also expressed as 5 line / mm or 5 lp / mm. “Line / cm”, “lp / mm”, etc. indicate the number of absorbers per unit length (1 cm or 1 mm) or the number of absorbers / transmitters (50 or 5) in the case of the grid 105. Yes.
  • the moire pattern of the grid (5 lp / mm) is folded back at 3.33 lp / mm. 1.67 lp / mm.
  • the conventional X-ray diagnostic apparatus 101 having such a configuration has the following problems. That is, even if the moire pattern is removed from the image in which the moire pattern of the grid 105 is reflected by the moire pattern removal processing unit 121, there is a problem that the moire pattern remains in the image after the removal.
  • the present invention has been made in view of such circumstances, and an object thereof is to provide an X-ray diagnostic apparatus and an X-ray diagnostic program for removing a moire pattern remaining in a moire pattern removal image.
  • FIG. 13 shows an example of a one-dimensional Fourier transform result of an image. Specifically, a grid 105 of 50 lines / cm is used for an X-ray detector (FPD) 104 having a pixel pitch of 150 ⁇ m.
  • FPD X-ray detector
  • the peak frequency of the first harmonic of the moire pattern of the grid 105 is located at about 1.84 lp / mm, and the peak frequency of the second harmonic is located at about 2.98 lp / mm.
  • This frequency component is the third harmonic.
  • the third harmonic is three times (integer multiple) of the first harmonic and has a frequency of 5.52 lp / mm.
  • the frequency of 5.52 lp / mm is folded back at 3.33 lp / mm and appears at a position of 1.14 lp / mm.
  • the second harmonic wave has a frequency of 3.68 lp / mm which is twice that of the first harmonic wave, and is folded back and appears at a position of 2.98 lp / mm.
  • the third harmonic of the moire pattern was generated in the image after the moire pattern was removed due to the density error of the grid 105.
  • the grid density is 50 lines / cm, there are those with 45 to 55 lines / cm. This is because the grid 105 has a grid density tolerance of ⁇ 10% in the Japanese Industrial Standard (JIS).
  • FIG. 14 is a diagram showing the relationship between the grid error and the peak frequency of the harmonics of the moire pattern.
  • the vertical axis indicates the peak frequency (lp / mm) of the first to third harmonics
  • the X-ray detector (FPD) 104 is assumed to have a pixel pitch of 150 ⁇ m. In the case of the moire pattern shown in FIG.
  • the actual grid density is 48.2 line / cm
  • the first harmonic and the third harmonic are almost overlapped, and if the first harmonic is removed, Third harmonics are also removed.
  • a shift occurs in the peak frequencies of the first harmonic and the third harmonic. Therefore, the third frequency of the moire pattern remains in the image after removal.
  • the third harmonic on the low frequency side is removed particularly when the third harmonic is located on the low frequency side with respect to the other first and second harmonics. Then, information representing a relatively large shape in the image is removed at the same time as the third harmonic, and the visibility is greatly deteriorated. Further, if the third harmonic on the high frequency side is removed, information representing a fine shape or an edge in the image is removed at the same time as the third harmonic.
  • the third harmonic is a frequency component having a small spectral intensity, and therefore, the third harmonic is conspicuous in a flat region where the change in pixel value (luminance change) is small in the image, and is recognized outside the flat region. I found out that it was not.
  • an X-ray diagnostic apparatus includes an X-ray irradiation unit that irradiates a subject with X-rays, an X-ray detector that detects X-rays transmitted through the subject, and the X-rays of the X-ray detector.
  • an X-ray diagnostic apparatus that is arranged on the incident side and that removes scattered X-rays and acquires an image based on the output of the X-ray detector, the moiré pattern of the grid reflected in the image
  • a peak frequency detector for detecting a peak frequency of the first harmonic, and removing the first harmonic and the second harmonic of the moire pattern from the image based on the peak frequency of the first harmonic
  • a first and second harmonic removal processing unit for obtaining a second harmonic removal image
  • a peak frequency of the third harmonic of the moire pattern based on the peak frequency of the first harmonic, and the third harmonic 3rd harmonic extraction field for extracting
  • a third harmonic extraction filter creating unit for creating a third harmonic
  • a third harmonic extracting unit for extracting the third harmonic from the first and second harmonic removed images based on the third harmonic extraction filter
  • a third harmonic subtracting unit that subtracts the extracted third harmonic from the first and second harmonic removed images to obtain a third harmonic removed image.
  • the third harmonic extraction filter creation unit calculates the peak frequency of the third harmonic of the moire pattern based on the peak frequency of the first harmonic, and the third harmonic A third harmonic extraction filter for extracting.
  • the third harmonic extraction unit extracts the third harmonic from the first and second harmonic removed images based on the third harmonic extraction filter.
  • the third harmonic subtraction unit subtracts the extracted third harmonic from the first and second harmonic removal images to obtain a third harmonic removal image. Thereby, the third harmonic of the moire pattern remaining in the moire pattern removal image can be removed.
  • the X-ray diagnostic apparatus includes a flatness calculation unit that calculates flatness in the first and second harmonic removal images and obtains flatness information, and a flat region of the flatness information. It is preferable to include an output selection unit that outputs a third harmonic removal image and outputs the first and second harmonic removal images in a region other than the flat region of the flatness information. Since the third harmonic of the moiré pattern is a frequency component having a small spectral intensity, it is conspicuous in a flat area where the luminance change is small in the image, and is not recognized outside the flat area.
  • the output selection unit outputs the third harmonic removal image in the flat region out of the flatness information obtained by calculating the flatness by the flatness calculation unit, and the first in the region other than the flatness of the flatness information. -Output a second harmonic removal image. That is, the third harmonic removal image is output only in the flat region where the third harmonic of the moire pattern is conspicuous in the first and second harmonic removal images. Accordingly, it is possible to suppress the deterioration of the visibility of the image and to remove the third harmonic moire pattern remaining in the moire pattern removed image.
  • the flatness calculation unit calculates a luminance change amount between the target pixel to be processed and its surrounding pixels, and acquires flatness information. It is preferable to provide. Thereby, the luminance change amount calculation unit acquires flatness information indicating the flatness in the first and second harmonic removed images by calculating the luminance change amount between the target pixel and the surrounding pixels. Can do.
  • the luminance change amount calculation unit obtains flatness information by subtracting an average value of the target pixel and surrounding pixels from a pixel value of the target pixel. It is preferable. Thereby, the luminance change amount calculation unit can acquire the flatness information by calculating the luminance change amount by subtracting the addition average value of the target pixel and surrounding pixels from the pixel value of the target pixel. .
  • the flatness calculation unit divides the luminance change amount of each pixel of the flatness information by the pixel value of the corresponding first and second harmonic removal images. It is preferable to provide a normalization processing unit that performs normalization processing. That is, a pixel having a high pixel value has a large luminance change amount, and a pixel having a low pixel value tends to have a small luminance change amount. Even if the brightness change amount is the same, the weight change is different between the brightness change amount of the high brightness pixel and the brightness change amount of the low brightness pixel. Therefore, by dividing the luminance change amount of each pixel by the pixel value of the corresponding pixel, the luminance change amounts of the high luminance pixel and the low luminance pixel can be compared on the same basis.
  • the flatness calculation unit performs an average of luminance change amounts for each region obtained by dividing the flatness information by a plurality of pixels, and divides the average value. It is preferable to include a smoothing processing unit that applies the smoothing to each pixel in the region. By smoothing the luminance change amount, output selection by the output selection unit can be performed in an arbitrary region. In addition, a region determined as a flat region scattered in a region other than the flat region can be processed as a region other than the flat region. If the third harmonic wave-removed image is output to a flat region scattered in a region other than the flat region, there is a risk that visibility may be deteriorated.
  • the X-ray diagnostic apparatus further includes a low-pass filter processing unit that performs one-dimensional low-pass filter processing in a direction along the third harmonic fringes extracted by the third harmonic extraction unit. It is preferable. It is possible to remove a pattern in an oblique direction with respect to the third harmonic fringe extracted in addition to the third harmonic of the moire pattern.
  • An X-ray diagnostic program is an X-ray diagnostic program for causing a computer to execute a process of processing an image taken through a grid that removes scattered X-rays.
  • a step of detecting a peak frequency of the first harmonic of the reflected moire pattern of the grid, and removing a first harmonic and a second harmonic of the moire pattern from the image based on the peak frequency of the first harmonic And acquiring the first and second harmonic removal images, calculating the peak frequency of the third harmonic of the moire pattern based on the peak frequency of the first harmonic, and extracting the third harmonic
  • a third harmonic extraction filter for extracting the third harmonic from the first and second harmonic removal images based on the third harmonic extraction filter, extracted Subtracting the third harmonic from the first and second harmonic removed images to obtain a third harmonic removed image, and causing a computer to execute these steps. is there.
  • the peak frequency of the third harmonic of the moire pattern is calculated based on the peak frequency of the first harmonic of the moire pattern, and the third harmonic is extracted.
  • Create a third harmonic extraction filter Next, the third harmonic is extracted from the first and second harmonic removed images based on the third harmonic extraction filter. Then, the extracted third harmonic is subtracted from the first and second harmonic removed images to obtain a third harmonic removed image. Thereby, the third harmonic of the moire pattern remaining in the moire pattern removal image can be removed.
  • the peak frequency of the third harmonic of the moire pattern is calculated based on the peak frequency of the first harmonic of the moire pattern, and the third harmonic is calculated.
  • a third harmonic extraction filter for extraction is created.
  • the third harmonic is extracted from the first and second harmonic removed images based on the third harmonic extraction filter.
  • the extracted third harmonic is subtracted from the first and second harmonic removed images to obtain a third harmonic removed image.
  • the third harmonic of the moire pattern remaining in the moire pattern removal image can be removed.
  • FIG. 1 It is a block diagram which shows schematic structure of the X-ray diagnostic apparatus which concerns on an Example. It is a figure where it uses for description of a moire pattern. It is a figure where it uses for description of the extraction filter preparation of the 3rd harmonic of a moire pattern. It is a figure which shows an example of the 3rd harmonic of a moire pattern.
  • A) is a figure where it uses for description of a vertical low-pass filter (LPF) process
  • (b) is a figure which shows an example of a vertical low-pass filter (LPF).
  • LPF vertical low-pass filter
  • FIG. 1 is a block diagram illustrating a schematic configuration of an X-ray diagnostic apparatus according to an embodiment
  • FIG. 2 is a diagram for explaining a moire pattern
  • FIG. 3 is a diagram for explaining the creation of an extraction filter for the third harmonic of the moire pattern
  • FIG. 4 is a diagram illustrating an example of the third harmonic of the moire pattern.
  • FIG. 5A is a diagram for explaining the vertical low-pass filter (LPF) processing
  • FIG. 5B is a diagram illustrating an example of the vertical low-pass filter (LPF).
  • FIG. 6 is a diagram for explaining the third harmonic removal of the moire pattern
  • FIG. 7 is a diagram for explaining the calculation of the luminance change amount.
  • FIG. 8 is a diagram for explaining the smoothing process
  • FIG. 9 is a diagram for explaining the output selection.
  • the X-ray diagnostic apparatus 1 is arranged to face the X-ray tube 3 with the subject M sandwiched between the top plate 2 on which the subject M is placed, the X-ray tube 3 that irradiates the subject M with X-rays, and the subject M. And an X-ray detector 4 for detecting X-rays transmitted through the subject M, and a grid 5 disposed on the X-ray incident side of the X-ray detector 4 for removing scattered X-rays.
  • the X-ray tube 3 corresponds to the X-ray irradiation unit in the present invention.
  • the X-ray tube 3 is subjected to control necessary for X-ray irradiation by the X-ray tube control unit 6.
  • the X-ray tube controller 6 has a high voltage generator 7 that generates the tube voltage and tube current of the X-ray tube 3.
  • the X-ray tube control unit 6 irradiates X-rays from the X-ray tube 3 in accordance with X-ray conditions such as a tube voltage, a tube current, and an irradiation time set by an input unit 12 described later.
  • the X-ray detector 4 detects X-rays emitted from the X-ray tube 3 and outputs an X-ray detection signal that is a signal corresponding to the X-ray intensity distribution.
  • the X-ray detector 4 is configured by a flat panel X-ray detector (FPD), an image intensifier, or the like.
  • the X-ray detector 4 is configured with, for example, a 0.15 um pixel pitch and has a resolution of 3.33 lp / mm.
  • the grid 5 is configured with a grid density of 50 lines / cm (5 lp / mm), for example.
  • the A / D converter 8 and the image processing unit 9 are provided in the subsequent stage of the X-ray detector 4.
  • the A / D converter 8 converts the analog X-ray detection signal output from the X-ray detector 4 into a digital signal.
  • the image processing unit 9 performs various processes on the image based on the digital X-ray detection signal output from the X-ray detector 4.
  • the X-ray diagnostic apparatus 1 includes a main control unit 10 that comprehensively controls each component of the apparatus 1, a display unit 11 that displays an image processed by the image processing unit 9, and an operator that performs input settings and An input unit 12 for performing various operations and a storage unit 13 for storing processed images are provided.
  • the main control unit 10 is composed of a central processing unit (CPU) and the like, and executes various programs. For example, the main control unit 10 performs control to move the top 2, the X-ray tube 3 or the X-ray detector 4 on which the subject M is placed to a predetermined position.
  • the display unit 11 includes a monitor or the like.
  • the input unit 12 includes a keyboard, a mouse, and the like.
  • the storage unit 13 includes a storage medium such as a ROM (Read-only Memory), a RAM (Random-Access Memory), or a hard disk.
  • the image processing unit 9 includes a first and second harmonic removal processing unit 21 and a third harmonic removal for removing the moire pattern of the grid 5 reflected in the image acquired based on the output of the X-ray detector 4. And a processing unit 23.
  • the first and second harmonic removal processing unit 21 includes a peak frequency detection unit 25 that detects the peak frequency of the first harmonic of the moire pattern of the grid 5 reflected in the image.
  • the first and second harmonic removal processing unit 21 removes the first and second harmonics of the moire pattern from the image based on the peak frequency of the first harmonic to remove the first and second harmonic removal images. To get.
  • the first / second harmonic removal processing unit 21 removes the first harmonic and the second harmonic of the moire pattern of the grid 5 reflected in the image.
  • the first and second harmonic removal processing unit 21 first applies a horizontal direction (a direction orthogonal to the fringes of the moire pattern) to the image in which the moire pattern (code mp) is reflected.
  • a one-dimensional Fourier transform is performed on the pixel row (reference L1).
  • the one-dimensional Fourier transform is performed in order from, for example, the pixel row at the upper end of the image.
  • the peak frequency detector 25 detects the peak frequency indicating the first harmonic of the moire pattern from the Fourier transform result (FIG. 13).
  • the peak frequency detection unit 25 detects a predetermined peak frequency set from a plurality of peak frequencies detected in order from the upper end of the image. Since the second harmonic is twice (integer multiple) the first harmonic, it is calculated from the peak frequency of the first harmonic. Based on the acquired peak frequencies of the first harmonic and the second harmonic, a filter for extracting the first harmonic and the second harmonic of the moire pattern is created. The first harmonic and the second harmonic are extracted using this filter. Then, the first and second harmonics are removed from the image by subtracting the extracted first and second harmonics from the image in which the moire pattern is reflected. A harmonic removal image is acquired. In FIG. 2, the symbol f indicates a flat region of the subject M in the image.
  • the third harmonic removal processing unit 23 includes a storage unit 27, a third harmonic extraction filter creation unit 29, a third harmonic extraction unit 31, a vertical low-pass filter processing unit (vertical LPF processing unit). 33, a third harmonic subtracting unit 35, a flatness calculating unit 37, and an output selecting unit 39.
  • the first and second harmonic removal images from which the first and second harmonics of the moire pattern have been removed by the first and second harmonic removal processing unit 21 are stored in the storage unit 27.
  • the third harmonic extraction filter creation unit 29 calculates the peak frequency of the third harmonic of the moire pattern based on the peak frequency of the first harmonic, and extracts the third harmonic. Create The peak frequency of the first harmonic detected by the peak frequency detector 25 is given as the peak frequency of the first harmonic.
  • the third harmonic is folded at a position of 3.33 lp / mm.
  • the case where the peak frequency of the first harmonic is 1.72 lp / mm is when the grid error (ratio) is 0.99 in FIG.
  • the third harmonic extraction filter is created by giving a predetermined band to the calculated peak frequency of the third harmonic.
  • the third harmonic extraction unit 31 extracts the third harmonic from the first and second harmonic removal images based on the third harmonic extraction filter.
  • the first and second harmonic removed images are given by the storage unit 27.
  • a filter process using a third harmonic extraction filter is performed on the pixel rows in the horizontal direction of the first and second harmonic removed images. This extraction processing is performed, for example, one column at a time in order from the upper end of the first and second harmonic removed images.
  • the third harmonic image of the moire pattern as shown in FIG. 4 is extracted.
  • a pattern in an oblique direction is extracted in addition to the third harmonic of the moire pattern, as indicated by reference numeral n in FIG.
  • n the third harmonic of the moire pattern
  • the vertical LPF processing unit 33 performs one-dimensional low-pass filter (LPF) processing in the vertical direction along the third harmonic fringes extracted by the third harmonic extraction unit 31.
  • the fringes of the third harmonic of the moiré pattern are indicated by reference numeral mp3 in FIG. That is, LPF processing is performed in the vertical direction along the third harmonic fringes.
  • the filter process shown in FIG. 5B is performed on each pixel of the extracted third harmonic image.
  • the filter illustrated in FIG. 5B is a vertical pixel row, and is performed using, for example, 17 pixels including 8 pixels above and below the target pixel to be processed.
  • the vertical LPF processing unit 33 corresponds to the low pass filter processing unit in the present invention.
  • the third harmonic fringes of the moire pattern are arranged in the horizontal direction, but when the fringes are arranged in the vertical direction, a one-dimensional low-pass filter is formed in the horizontal direction along the fringes. (LPF) processing is performed.
  • the third harmonic subtracting unit 35 subtracts the extracted third harmonic from the first and second harmonic removed images to obtain a third harmonic removed image. That is, as shown in FIG. 6, the third harmonic subtraction unit 35 subtracts the third harmonic image P2 extracted from the first and second harmonic removal images P1. Thereby, the 3rd harmonic removal image P3 is acquired.
  • the flatness calculation unit 37 calculates flatness in the first and second harmonic removal images and acquires flatness information.
  • the calculation of the flatness information by the flatness calculation unit 37 is performed in parallel with the process of acquiring the third harmonic removal image for extracting and removing the third harmonic.
  • the flatness information is information for identifying a flat area (flat area) and a non-flat structure area (an area other than the flat area) in the first and second harmonic removal images. For example, in imaging of the leg of the subject M, a flat portion includes a meat portion, and a region other than the flat region includes a bone portion.
  • the flatness calculation unit 37 includes a luminance change amount calculation unit 41, a normalization processing unit 43, and a smoothing processing unit 45.
  • the luminance change amount calculation unit 41 calculates the luminance change amount between the target pixel to be processed and its surrounding pixels, and acquires flatness information. Specifically, the luminance change amount calculation unit 41 obtains flatness information by subtracting the addition average value of the target pixel and the surrounding pixels from the pixel value of the target pixel. That is, as shown in FIG. 7, the luminance change amount calculation unit 41 performs an average of 21 ⁇ 21 pixels of the target pixel, the target pixel, and, for example, 10 pixels in the upper, lower, left, and right directions, and calculates the average from the pixel value of the target pixel. Subtract the value.
  • a pixel having a high pixel value has a large amount of change in luminance (a variation amount of the pixel value), and a pixel having a low pixel value tends to have a small amount of change in luminance. For example, it is assumed that there is a target pixel having a pixel value of 1000 and a target pixel having a pixel value of 10. When both of the luminance change amounts are 2, the luminance change amounts are the same value.
  • the luminance change amount between the high luminance pixel and the low luminance pixel can be compared on the same basis.
  • the smoothing processing unit 45 performs addition averaging of the luminance change amount for each region obtained by dividing the flatness information by a plurality of pixels, and smoothes the addition average value by applying it to each pixel of the divided region.
  • the smoothing processing unit 45 performs addition averaging, for example, by 16 ⁇ 16 pixels and reduces it to one pixel. Then, reduced flatness information P5 of 4 ⁇ 4 pixels is acquired. Then, the flatness information P6 that has been smoothed is obtained by returning (enlarging) the flatness information P5 of 4 ⁇ 4 pixels to the original size of 64 ⁇ 64 pixels.
  • the output selection unit 39 outputs the third harmonic removal image in the flat region of the flatness information, and outputs the first and second harmonic removal images in the region other than the flat region of the flatness information. That is, as shown in FIG. 9, the output selection unit 39 outputs either one of the first / second harmonic removal image P1 and the third harmonic removal image P3 based on the flatness information P6. Select whether or not. If the luminance change amount subjected to the normalization process and the smoothing process in the flatness information P6 is equal to or greater than the threshold value, it is identified as a flat area, and if it is less than the threshold value, it is identified as an area other than the flat area.
  • the output selection unit is in units of the 16 ⁇ 16 pixel area. 39 is selected and output. As shown in FIG. 9, in the flat region (symbol f) of the flatness information P6, the third harmonic removal image P3 is output. On the other hand, in regions other than the flat region of the flatness information P6, the first and second harmonic removed images P1 are output.
  • Step S01 X-ray Imaging
  • the X-ray diagnostic apparatus 1 uses an X-ray tube according to X-ray conditions such as tube voltage, tube current and irradiation time set by the X-ray tube controller 6. 3 is irradiated with X-rays. X-rays irradiated from the X-ray tube 3 pass through the subject M and enter the X-ray detector 4. The X-ray detector 4 outputs an X-ray detection signal corresponding to the X-ray intensity distribution.
  • the A / D converter 8 converts the analog X-ray detection signal output from the X-ray detector 4 into a digital signal. An image obtained by photographing the subject M can be acquired based on the digitally converted X-ray detection signal.
  • Step S02 First / Second Harmonic Elimination Processing
  • the peak frequency detection unit 25 detects the peak frequency of the first harmonic of the moiré pattern of the grid reflected in the image.
  • the first and second harmonic removal processing unit 21 removes the first and second harmonics of the moire pattern from the image based on the peak frequency of the first harmonic to remove the first and second harmonic removal images. To get.
  • Step S03 Creation of Third Harmonic Extraction Filter
  • the third harmonic extraction filter creation unit 29 calculates the peak frequency of the third harmonic of the moire pattern based on the peak frequency of the first harmonic, and the third harmonic.
  • a third harmonic extraction filter for extracting waves is created (FIG. 3).
  • the third harmonic extraction unit 31 extracts the third harmonic from the first and second harmonic removed images based on the third harmonic extraction filter (FIG. 5A )).
  • the extraction of the third harmonic is performed, for example, by performing a filter process using the third harmonic extraction filter for each column in order from the upper end of the horizontal pixel column of the first and second harmonic removed images.
  • the vertical LPF processing unit 33 is a one-dimensional low-pass filter in the vertical direction along the stripes of the third harmonic (FIG. 5A) extracted by the third harmonic extraction unit 31.
  • (LPF) processing is performed.
  • the filter shown in FIG. 5B that is, 17 pixels including the upper and lower 8 pixels of the target pixel to be processed is used.
  • a new pixel value of the target pixel is set.
  • Step S06 Third harmonic removal processing (subtraction) As shown in FIG. 6, the third harmonic subtraction unit 35 subtracts the extracted third harmonic P2 from the first / second harmonic removal image P1 to obtain a third harmonic removal image P3.
  • Steps S13 to S15 for calculating the flatness in the first and second harmonic removal images and acquiring the flatness information are performed.
  • Step S13 Calculation of luminance change amount
  • the normalization processing unit 43 performs normalization processing by dividing the luminance change amount of each pixel of the flatness information by the pixel value of the corresponding first and second harmonic removed images. That is, the luminance change amount C1 is divided by the pixel value B of the target pixel used in the calculation of the luminance change amount C1.
  • the smoothing processing unit 45 performs addition averaging of the luminance change amount for each region obtained by dividing the flatness information by a plurality of pixels, and applies the addition average value to each pixel of the divided region. To smooth. That is, as shown in FIG. 8, the smoothing processing unit 45 performs addition averaging with, for example, 16 ⁇ 16 pixels, and reduces it to one pixel. Then, when returning to the original size, the luminance change amount obtained by adding and averaging all the pixels is given to the original size of 16 ⁇ 16 pixels.
  • Step S07 Output Selection
  • the output selection unit 39 outputs the third harmonic removal image P3 in the flat region f of the flatness information P6, and a region other than the flat region f of the flatness information P6. Then, the first and second harmonic removed images P1 are output. Thus, in the flat region f, the output selection image P7 from which the third harmonic of the moire pattern has been removed is acquired.
  • Step S08 Display / Save The output selection image P7 selected for output in step S07 is subjected to other necessary processing in the image processing section 9.
  • the output selection image P ⁇ b> 7 is displayed on the display unit 11 and stored in the storage unit 13.
  • the third harmonic extraction filter creation unit 29 calculates the peak frequency of the third harmonic of the moire pattern based on the peak frequency of the first harmonic, A third harmonic extraction filter for extracting the third harmonic is created.
  • the third harmonic extraction unit 31 extracts the third harmonic from the first and second harmonic removal images based on the third harmonic extraction filter.
  • the 3rd harmonic subtraction part 35 subtracts the extracted 3rd harmonic from a 1st, 2nd harmonic removal image, and acquires a 3rd harmonic removal image. Thereby, the third harmonic of the moire pattern remaining in the first and second harmonic removal images can be removed.
  • the X-ray diagnostic apparatus 1 calculates the flatness in the first and second harmonic removal images and acquires the flatness information, and the third harmonic in the flat region of the flatness information.
  • An output selection unit 39 is provided that outputs the removed image and outputs the first and second harmonic removed images in a region other than the flat region of the flatness information. Since the third harmonic of the moiré pattern is a frequency component having a small spectral intensity, it is conspicuous in a flat area where the luminance change is small in the image, and is not recognized outside the flat area. Therefore, the output selection unit 39 outputs the third harmonic removal image in the flat region among the flatness information obtained by calculating the flatness by the flatness calculation unit 37, and in the region other than the flatness of the flatness information.
  • an image having no remaining moire pattern can be obtained without deteriorating the visibility except in a flat region.
  • the luminance change is small, and even if the third harmonic of the moire pattern is removed, it is relatively difficult to be affected by the visibility deterioration. Further, since the third harmonic of the moire pattern in the flat region is removed, it is possible to improve the S / N ratio in the flat region.
  • the luminance change amount calculation unit 41 acquires flatness information indicating the flatness in the first and second harmonic removed images by calculating the luminance change amount between the target pixel and the surrounding pixels. Can do. Further, the luminance change amount calculation unit 41 can acquire the flatness information by calculating the luminance change amount by subtracting the addition average value of the target pixel and the surrounding pixels from the pixel value of the target pixel. .
  • the normalization processing unit 43 performs normalization processing by dividing the luminance change amount of each pixel of the flatness information by the pixel value of the corresponding first / second harmonic removed image. That is, a pixel having a high pixel value has a large luminance change amount, and a pixel having a low pixel value tends to have a small luminance change amount. Even if the brightness change amount is the same, the weight change is different between the brightness change amount of the high brightness pixel and the brightness change amount of the low brightness pixel. Therefore, by dividing the luminance change amount of each pixel by the pixel value of the corresponding pixel, the luminance change amounts of the high luminance pixel and the low luminance pixel can be compared on the same basis.
  • the smoothing processing unit 45 performs addition averaging of the luminance change amount for each region obtained by dividing the flatness information by a plurality of pixels (for example, 16 ⁇ 16 pixels), and each pixel of the region obtained by dividing the addition average value. (16 ⁇ 16 pixels) is applied for smoothing.
  • Image output selection by the output selection unit 39 is not performed on a pixel basis, but can be performed in an arbitrary region by smoothing.
  • a region determined as a flat region scattered in a region other than the flat region can be processed as a region other than the flat region. If the third harmonic wave-removed image is output to a flat region scattered in a region other than the flat region, there is a risk that visibility may be deteriorated.
  • the vertical LPF processing unit 33 performs one-dimensional LPF processing in the vertical direction along the third harmonic fringes extracted by the third harmonic extraction unit 31. It is possible to remove a pattern in an oblique direction with respect to the third harmonic fringe extracted in addition to the third harmonic of the moire pattern.
  • the present invention is not limited to the above embodiment, and can be modified as follows.
  • the image processing unit 9 includes a control unit configured by a CPU or the like for executing a program, and a storage unit configured by a storage medium such as a ROM or RAM that stores the program or the like. May be provided.
  • the storage unit may store a program for the operations of steps S01 to S08 and S13 to S15, and the program may be executed by the control unit. In this case, operations necessary for this program are input through the input unit 12, and the output selection image P ⁇ b> 7 is displayed on the display unit 11.
  • each of steps S01 to S08 and S13 to S15 is not limited to the image processing unit 9, and the operation program is stored in the storage unit 13 and executed by the main control unit 10. You may do it. In this case, operations necessary for this program are input through the input unit 12, and the output selection image P ⁇ b> 7 is displayed on the display unit 11.
  • the operation program may be executed on a personal computer connected to the X-ray diagnostic apparatus 1 in a network system such as a LAN.

Landscapes

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

Abstract

　Ｘ線診断装置は、画像に映りこんだ前記グリッドのモアレパターンの第１高調波のピーク周波数を検出するピーク周波数検出部と、第１高調波のピーク周波数に基づいて画像からモアレパターンの第１高調波と第２高調波を除去して第１・第２高調波除去画像を取得する第１・第２高調波除去処理部と、第１高調波のピーク周波数に基づいてモアレパターンの第３高調波のピーク周波数を算出し、第３高調波を抽出するための第３高調波抽出フィルタを作成する第３高調波抽出フィルタ作成部と、第３高調波抽出フィルタに基づいて第１・第２高調波除去画像から第３高調波を抽出する第３高調波抽出部と、抽出された第３高調波を第１・第２高調波除去画像から減算して第３高調波除去画像を取得する第３高調波減算部とを備えている。

Description

Ｘ線診断装置およびＸ線診断用プログラム
　本発明は、Ｘ線画像に映りこんだグリッドのモアレパターンを除去するＸ線診断装置およびＸ線診断用プログラムに関する。
　従来のＸ線診断装置は、図１１に示すように、Ｘ線を被検体Ｍに照射するＸ線管１０３と、被検体Ｍを透過したＸ線を検出するＸ線検出器（例えば、フラットパネル型Ｘ線検出器（ＦＰＤ））１０４と、Ｘ線検出器１０４のＸ線入射側に設けられ、被検体Ｍを透過するとき等に散乱するＸ線を除去する散乱線除去グリッド（以下、「グリッド」と略す）１０５とを備えている。グリッド１０５は、Ｘ線を吸収する吸収体（例えば鉛）と、Ｘ線を透過する透過体（例えばアルミニウムや空気）とが交互に配設して構成される。このグリッド１０５によって、散乱して斜め方向から入射するＸ線は吸収体で吸収され、透過体を透過したＸ線のみがＸ線検出器１０４で検出される。それにより、鮮明な画像を取得することができる。
　しかしながら、従来のＸ線診断装置１０１では、Ｘ線検出器１０４の解像度とグリッド１０５の密度との違いにより、取得した画像にグリッド１０５のモアレパターンが映りこんでしまう問題がある。このグリッド１０５のモアレパターンの対策については、種々の提案がなされている（例えば、特許文献１参照）。
　また、図１１に示す従来のＸ線診断装置１０１では、モアレパターン除去処理部１２１を備えて、画像に映りこんだグリッド１０５のモアレパターンの第１高調波と第２高調波を除去している。まず、モアレパターンが映りこんだ画像に対して１次元のフーリエ変換を行い、そのフーリエ変換結果からモアレパターンの第１高調波を示すピーク周波数を検出する。第２高調波は、第１高調波の２倍（整数倍）であるので、第１高調波のピーク周波数により算出される。取得された第１高調波と第２高調波のピーク周波数に基づき、モアレパターンの第１高調波と第２高調波を抽出するためのフィルタを作成し、このフィルタを用いて第１高調波と第２高調波を抽出する。そして、モアレパターンが映りこんだ画像から第１高調波と第２高調波を減算することでモアレパターンの第１高調波と第２高調波を除去する。
　なお、Ｘ線検出器１０４の解像度は、ナイキスト周波数Ｎｙで表すことができる。ナイキスト周波数Ｎｙは、ピクセルピッチδ（mm）でサンプリングできる空間周波数の限界を示す。ナイキスト周波数Ｎｙは、次の式（１）から求められる。
　Ｎｙ＝１／（２×δ）…（１）
　例えば、Ｘ線検出器１０４は、ピクセルピッチδが１５０umで構成されているとすると、ナイキスト周波数Ｎｙは、Ｎｙ＝１／（２×０．１５）＝３．３３…lp/mmとなる。すなわちＸ線検出器１０４の解像度は、３．３３…lp/mmである。なお、Ｘ線検出器１０４の解像度は、以下適宜「３．３３lp/mm」と略する。
　また、グリッド１０５の密度は、例えば５０line/cmで構成される。この５０line/cmは、５line/mmや５lp/mmなどとも表される。「line/cm」や「lp/mm」等は、グリッド１０５の場合、単位長さ（１cmや1mm）当たり吸収体の本数や吸収体と透過体の組の本数（５０や５）を示している。解像度が３．３３lp/mmのＸ線検出器と、密度が５lp/mm（５０line/cm）のグリッドとを用いる場合、グリッド（５lp/mm）のモアレパターンは、３．３３lp/mmで折り返して、１．６７lp/mmの位置に生じる。すなわち、３．３３…lp/mm－（５lp/mm－３．３３…lp/mm）＝１．６７（１．６６６…）lp/mmの位置にグリッドのモアレパターンの第１高調波が現れる（図１２）。
特開２００２－１５２４６７号公報
　しかしながら、このような構成を有する従来のＸ線診断装置１０１では、次のような問題がある。すなわち、モアレパターン除去処理部１２１により、グリッド１０５のモアレパターンが映りこんだ画像からモアレパターンを除去しても、除去後の画像にモアレパターン残りが存在してしまうという問題がある。
　本発明は、このような事情に鑑みてなされたものであって、モアレパターン除去画像に残るモアレパターンを除去するＸ線診断装置およびＸ線診断用プログラムを提供することを目的とする。
　発明者は、上記の問題を解決するために鋭意研究した結果、次のような知見を得た。まず、モアレパターン除去画像に残るモアレパターンは、第３高調波であることがわかった。図１３は、画像の１次元のフーリエ変換結果の一例を示す。なお、具体的には、１５０umのピクセルピッチのＸ線検出器（ＦＰＤ）１０４に対して、５０line/cmのグリッド１０５を用いている。図１３において、グリッド１０５のモアレパターンの第１高調波のピーク周波数は約１．８４lp/mmに位置し、第２高調波のピーク周波数は約２．９８lp/mmに位置する。そして、約１．１４lp/mmに位置するピーク周波数を有する周波数成分が存在する。この周波数成分が第３高調波である。すなわち、第３高調波は、第１高調波の３倍（整数倍）であり、５．５２lp/mmの周波数を有する。５．５２lp/mmの周波数は、３．３３lp/mmで折り返されて１．１４lp/mmの位置に現れる。なお、第２高調波は、第１高調波の２倍の３．６８lp/mmの周波数を有し、折り返されて２．９８lp/mmの位置に現れる。
　また、グリッド１０５の密度の誤差により、モアレパターン除去後の画像にモアレパターンの第３高調波が発生することがわかった。例えば、グリッド密度が５０line/cmの場合であっても４５～５５line/cmのものが存在する。なぜなら、グリッド１０５には、日本工業規格(JIS: Japanese Industrial Standard)において、グリッド密度の許容差が±１０％の範囲で認められているからである。
　図１４は、グリッド誤差とモアレパターンの高調波のピーク周波数との関係を示す図である。具体的には、図１４において、縦軸には第１～３高調波のピーク周波数（lp/mm）を示し、横軸には、５０line/cmのグリッドと実際に誤差があるグリッドとの比を示す。例えば、実際のグリッド密度が５１line/cmの場合、グリッド誤差（比）は、５１／５０＝１．０２である。なお、Ｘ線検出器（ＦＰＤ）１０４はピクセルピッチが１５０umで構成されているものとする。なお、上述の図１３に示すモアレパターンの場合は、実際のグリッド密度は４８．２line/cmであり、グリッド誤差（比）は、４８．２／５０＝０．９６である（符号Ｅ）。図１４において、５０line/cmのグリッド１０５のライン数に誤差が無い場合（符号Ｒ）、第１高調波と第３高調波とは、ほぼ重なっており、第１高調波を除去すれば、同時に第３高調波も除去される。しかしながら、グリッド密度に誤差がある場合は、第１高調波と第３高調波のピーク周波数にずれが生じる。そのため、モアレパターンの第３周波数が除去後の画像に残ってしまう。
　また、モアレパターンの第１高調波と第２高調波を除去して、更に第３高調波を除去すると、画像内のより多くの情報が同時に除去されることになる。そのため、単に第３高調波を除去しただけでは、画像（画像内の被写体）を明確に認識できるか否かの画像の視認性を劣化させてしまう。図１３や図１４に示すように、第３高調波が他の第１高調波や第２高調波に対して低周波数側に位置する場合など、特に、低周波側の第３高調波を除去すると、画像内の比較的大きな形状を表す情報が第３高調波と同時に除去されてしまい、視認性の劣化が大きい。また、高周波数側の第３高調波を除去すると、画像内の細かな形状やエッジを表す情報が第３高調波と同時に除去されてしまう。また、第３高調波は、図１３に示すように、スペクトル強度が小さい周波数成分であるので、画像内では、画素値の変化（輝度変化）の少ない平坦領域に目立ち、その平坦領域以外では認識されないことがわかった。
　このような知見に基づく本発明は、次のような構成をとる。すなわち、本発明に係るＸ線診断装置は、被検体にＸ線を照射するＸ線照射部と、被検体を透過したＸ線を検出するＸ線検出器と、前記Ｘ線検出器のＸ線入射側に配置され、散乱したＸ線を除去するグリッドとを備え、前記Ｘ線検出器の出力に基づいて画像を取得するＸ線診断装置において、前記画像に映りこんだ前記グリッドのモアレパターンの第１高調波のピーク周波数を検出するピーク周波数検出部と、前記第１高調波のピーク周波数に基づいて前記画像からモアレパターンの第１高調波と第２高調波を除去して第１・第２高調波除去画像を取得する第１・第２高調波除去処理部と、前記第１高調波のピーク周波数に基づいてモアレパターンの第３高調波のピーク周波数を算出し、前記第３高調波を抽出するための第３高調波抽出フィルタを作成する第３高調波抽出フィルタ作成部と、前記第３高調波抽出フィルタに基づいて前記第１・第２高調波除去画像から前記第３高調波を抽出する第３高調波抽出部と、抽出された前記第３高調波を前記第１・第２高調波除去画像から減算して第３高調波除去画像を取得する第３高調波減算部と、を備えていることを特徴とするものである。
　本発明に係るＸ線診断装置によれば、第３高調波抽出フィルタ作成部は、第１高調波のピーク周波数に基づいてモアレパターンの第３高調波のピーク周波数を算出し、第３高調波を抽出するための第３高調波抽出フィルタを作成する。第３高調波抽出部は、第３高調波抽出フィルタに基づいて第１・第２高調波除去画像から第３高調波を抽出する。第３高調波減算部は、抽出された第３高調波を第１・第２高調波除去画像から減算して第３高調波除去画像を取得する。それにより、モアレパターン除去画像に残っているモアレパターンの第３高調波を除去することができる。
　また、本発明に係るＸ線診断装置は、前記第１・第２高調波除去画像内の平坦度を算出して平坦度情報を取得する平坦度算出部と、前記平坦度情報の平坦領域では第３高調波除去画像を出力し、前記平坦度情報の平坦領域以外の領域では前記第１・第２高調波除去画像を出力する出力選択部とを備えていることが好ましい。モアレパターンの第３高調波は、スペクトル強度が小さい周波数成分であるので、画像内では、輝度変化の少ない平坦領域に目立ち、その平坦領域以外では認識されない。そのため、出力選択部は、平坦度算出部で平坦度を算出して取得した平坦度情報のうち、平坦領域では第３高調波除去画像を出力し、平坦度情報の平坦以外の領域では第１・第２高調波除去画像を出力する。すなわち、第１・第２高調波除去画像にてモアレパターンの第３高調波が目立つ平坦領域のみ、第３高調波除去画像を出力させている。それにより、画像の視認性の劣化を抑制するとともにモアレパターン除去画像に残っていた第３高調波のモアレパターンを除去することができる。
　また、本発明に係るＸ線診断装置において、前記平坦度算出部は、処理対象である注目画素とその周囲の画素との輝度変化量を算出して平坦度情報を取得する輝度変化量算出部を備えていることが好ましい。それにより、輝度変化量算出部は、第１・第２高調波除去画像内の平坦度を示す平坦度情報を、注目画素とその周囲の画素との輝度変化量を算出することで取得することができる。
　また、本発明に係るＸ線診断装置において、前記輝度変化量算出部は、前記注目画素の画素値から前記注目画素とその周囲の画素との加算平均値を減算して平坦度情報を取得することが好ましい。それにより、輝度変化量算出部は、平坦度情報を、注目画素の画素値から注目画素とその周囲の画素との加算平均値を減算して輝度変化量を算出することで取得することができる。
　また、本発明に係るＸ線診断装置において、前記平坦度算出部は、前記平坦度情報の各画素の輝度変化量をそれぞれ対応する前記第１・第２高調波除去画像の画素値で除算して正規化処理する正規化処理部を備えていることが好ましい。すなわち、画素値が高い画素では輝度変化量が大きく、画素値が低い画素では輝度変化量が小さい傾向にある。輝度変化量が同じであっても高輝度画素の輝度変化量と低輝度画素の輝度変化量とでは、その重みが異なる。そのため、各画素の輝度変化量をそれぞれ対応する画素の画素値で除算することにより、高輝度画素と低輝度画素との輝度変化量を同じ基準で比較することができる。
　また、本発明に係るＸ線診断装置において、前記平坦度算出部は、前記平坦度情報を複数個の画素で分割した領域ごとに輝度変化量の加算平均を行い、その加算平均値を分割した領域の各画素に適用させて平滑化する平滑化処理部を備えていることが好ましい。輝度変化量を平滑化することにより任意の領域で出力選択部による出力選択を行うことができる。また、平坦領域以外の領域内に点在する平坦領域と判定された領域を平坦領域以外の領域として処理することができる。平坦領域以外の領域内に点在する平坦領域まで第３高調波除去画像を出力させると視認性劣化を招くおそれがある。
　また、本発明に係るＸ線診断装置は、前記第３高調波抽出部で抽出された前記第３高調波の縞目に沿った方向に１次元のローパスフィルタ処理を行うローパスフィルタ処理部を備えていることが好ましい。モアレパターンの第３高調波の他に抽出される第３高調波の縞目に対して斜め方向のパターンを除去することができる。
　また、本発明に係るＸ線診断用プログラムは、散乱したＸ線を除去するグリッドを介して撮影された画像の処理工程をコンピュータに実行させるためのＸ線診断用プログラムであって、前記画像に映りこんだ前記グリッドのモアレパターンの第１高調波のピーク周波数を検出する工程と、前記第１高調波のピーク周波数に基づいて前記画像からモアレパターンの第１高調波と第２高調波を除去して第１・第２高調波除去画像を取得する工程と、前記第１高調波のピーク周波数に基づいてモアレパターンの第３高調波のピーク周波数を算出し、前記第３高調波を抽出するための第３高調波抽出フィルタを作成する工程と、前記第３高調波抽出フィルタに基づいて前記第１・第２高調波除去画像から前記第３高調波を抽出する工程と、抽出された前記第３高調波を前記第１・第２高調波除去画像から減算して第３高調波除去画像を取得する工程と、を備え、これらの工程をコンピュータに実行させることを特徴とするものである。
　本発明に係るＸ線診断用プログラムによれば、モアレパターンの第１高調波のピーク周波数に基づいてモアレパターンの第３高調波のピーク周波数を算出し、第３高調波を抽出するための第３高調波抽出フィルタを作成する。次に、第３高調波抽出フィルタに基づいて第１・第２高調波除去画像から第３高調波を抽出する。そして、抽出された第３高調波を第１・第２高調波除去画像から減算して第３高調波除去画像を取得する。それにより、モアレパターン除去画像に残っているモアレパターンの第３高調波を除去することができる。
　本発明に係るＸ線診断装置およびＸ線診断用プログラムによれば、モアレパターンの第１高調波のピーク周波数に基づいてモアレパターンの第３高調波のピーク周波数を算出し、第３高調波を抽出するための第３高調波抽出フィルタを作成する。次に、第３高調波抽出フィルタに基づいて第１・第２高調波除去画像から第３高調波を抽出する。そして、抽出された第３高調波を第１・第２高調波除去画像から減算して第３高調波除去画像を取得する。それにより、モアレパターン除去画像に残っているモアレパターンの第３高調波を除去することができる。
実施例に係るＸ線診断装置の概略構成を示すブロック図である。 モアレパターンの説明に供する図である。 モアレパターンの第３高調波の抽出フィルタ作成の説明に供する図である。 モアレパターンの第３高調波の一例を示す図である。 （ａ）は縦ローパスフィルタ（ＬＰＦ）処理の説明に供する図であり、（ｂ）は縦ローパスフィルタ（ＬＰＦ）の一例を示す図である。 モアレパターンの第３高調波除去の説明に供する図である。 輝度変化量の算出の説明に供する図である。 平滑化処理の説明に供する図である。 出力選択の説明に供する図である。 実施例のＸ線診断装置の動作を示すフローチャートである。 従来のＸ線診断装置の概略構成を示すブロック図である。 Ｘ線検出器とグリッドとによるモアレパターンの周波数の説明に供する図である。 画像の１次元のフーリエ変換結果の一例を示す。 グリッド誤差とモアレパターンの高調波のピーク周波数との関係を示す図である。
　１　　…　Ｘ線診断装置
　３　　…　Ｘ線管
　４　　…　Ｘ線検出器
　５　　…　散乱線除去グリッド
　９　　…　画像処理部
　１０　…　主制御部
　２１　…　第１・第２高調波除去処理部
　２３　…　第３高調波除去処理部
　２５　…　ピーク周波数検出部
　２９　…　第３高調波抽出フィルタ作成部
　３１　…　第３高調波抽出部
　３３　…　縦ローパスフィルタ処理部（縦ＬＰＦ処理部）
　３５　…　第３高調波減算部
　３７　…　平坦度算出部
　３９　…　出力選択部
　４１　…　輝度変化量算出部
　４３　…　正規化処理部
　４５　…　平滑化処理部
　以下、図面を参照して本発明の実施例を説明する。図１は、実施例に係るＸ線診断装置の概略構成を示すブロック図であり、図２は、モアレパターンの説明に供する図である。図３は、モアレパターンの第３高調波の抽出フィルタ作成の説明に供する図であり、図４は、モアレパターンの第３高調波の一例を示す図である。図５（ａ）は、縦ローパスフィルタ（ＬＰＦ）処理の説明に供する図であり、図５（ｂ）は縦ローパスフィルタ（ＬＰＦ）の一例を示す図である。図６は、モアレパターンの第３高調波除去の説明に供する図であり、図７は、輝度変化量の算出の説明に供する図である。図８は平滑化処理の説明に供する図であり、図９は出力選択の説明に供する図である。
　図１を参照する。Ｘ線診断装置１は、被検体Ｍを載置する天板２と、被検体ＭにＸ線を照射するＸ線管３と、被検体Ｍを挟んでＸ線管３と対向して配置され、被検体Ｍを透過したＸ線を検出するＸ線検出器４と、Ｘ線検出器４のＸ線入射側に配置され、散乱したＸ線を除去するグリッド５とを備えている。なお、Ｘ線管３は本発明におけるＸ線照射部に相当する。
　Ｘ線管３は、Ｘ線管制御部６によりＸ線照射に必要な制御が実行される。Ｘ線管制御部６は、Ｘ線管３の管電圧や管電流を発生させる高電圧発生部７を有している。Ｘ線管制御部６は、後述する入力部１２で設定された管電圧や管電流や照射時間等のＸ線条件に応じてＸ線管３からＸ線を照射させている。
　Ｘ線検出器４は、Ｘ線管３から照射されたＸ線を検出し、Ｘ線強度分布に応じた信号であるＸ線検出信号を出力する。Ｘ線検出器４は、フラットパネル型Ｘ線検出器（ＦＰＤ）やイメージインテンシファイア等で構成される。Ｘ線検出器４は、例えば、０．１５umピクセルピッチで構成され、３．３３lp/mmの解像度を持つ。また、グリッド５は、例えば、５０line/cm（５lp/mm）のグリッド密度で構成される。
　Ｘ線検出器４の後段には、順番にＡ／Ｄ変換器８と画像処理部９とが設けられている。Ａ／Ｄ変換器８は、Ｘ線検出器４から出力されたアナログのＸ線検出信号をデジタル信号に変換する。画像処理部９は、Ｘ線検出器４から出力されたデジタルのＸ線検出信号に基づいた画像に対して種々の処理を行う。また、Ｘ線診断装置１は、この装置１の各構成を統括的に制御する主制御部１０と、画像処理部９で処理された画像を表示する表示部１１と、操作者が入力設定や各種操作を行う入力部１２と、処理された画像を記憶する記憶部１３とを備えている。
　主制御部１０は、中央演算処理装置（ＣＰＵ）などで構成され、各種プログラムを実行するようになっている。主制御部１０は、例えば、被検体Ｍが載置された天板２、Ｘ線管３またはＸ線検出器４を所定の位置に移動させる制御を行う。表示部１１は、モニタ等で構成される。入力部１２は、キーボードやマウス等で構成される。記憶部１３は、ＲＯＭ（Read-only Memory）、ＲＡＭ（Random-Access Memory）またはハードディスク等の記憶媒体で構成される。
　画像処理部９は、Ｘ線検出器４の出力に基づいて取得した画像に映りこんだグリッド５のモアレパターンを除去するための第１・第２高調波除去処理部２１と第３高調波除去処理部２３とを備えている。第１・第２高調波除去処理部２１は、画像に映りこんだグリッド５のモアレパターンの第１高調波のピーク周波数を検出するピーク周波数検出部２５を備えている。第１・第２高調波除去処理部２１は、第１高調波のピーク周波数に基づいて画像からモアレパターンの第１高調波と第２高調波を除去して第１・第２高調波除去画像を取得する。
　第１・第２高調波除去処理部２１は、画像に映りこんだグリッド５のモアレパターンの第１高調波と第２高調波を除去する。第１・第２高調波除去処理部２１は、図２に示すように、まず、モアレパターン（符号ｍｐ）が映りこんだ画像に対して、横方向（モアレパターンの縞目と直交する方向）の画素列（符号Ｌ１）に１次元のフーリエ変換を行う。１次元のフーリエ変換は、例えば画像の上端の画素列から順番に行う。ピーク周波数検出部２５は、そのフーリエ変換結果からモアレパターンの第１高調波を示すピーク周波数を検出する（図１３）。なお、ピーク周波数検出部２５は、画像の上端から順番に検出された複数のピーク周波数から設定された所定のピーク周波数を検出する。第２高調波は、第１高調波の２倍（整数倍）であるので、第１高調波のピーク周波数により算出される。取得された第１高調波と第２高調波のピーク周波数に基づき、モアレパターンの第１高調波と第２高調波を抽出するためのフィルタを作成する。このフィルタを用いて第１高調波と第２高調波を抽出する。そして、抽出された第１高調波と第２高調波をモアレパターンが映りこんだ画像から減算することで、画像からモアレパターンの第１高調波と第２高調波を除去した第１・第２高調波除去画像を取得する。なお、図２において、符号ｆは画像内における被検体Ｍの平坦領域を示す。
　第３高調波除去処理部２３は、図１に示すように、記憶部２７、第３高調波抽出フィルタ作成部２９、第３高調波抽出部３１、縦ローパスフィルタ処理部（縦ＬＰＦ処理部）３３、第３高調波減算部３５、平坦度算出部３７および出力選択部３９を備えている。第１・第２高調波除去処理部２１によってモアレパターンの第１高調波と第２高調波が除去された第１・第２高調波除去画像は、記憶部２７に記憶される。
　第３高調波抽出フィルタ作成部２９は、第１高調波のピーク周波数に基づいてモアレパターンの第３高調波のピーク周波数を算出し、第３高調波を抽出するための第３高調波抽出フィルタを作成する。第１高調波のピーク周波数は、ピーク周波数検出部２５で検出された第１高調波のピーク周波数が与えられる。モアレパターンの第３高調波は、モアレパターンの第１高調波の３倍であるが、第３高調波がナイキスト周波数Ｎｙの３．３３lp/mmより大きい場合は、３．３３lp/mmの位置で折り返された位置が第３高調波のピーク周波数位置となる。すなわち、第１高調波のピーク周波数が１．７２lp/mmの場合、第３高調波のピーク周波数は１．７２×３＝５．１６lp/mmとなる。このとき第３高調波がナイキスト周波数Ｎｙの３．３３lp/mmより大きいので、３．３３lp/mmの位置で折り返される。第３高調波のピーク周波数は、３．３３－（５．１６－３．３３）＝１．５０lp/mmとなる。なお、第１高調波のピーク周波数が１．７２lp/mmの場合とは、図１４においてグリッド誤差（比）が０．９９のときである。第３高調波抽出フィルタは、図３に示すように、算出された第３高調波のピーク周波数に所定の帯域を持たせて作成される。
　第３高調波抽出部３１は、第３高調波抽出フィルタに基づいて第１・第２高調波除去画像から第３高調波を抽出する。第１・第２高調波除去画像は、記憶部２７により与えられる。第３高調波の抽出は、例えば１・第２高調波除去画像の横方向の画素列に対して、第３高調波抽出フィルタを用いたフィルタ処理を施す。この抽出処理を第１・第２高調波除去画像の例えば上端から順番に１列ずつ行う。これにより、図４に示すようなモアレパターンの第３高調波の画像が抽出される。
　しかしながら、抽出された第３高調波の画像には、図５（ａ）の符号ｎに示すように、モアレパターンの第３高調波の他に斜め方向のパターンが抽出される。これは、第３高調波抽出フィルタの周波数と同じ周波数の第１・第２高調波除去画像内のパターンが抽出されることが原因で起こる。すなわち、横方向の画素列Ｌ２において、周期がｄ１のモアレパターンの第３高調波が抽出されるとともに、第１・第２高調波除去画像内の同じ周期（ｄ２＝ｄ１）のパターンが抽出される。
　そこで、縦ＬＰＦ処理部３３は、第３高調波抽出部３１で抽出された第３高調波の縞目に沿った縦方向に１次元のローパスフィルタ（ＬＰＦ）処理を行う。モアレパターンの第３高調波の縞目を図５（ａ）の符号ｍｐ３に示す。すなわち、第３高調波の縞目に沿った縦方向にＬＰＦ処理を行う。縦方向のＬＰＦ処理は、抽出された第３高調波の画像の各画素に対して、図５（ｂ）に示すフィルタでフィルタ処理を行う。図５（ｂ）に示すフィルタは、縦方向の画素列であり、例えば、処理対象の注目画素の上下８画素を含む１７画素を用いて行う。この１７画素の処理前の画素値の加算平均を算出し、算出した加算平均値を注目画素の新たな画素値とする。この処理を各画素で行うことにより、画像の縦方向の画素が平滑化されて図５（ａ）に示す符号ｎの斜め方向のパターンを除去される。そして、図４に示すモアレパターンの第３高調波の画像を取得される。なお、縦ＬＰＦ処理部３３は本発明におけるローパスフィルタ処理部に相当する。なお、図５（ａ）において、モアレパターンの第３高調波の縞目が横方向に並んでいるが、縞目が縦方向に並ぶ場合、縞目に沿った横方向に１次元のローパスフィルタ（ＬＰＦ）処理を行う。
　第３高調波減算部３５は、抽出された第３高調波を第１・第２高調波除去画像から減算して第３高調波除去画像を取得する。すなわち、第３高調波減算部３５は、図６に示すように、第１・第２高調波除去画像Ｐ１から抽出された第３高調波の画像Ｐ２を減算する。これにより、第３高調波除去画像Ｐ３が取得される。
　平坦度算出部３７は、図１に示すように、第１・第２高調波除去画像内の平坦度を算出して平坦度情報を取得する。この平坦度算出部３７による平坦度情報の算出は、第３高調波を抽出して除去する第３高調波除去画像を取得する処理と並行して行われる。平坦度情報は、第１・第２高調波除去画像の平坦な領域（平坦領域）と平坦ではない構造物の領域（平坦領域以外の領域）識別するための情報である。例えば、被検体Ｍの脚部の撮影においては、平坦領域としては肉部分、平坦領域以外の領域としては骨部分が挙げられる。また、平坦度算出部３７は、輝度変化量算出部４１、正規化処理部４３および平滑化処理部４５を備えている。
　輝度変化量算出部４１は、処理対象である注目画素とその周囲の画素との輝度変化量を算出して平坦度情報を取得する。具体的には、輝度変化量算出部４１は、注目画素の画素値から注目画素とその周囲の画素との加算平均値を減算して平坦度情報を取得する。すなわち、輝度変化量算出部４１は、図７に示すように、注目画素とその注目画素と例えば上下左右１０画素との２１×２１画素で加算平均を行い、注目画素の画素値からその加算平均値を減算する。すなわち、例えば２１×２１画素の画素値の加算平均値をＡとし、注目画素の画素値をＢとすると輝度変化量Ｃ１は、Ｃ１＝｜Ｂ－Ａ｜で表すことができる。この処理を第１・第２高調波除去画像の各画素で行う。このように輝度変化量Ｃ１で表された新たな画像を平坦度情報として取得する。
　正規化処理部４３は、平坦度情報の各画素の輝度変化量をそれぞれ対応する第１・第２高調波除去画像の画素値で除算して正規化処理する。すなわち、正規化処理後の輝度変化量Ｃ２は、Ｃ２＝｜Ｂ－Ａ｜／Ｂで表される。画素値が高い画素では輝度変化量（画素値の変動量）が大きく、画素値が低い画素では輝度変化量が小さい傾向にある。また、例えば、画素値が１０００の注目画素と画素値が１０の注目画素があるものとする。両方の輝度変化量が２である場合、輝度変化量が同じ値である。しかしながら、画素値が１０００のときにおける輝度変化量２と、画素値が１０のときにおける輝度変化量２とでは、２／１０００＝０．００２と、２／１０＝０．２と、重みが異なる。そのため、高輝度画素と低輝度画素との輝度変化量を同じ基準で比較できるようにしている。
　平滑化処理部４５は、平坦度情報を複数個の画素で分割した領域ごとに輝度変化量の加算平均を行い、その加算平均値を分割した領域の各画素に適用させて平滑化する。平滑化処理部４５は、図８に示すように、例えば平坦度情報Ｐ４が６４×６４画素で構成される場合、例えば１６×１６画素で加算平均を行って１画素に縮小する。すると、４×４画素の縮小された平坦度情報Ｐ５が取得される。そして、４×４画素の平坦度情報Ｐ５を元の６４×６４画素の大きさに戻すこと（拡大）で平滑化処理された平坦度情報Ｐ６を得る。なお、平坦度情報Ｐ５の１画素から平坦度情報Ｐ６の１６×１６画素の大きさに戻す際に、平坦度情報Ｐ６の１６×１６画素の全ての輝度変化量は、平坦度情報Ｐ５の１画素の輝度変化量が与えられる。
　出力選択部３９は、平坦度情報の平坦領域では第３高調波除去画像を出力し、平坦度情報の平坦領域以外の領域では第１・第２高調波除去画像を出力する。すなわち、出力選択部３９は、図９に示すように、平坦度情報Ｐ６に基づいて、第１・第２高調波除去画像Ｐ１および第３高調波除去画像Ｐ３のいずれか一方の画像を出力するのかを選択する。平坦度情報Ｐ６における正規化処理および平滑化処理された輝度変化量がしきい値以上であれば、平坦領域と識別され、しきい値未満であれば、平坦領域以外の領域と識別される。平滑化処理部４５により１６×１６画素で加算平均を行った分割した領域の各画素（１６×１６画素）は、同じ輝度変化量であるので、その１６×１６画素の領域単位で出力選択部３９による選択出力をするようになっている。図９に示すように、平坦度情報Ｐ６の平坦領域（符号ｆ）では、第３高調波除去画像Ｐ３が出力される。一方、平坦度情報Ｐ６の平坦領域以外の領域では、第１・第２高調波除去画像Ｐ１が出力される。
　次に図１０を参照して、本実施例のＸ線診断装置１の動作をフローチャートに沿って説明する。
　〔ステップＳ０１〕Ｘ線撮影
　Ｘ線診断装置１は、図１に示すように、Ｘ線管制御部６で設定された管電圧や管電流や照射時間等のＸ線条件に応じてＸ線管３からＸ線を照射させる。Ｘ線管３から照射されたＸ線は、被検体Ｍを透過してＸ線検出器４に入射する。Ｘ線検出器４は、Ｘ線強度分布に応じたＸ線検出信号を出力する。Ａ／Ｄ変換器８は、Ｘ線検出器４から出力されたアナログのＸ線検出信号をデジタル信号に変換する。このデジタル変換されたＸ線検出信号に基づいて被検体Ｍを撮影した画像を取得することができる。
　〔ステップＳ０２〕第１・第２高調波除去処理
　ピーク周波数検出部２５は、画像に映りこんだ前記グリッドのモアレパターンの第１高調波のピーク周波数を検出する。第１・第２高調波除去処理部２１は、第１高調波のピーク周波数に基づいて画像からモアレパターンの第１高調波と第２高調波を除去して第１・第２高調波除去画像を取得する。
　〔ステップＳ０３〕第３高調波抽出フィルタの作成
　第３高調波抽出フィルタ作成部２９は、第１高調波のピーク周波数に基づいてモアレパターンの第３高調波のピーク周波数を算出し、第３高調波を抽出するための第３高調波抽出フィルタを作成する（図３）。
　〔ステップＳ０４〕第３高調波の抽出
　第３高調波抽出部３１は、第３高調波抽出フィルタに基づいて第１・第２高調波除去画像から第３高調波を抽出する（図５（ａ））。第３高調波の抽出は、例えば、１・第２高調波除去画像の横方向の画素列の上端から順番に１列ずつに第３高調波抽出フィルタを用いたフィルタ処理を施して行われる。
　〔ステップＳ０５〕縦ＬＰＦ処理
　縦ＬＰＦ処理部３３は、第３高調波抽出部３１で抽出された第３高調波（図５（ａ））の縞目に沿った縦方向に１次元のローパスフィルタ（ＬＰＦ）処理を行う。縦方向のＬＰＦ処理は、例えば、図５（ｂ）に示すフィルタ、すなわち、処理対象の注目画素の上下８画素を含む１７画素を用い、この１７画素の処理前の画素値の加算平均値を注目画素の新たな画素値とする。この処理を各画素で行うことにより、画像の縦方向の画素が平滑化されて図５（ａ）に示す符号ｎの斜め方向のパターンを除去される。そのため、図４に示すようなモアレパターンの第３高調波の画像を取得することができる。
　〔ステップＳ０６〕第３高調波除去処理（減算）
　第３高調波減算部３５は、図６に示すように、抽出された第３高調波Ｐ２を第１・第２高調波除去画像Ｐ１から減算して第３高調波除去画像Ｐ３を取得する。
　また、第３高調波除去画像を取得するステップＳ０３～Ｓ０６と並行して、第１・第２高調波除去画像内の平坦度を算出して平坦度情報を取得するステップＳ１３～Ｓ１５を行う。
　〔ステップＳ１３〕輝度変化量の算出
　輝度変化量算出部４１は、注目画素の画素値から注目画素とその周囲の画素との加算平均値を減算して平坦度情報を取得する。すなわち、輝度変化量算出部４１は、図７に示すように、注目画素を含む２１×２１画素の画素値の加算平均値をＡとし、注目画素の画素値をＢとして、輝度変化量Ｃ１＝｜Ｂ－Ａ｜を算出する。この処理を各画素で行うことにより、新たな画像として平坦度情報として取得する。
　〔ステップＳ１４〕正規化処理
　正規化処理部４３は、平坦度情報の各画素の輝度変化量をそれぞれ対応する第１・第２高調波除去画像の画素値で除算して正規化処理する。すなわち、輝度変化量Ｃ１の算出で用いた注目画素の画素値Ｂで輝度変化量Ｃ１を除算する。正規化処理後の輝度変化量Ｃ２は、Ｃ２＝｜Ｂ－Ａ｜／Ｂで算出される。
　〔ステップＳ１５〕平滑化処理
　平滑化処理部４５は、平坦度情報を複数個の画素で分割した領域ごとに輝度変化量の加算平均を行い、その加算平均値を分割した領域の各画素に適用させて平滑化する。すなわち、平滑化処理部４５は、図８に示すように、例えば１６×１６画素で加算平均を行って１画素に縮小する。そして、元の大きさに戻す際に、元の大きさの１６×１６画素には、全画素に加算平均された輝度変化量が与えられる。
　〔ステップＳ０７〕出力選択
　出力選択部３９は、図９に示すように、平坦度情報Ｐ６の平坦領域ｆでは第３高調波除去画像Ｐ３を出力し、平坦度情報Ｐ６の平坦領域ｆ以外の領域では第１・第２高調波除去画像Ｐ１を出力する。このように平坦領域ｆではモアレパターンの第３高調波が除去された出力選択画像Ｐ７が取得される。
　〔ステップＳ０８〕表示・保存
　ステップＳ０７で出力選択された出力選択画像Ｐ７は、画像処理部９でその他必要な処理が行われる。そして、出力選択画像Ｐ７は、表示部１１に表示され、また、記憶部１３に記憶される。
　本実施例に係るＸ線診断装置１によれば、第３高調波抽出フィルタ作成部２９は、第１高調波のピーク周波数に基づいてモアレパターンの第３高調波のピーク周波数を算出し、第３高調波を抽出するための第３高調波抽出フィルタを作成する。第３高調波抽出部３１は、第３高調波抽出フィルタに基づいて第１・第２高調波除去画像から第３高調波を抽出する。第３高調波減算部３５は、抽出された第３高調波を第１・第２高調波除去画像から減算して第３高調波除去画像を取得する。それにより、第１・第２高調波除去画像に残っているモアレパターンの第３高調波を除去することができる。
　また、Ｘ線診断装置１は、第１・第２高調波除去画像内の平坦度を算出して平坦度情報を取得する平坦度算出部３７と、平坦度情報の平坦領域では第３高調波除去画像を出力し、平坦度情報の平坦領域以外の領域では第１・第２高調波除去画像を出力する出力選択部３９とを備えている。モアレパターンの第３高調波は、スペクトル強度が小さい周波数成分であるので、画像内では、輝度変化の少ない平坦領域に目立ち、その平坦領域以外では認識されない。そのため、出力選択部３９は、平坦度算出部３７で平坦度を算出して取得した平坦度情報のうち、平坦領域では第３高調波除去画像を出力し、平坦度情報の平坦以外の領域では第１・第２高調波除去画像を出力する。すなわち、第１・第２高調波除去画像にてモアレパターンの第３高調波が目立つ領域は平坦領域であるので、平坦領域とそれ以外の領域とを識別して、平坦領域のみに第３高調波除去画像を選択して出力させている。それにより、画像の視認性の劣化を抑制するとともに第１・第２高調波除去画像に残っていた第３高調波のモアレパターンを除去することができる。
　すなわち、平坦領域以外では視認性を劣化させず、モアレパターン残りのない画像を取得することができる。なお、平坦領域では輝度変化が少なくモアレパターンの第３高調波を除去しても視認性劣化の影響を比較的に受けにくい。また、平坦領域のモアレパターンの第３高調波を除去しているので、平坦領域のＳ／Ｎ比を改善することも可能である。
　また、輝度変化量算出部４１は、第１・第２高調波除去画像内の平坦度を示す平坦度情報を、注目画素とその周囲の画素との輝度変化量を算出することで取得することができる。さらに、輝度変化量算出部４１は、平坦度情報を、注目画素の画素値から注目画素とその周囲の画素との加算平均値を減算して輝度変化量を算出することで取得することができる。
　また、正規化処理部４３は、平坦度情報の各画素の輝度変化量をそれぞれ対応する第１・第２高調波除去画像の画素値で除算して正規化処理する。すなわち、画素値が高い画素では輝度変化量が大きく、画素値が低い画素では輝度変化量が小さい傾向にある。輝度変化量が同じであっても高輝度画素の輝度変化量と低輝度画素の輝度変化量とでは、その重みが異なる。そのため、各画素の輝度変化量をそれぞれ対応する画素の画素値で除算することにより、高輝度画素と低輝度画素との輝度変化量を同じ基準で比較することができる。
　また、平滑化処理部４５は、平坦度情報を複数個の画素（例えば１６×１６画素）で分割した領域ごとに輝度変化量の加算平均を行い、その加算平均値を分割した領域の各画素（１６×１６画素）に適用させて平滑化している。出力選択部３９による画像の出力選択を画素単位で行うのではなく、平滑化により任意の領域で行うことができる。また、平坦領域以外の領域内に点在する平坦領域と判定された領域を平坦領域以外の領域として処理することができる。平坦領域以外の領域内に点在する平坦領域まで第３高調波除去画像を出力させると視認性劣化を招くおそれがある。
　また、縦ＬＰＦ処理部３３は、第３高調波抽出部３１で抽出された第３高調波の縞目に沿った縦方向に１次元のＬＰＦ処理を行う。モアレパターンの第３高調波の他に抽出される第３高調波の縞目に対して斜め方向のパターンを除去することができる。
　本発明は、上記実施形態に限られることはなく、下記のように変形実施することができる。
　（１）上述した実施例において、画像処理部９は、プログラムを実行させるためのＣＰＵ等で構成される制御部と、プログラム等を記憶するＲＯＭやＲＡＭ等の記憶媒体で構成される記憶部とを備えてもよい。記憶部には、各ステップＳ０１～Ｓ０８，Ｓ１３～Ｓ１５の動作のプログラムを記憶させて、そのプログラムを制御部で実行するようにしてもよい。この場合、このプログラムに必要な操作は、入力部１２で入力され、出力選択画像Ｐ７は表示部１１に表示される。
　（２）上述した実施例において、画像処理部９に限定されず、各ステップＳ０１～Ｓ０８，Ｓ１３～Ｓ１５の動作を、その動作のプログラムを記憶部１３に記憶し、主制御部１０で実行するようにしてもよい。この場合、このプログラムに必要な操作は入力部１２で入力され、出力選択画像Ｐ７は表示部１１に表示される。また、その動作のプログラムは、ＬＡＮ等のネットワークシステムでＸ線診断装置１と接続されたパソコン上でも実行できるようにしてもよい。

Claims (8)

  1. 　被検体にＸ線を照射するＸ線照射部と、
    　被検体を透過したＸ線を検出するＸ線検出器と、
    　前記Ｘ線検出器のＸ線入射側に配置され、散乱したＸ線を除去するグリッドとを備え、
    　前記Ｘ線検出器の出力に基づいて画像を取得するＸ線診断装置において、
    　前記画像に映りこんだ前記グリッドのモアレパターンの第１高調波のピーク周波数を検出するピーク周波数検出部と、
    　前記第１高調波のピーク周波数に基づいて前記画像からモアレパターンの第１高調波と第２高調波を除去して第１・第２高調波除去画像を取得する第１・第２高調波除去処理部と、
    　前記第１高調波のピーク周波数に基づいてモアレパターンの第３高調波のピーク周波数を算出し、前記第３高調波を抽出するための第３高調波抽出フィルタを作成する第３高調波抽出フィルタ作成部と、
    　前記第３高調波抽出フィルタに基づいて前記第１・第２高調波除去画像から前記第３高調波を抽出する第３高調波抽出部と、
    　抽出された前記第３高調波を前記第１・第２高調波除去画像から減算して第３高調波除去画像を取得する第３高調波減算部と、
    　を備えていることを特徴とするＸ線診断装置。
  2. 　請求項１に記載のＸ線診断装置において、
    　前記第１・第２高調波除去画像内の平坦度を算出して平坦度情報を取得する平坦度算出部と、
    　前記平坦度情報の平坦領域では第３高調波除去画像を出力し、前記平坦度情報の平坦領域以外の領域では前記第１・第２高調波除去画像を出力する出力選択部とを備えていることを特徴とするＸ線診断装置。
  3. 　請求項２に記載のＸ線診断装置において、
    　前記平坦度算出部は、処理対象である注目画素とその周囲の画素との輝度変化量を算出して平坦度情報を取得する輝度変化量算出部を備えていることを特徴とするＸ線診断装置。
  4. 　請求項３に記載のＸ線診断装置において、
    　前記輝度変化量算出部は、前記注目画素の画素値から前記注目画素とその周囲の画素との加算平均値を減算して平坦度情報を取得することを特徴とするＸ線診断装置。
  5. 　請求項３または４に記載のＸ線診断装置において、
    　前記平坦度算出部は、前記平坦度情報の各画素の輝度変化量をそれぞれ対応する前記第１・第２高調波除去画像の画素値で除算して正規化処理する正規化処理部を備えていることを特徴とするＸ線診断装置。
  6. 　請求項３から５のいずれかに記載のＸ線診断装置において、
    　前記平坦度算出部は、前記平坦度情報を複数個の画素で分割した領域ごとに輝度変化量の加算平均を行い、その加算平均値を分割した領域の各画素に適用させて平滑化する平滑化処理部を備えていることを特徴とするＸ線診断装置。
  7. 　請求項１から６のいずれかに記載のＸ線診断装置において、
    　前記第３高調波抽出部で抽出された前記第３高調波の縞目に沿った方向に１次元のローパスフィルタ処理を行うローパスフィルタ処理部を備えていることを特徴とするＸ線診断装置。
  8. 　散乱したＸ線を除去するグリッドを介して撮影された画像の処理工程をコンピュータに実行させるためのＸ線診断用プログラムであって、
    　前記画像に映りこんだ前記グリッドのモアレパターンの第１高調波のピーク周波数を検出する工程と、
    　前記第１高調波のピーク周波数に基づいて前記画像からモアレパターンの第１高調波と第２高調波を除去して第１・第２高調波除去画像を取得する工程と、
    　前記第１高調波のピーク周波数に基づいてモアレパターンの第３高調波のピーク周波数を算出し、前記第３高調波を抽出するための第３高調波抽出フィルタを作成する工程と、
    　前記第３高調波抽出フィルタに基づいて前記第１・第２高調波除去画像から前記第３高調波を抽出する工程と、
    　抽出された前記第３高調波を前記第１・第２高調波除去画像から減算して第３高調波除去画像を取得する工程と、を備え、
    　これらの工程をコンピュータに実行させることを特徴とするＸ線診断用プログラム。
PCT/JP2011/002377 2011-04-22 2011-04-22 Ｘ線診断装置およびｘ線診断用プログラム Ceased WO2012143979A1 (ja)

Priority Applications (5)

Application Number Priority Date Filing Date Title
EP11864008.5A EP2700359B1 (en) 2011-04-22 2011-04-22 X-ray diagnostic apparatus and x-ray diagnostic program
PCT/JP2011/002377 WO2012143979A1 (ja) 2011-04-22 2011-04-22 Ｘ線診断装置およびｘ線診断用プログラム
JP2013510742A JP5527481B2 (ja) 2011-04-22 2011-04-22 Ｘ線診断装置およびｘ線診断用プログラム
US14/113,176 US9232925B2 (en) 2011-04-22 2011-04-22 X-ray diagnostic apparatus and storage medium storing X-ray diagnostic program
CN201180070344.4A CN103491875B (zh) 2011-04-22 2011-04-22 X射线诊断装置以及莫尔条纹去除方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/JP2011/002377 WO2012143979A1 (ja) 2011-04-22 2011-04-22 Ｘ線診断装置およびｘ線診断用プログラム

Publications (1)

Publication Number Publication Date
WO2012143979A1 true WO2012143979A1 (ja) 2012-10-26

Family

ID=47041134

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2011/002377 Ceased WO2012143979A1 (ja) 2011-04-22 2011-04-22 Ｘ線診断装置およびｘ線診断用プログラム

Country Status (5)

Country Link
US (1) US9232925B2 (ja)
EP (1) EP2700359B1 (ja)
JP (1) JP5527481B2 (ja)
CN (1) CN103491875B (ja)
WO (1) WO2012143979A1 (ja)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2015208638A (ja) * 2014-04-30 2015-11-24 株式会社東芝 Ｘ線診断装置及び画像処理装置
WO2017138097A1 (ja) * 2016-02-09 2017-08-17 株式会社島津製作所 Ｘ線撮影装置
JP7551703B2 (ja) 2022-07-20 2024-09-17 キヤノン株式会社 画像処理装置、画像処理方法、およびプログラム

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP6169626B2 (ja) * 2014-03-10 2017-07-26 富士フイルム株式会社 放射線画像処理装置、方法およびプログラム
CN105832359B (zh) * 2016-03-22 2018-10-19 广州七喜医疗设备有限公司 一种自适应x线滤线栅栅影去除的方法
JP2024025163A (ja) * 2022-08-10 2024-02-26 キヤノン株式会社 放射線撮影装置及びその制御方法、放射線撮影システム、プログラム

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2002152467A (ja) 2000-08-28 2002-05-24 Fuji Photo Film Co Ltd 画像信号生成方法および装置並びにプログラム
JP2003319260A (ja) * 2002-04-25 2003-11-07 Fuji Photo Film Co Ltd 画像信号生成方法および装置
JP2005052553A (ja) * 2003-08-07 2005-03-03 Canon Inc 放射線画像処理方法及び装置並びにグリッドの選別方法及び装置
JP2010240259A (ja) * 2009-04-08 2010-10-28 Shimadzu Corp 放射線撮像装置
JP2011010873A (ja) * 2009-07-02 2011-01-20 Shimadzu Corp 放射線撮影装置

Family Cites Families (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
DE69522844T2 (de) * 1994-02-16 2002-05-16 Koninklijke Philips Electronics N.V., Eindhoven Bildverarbeitungsverfahren zur örtlichen Bestimmung des Zentrums und der Objekte, die sich gegen einen Hintergrund abheben, und Vorrichtung zur Durchführung des Verfahrens
EP1202555A3 (en) 2000-08-28 2004-12-29 Fuji Photo Film Co., Ltd. Image signal generating method, apparatus and program
JP4118535B2 (ja) * 2001-07-03 2008-07-16 株式会社日立メディコ Ｘ線検査装置
JP2003233818A (ja) * 2002-02-06 2003-08-22 Fuji Photo Film Co Ltd 周期的パターン判定方法および装置並びにプログラム
US7796792B2 (en) * 2005-06-29 2010-09-14 Agfa Healthcare, N.V. Method of identifying disturbing frequencies originating from the presence of an anti-scatter grid during acquisition of a radiation image
CN100510725C (zh) * 2006-11-14 2009-07-08 北京国药恒瑞美联信息技术有限公司 用于消除散射辐射影响的虚拟滤线栅成像方法及其系统
CN101743565B (zh) * 2007-03-06 2012-08-08 皇家飞利浦电子股份有限公司 图像序列的滤波
EP2196148B1 (en) * 2007-10-02 2015-03-04 Shimadzu Corporation Radiological image processing apparatus and radiological image processing program

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2002152467A (ja) 2000-08-28 2002-05-24 Fuji Photo Film Co Ltd 画像信号生成方法および装置並びにプログラム
JP2003319260A (ja) * 2002-04-25 2003-11-07 Fuji Photo Film Co Ltd 画像信号生成方法および装置
JP2005052553A (ja) * 2003-08-07 2005-03-03 Canon Inc 放射線画像処理方法及び装置並びにグリッドの選別方法及び装置
JP2010240259A (ja) * 2009-04-08 2010-10-28 Shimadzu Corp 放射線撮像装置
JP2011010873A (ja) * 2009-07-02 2011-01-20 Shimadzu Corp 放射線撮影装置

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
See also references of EP2700359A4

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2015208638A (ja) * 2014-04-30 2015-11-24 株式会社東芝 Ｘ線診断装置及び画像処理装置
WO2017138097A1 (ja) * 2016-02-09 2017-08-17 株式会社島津製作所 Ｘ線撮影装置
JP7551703B2 (ja) 2022-07-20 2024-09-17 キヤノン株式会社 画像処理装置、画像処理方法、およびプログラム

Also Published As

Publication number Publication date
US20140086472A1 (en) 2014-03-27
US9232925B2 (en) 2016-01-12
EP2700359A1 (en) 2014-02-26
JP5527481B2 (ja) 2014-06-18
JPWO2012143979A1 (ja) 2014-07-28
EP2700359B1 (en) 2015-07-15
CN103491875A (zh) 2014-01-01
CN103491875B (zh) 2015-09-30
EP2700359A4 (en) 2014-09-03

Similar Documents

Publication Publication Date Title
JP5527481B2 (ja) Ｘ線診断装置およびｘ線診断用プログラム
JP6105586B2 (ja) エネルギー高感度検出の微分位相コントラストイメージング
US9820713B2 (en) Image processing apparatus, radiation imaging system, control method, and storage medium
JP5848697B2 (ja) Ｘ線画像診断装置
JP6215011B2 (ja) Ｘ線診断装置
JP2013102951A (ja) 撮像装置および画像処理方法
WO2013035255A1 (ja) 画像処理装置およびそれを備えた放射線撮影装置
CN107203983A (zh) 用于减少x射线图像中的栅线伪影的方法及系统
JP2011136102A (ja) 放射線撮像装置
JP2019122616A (ja) 画像処理装置、画像処理方法、及び画像処理プログラム
JP5699842B2 (ja) 画質評価方法およびそれを用いたｘ線透視撮影装置
JP5821790B2 (ja) Ｘ線診断装置
JP5360397B2 (ja) 放射線撮影装置
JP2018029777A (ja) Ｘ線位相差撮像装置
JP2014176565A (ja) 画像処理装置、放射線撮影装置、画像処理方法、コンピュータプログラム及び記憶媒体
JP2018038638A (ja) 画像処理装置、方法およびプログラム
US20200229785A1 (en) A device and method for scatter correction in an x-ray image
US10335109B2 (en) Radiation phase-contrast imaging device
JP5077587B2 (ja) 放射線撮像装置
Strauss et al. Image quality dependence on image processing software in computed radiography
US8897418B2 (en) X-ray apparatus
JP5407937B2 (ja) Ｘ線撮影装置
JP6677263B2 (ja) Ｘ線撮影装置
US8867699B2 (en) Radiographic device
JP2002330344A (ja) 放射線画像処理装置、画像処理システム、放射線画像処理方法、記憶媒体、及びプログラム

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

Country of ref document: EP

Kind code of ref document: A1

ENP Entry into the national phase

Ref document number: 2013510742

Country of ref document: JP

Kind code of ref document: A

WWE Wipo information: entry into national phase

Ref document number: 14113176

Country of ref document: US

NENP Non-entry into the national phase

Ref country code: DE

WWE Wipo information: entry into national phase

Ref document number: 2011864008

Country of ref document: EP