WO2013136784A1 - 画像処理装置および方法並びにプログラム - Google Patents

画像処理装置および方法並びにプログラム Download PDF

Info

Publication number
WO2013136784A1
WO2013136784A1 PCT/JP2013/001636 JP2013001636W WO2013136784A1 WO 2013136784 A1 WO2013136784 A1 WO 2013136784A1 JP 2013001636 W JP2013001636 W JP 2013001636W WO 2013136784 A1 WO2013136784 A1 WO 2013136784A1
Authority
WO
WIPO (PCT)
Prior art keywords
solid sphere
order partial
image
partial differential
filtering
Prior art date
Application number
PCT/JP2013/001636
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 CN201380014282.4A priority Critical patent/CN104168820B/zh
Priority to EP13760382.5A priority patent/EP2826415B1/en
Publication of WO2013136784A1 publication Critical patent/WO2013136784A1/ja
Priority to US14/484,046 priority patent/US9183637B2/en
Priority to IN8510DEN2014 priority patent/IN2014DN08510A/en

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/12—Edge-based segmentation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/0002—Inspection of images, e.g. flaw detection
    • G06T7/0012—Biomedical image inspection
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/11—Region-based segmentation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/149—Segmentation; Edge detection involving deformable models, e.g. active contour models
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/168—Segmentation; Edge detection involving transform domain methods
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B2576/00—Medical imaging apparatus involving image processing or analysis
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/48—Other medical applications
    • A61B5/4887—Locating particular structures in or on the body
    • A61B5/489—Blood vessels
    • 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/10—Image acquisition modality
    • G06T2207/10072—Tomographic images
    • G06T2207/10088—Magnetic resonance imaging [MRI]
    • 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/10104—Positron emission tomography [PET]
    • 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/20016—Hierarchical, coarse-to-fine, multiscale or multiresolution image processing; Pyramid transform
    • 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/20024—Filtering details
    • 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/20048—Transform domain processing
    • 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/20048—Transform domain processing
    • G06T2207/20056—Discrete and fast Fourier transform, [DFT, FFT]
    • 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/20072—Graph-based image processing
    • 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/20112—Image segmentation details
    • G06T2207/20116—Active contour; Active surface; Snakes
    • 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/20112—Image segmentation details
    • G06T2207/20152—Watershed segmentation
    • 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
    • G06T2207/30008—Bone
    • 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
    • G06T2207/30061—Lung
    • 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
    • G06T2207/30088—Skin; Dermal
    • 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
    • G06T2207/30101—Blood vessel; Artery; Vein; Vascular
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16H—HEALTHCARE INFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR THE HANDLING OR PROCESSING OF MEDICAL OR HEALTHCARE DATA
    • G16H30/00—ICT specially adapted for the handling or processing of medical images
    • G16H30/40—ICT specially adapted for the handling or processing of medical images for processing medical images, e.g. editing

Definitions

  • the present invention relates to an image processing apparatus and method for determining the structure of a structure in an image, and a program for causing a computer to execute the image processing method.
  • Hessian analysis using a Hessian matrix has been proposed as a technique for extracting blood vessels and bones in medical images (see Non-Patent Document 1).
  • the Hessian analysis is performed by analyzing the eigenvalues of the Hessian matrix having the second-order partial differential coefficient calculated by using a second-order differential kernel for a predetermined filter such as a Gaussian kernel. Whether a point, a line, or a surface is determined.
  • Hesse analysis it is possible to discriminate blood vessels as linear structures and bones as planar structures.
  • Non-Patent Document 1 is that when another structure (peripheral structure) is located around the linear structure, a part of the peripheral structure is mistakenly identified as a linear structure. There was a case.
  • the method of Non-Patent Document 2 uses the filter proposed in Non-Patent Document 1 as a second-order partial differential of a function (solid sphere model function) representing a sphere (a solid sphere with 1 on the inside and 0 on the outside).
  • a function solid sphere model function representing a sphere (a solid sphere with 1 on the inside and 0 on the outside).
  • An object of the present invention is to prevent an erroneous determination of a structure that occurs in a contour portion of a structure having a radius of curvature larger than that of a solid sphere represented by a solid sphere model function.
  • the image processing apparatus includes a filtering unit that calculates a Hessian matrix by performing filtering using second-order partial differentiation of a function representing a solid sphere at each pixel position in an image, and a calculated Hessian.
  • An evaluator that discriminates structures included in the image using eigenvalues and eigenvectors obtained by eigenvalue analysis of the matrix, and the filtering unit includes the radius of the solid sphere at each pixel position in the image.
  • a correction unit that corrects one of the response waveforms appearing at two positions symmetrically spaced from the center of the solid sphere of the second-order partial differential in each direction of the function representing the real sphere is provided. It is characterized in that those were.
  • the image processing method includes a filtering step of calculating a Hessian matrix by performing filtering using second-order partial differentiation of a function representing a solid sphere at each pixel position in an image, and a calculated Hessian.
  • An image processing method comprising: an evaluation step for discriminating structures included in an image using eigenvalues and eigenvectors obtained by eigenvalue analysis of the matrix, wherein the filtering step is performed at each pixel position in the image.
  • the first partial differential vector is obtained by filtering using the first partial partial differential of the function representing the hollow sphere having the same radius as the radius of the solid sphere, and the obtained first partial partial vector is projected in the direction of the eigenvector Using the obtained value, one of the response waveforms of the second-order partial differential in each direction of the function representing the solid sphere appearing at two positions symmetrically spaced from the center of the solid sphere is obtained. It performs a correction to erase Chi.
  • the image processing method according to the first invention may be provided as a program for causing a computer to execute the image processing method.
  • the above-mentioned “hollow sphere having the same radius as that of the solid sphere” is not limited to the case where the spheres of the solid sphere and the hollow sphere strictly match each other, but also uses the first-order partial differential vector of the function representing the hollow sphere In the range where the effect of canceling the response waveform at one position among the response waveforms appearing at two positions symmetrically spaced from the center of the solid sphere in the second order partial differential of the function representing the solid sphere is obtained. If so, the case where the radius of the hollow sphere is larger or smaller than the radius of the solid sphere is included.
  • the radius of the hollow sphere and the radius of the solid sphere are preferably as equal as possible.
  • the difference between the radius of the hollow sphere and the radius of the solid sphere Is preferably 20% or less, and the difference between the radius of the hollow sphere and the radius of the solid sphere is preferably 10% or less.
  • the eigenvalues of the evaluation matrix are ⁇ 1 , ⁇ 2 , ⁇ 3
  • the corrected eigenvalues are ⁇ 1 ′ , ⁇ 2 ′ , ⁇ 3 ′
  • e 2 (x 2 , y 2 , z 2 )
  • e 3 (x 3 , y 3 , z 3 )
  • the function representing the hollow sphere is represented by the following formula (10).
  • x, y, and z are coordinates in a three-dimensional space
  • r is a polar coordinate representation thereof
  • R 4 is a radius of the hollow sphere.
  • An image processing apparatus includes a filtering unit that calculates a Hessian matrix by performing filtering using second-order partial differentiation of a function representing a solid sphere at each pixel position in an image, and a calculated Hessian.
  • An evaluator that discriminates a structure included in the image using eigenvalues and eigenvectors obtained by eigenvalue analysis of the matrix, and the filtering unit includes a first solid sphere at each pixel position in the image.
  • the first partial differential vector is calculated by performing filtering using the first partial differential of the function representing the second solid sphere having the second radius larger than the first radius which is the radius of the image, and in the image
  • the first partial differential vector is calculated by performing filtering using the first partial partial differentiation of the function representing the third solid sphere having the third radius smaller than the first radius at each pixel position of Solid sphere
  • a function that represents a solid sphere by using a value obtained by projecting the first-order partial differential vector of the function to be represented in the direction of the eigenvector and a value obtained by projecting the first-order partial differential vector of the function that represents the third solid sphere in the direction of the eigenvector.
  • a correction unit that performs correction to cancel the response waveform at one of the response waveforms that appear at two positions that are symmetrically spaced from the center of the solid sphere.
  • the response waveform at one position is a waveform in which one positive peak and one negative peak are adjacent to each other, and the second radius is positive from the center of the solid sphere.
  • the length from the center of the solid sphere to the negative peak coincides with the longer one, and the third radius is from the center of the solid sphere to the positive peak.
  • the length and the length from the center of the solid sphere to the negative peak it is preferable to match the shorter one.
  • the “second radius is the longer of the length from the center of the first solid sphere to the positive peak and the length from the center of the first solid sphere to the negative peak. Is the same as the length of the longer radius from the center of the solid sphere to the positive peak and the length from the center of the solid sphere to the negative peak (hereinafter referred to as the first radius).
  • the second radius is larger or smaller than the first length. If also be included.
  • the second radius and the first length are as equal as possible.
  • the difference between the second radius and the first length Is preferably 20% or less, and more preferably 10% or less.
  • the third radius is longer among the length from the center of the first solid sphere to the positive peak and the length from the center of the first solid sphere to the negative peak. Is equal to the length of the direction "means that the third radius is the shorter of the length from the center of the solid sphere to the positive peak and the length from the center of the solid sphere to the negative peak ( In the following, not only the case where the second length) is exactly the same, but also the first order of the function representing the third solid sphere and the weighted vector of the first partial differential vector of the function representing the third solid sphere One of response waveforms appearing at two positions symmetrically spaced from the center of the solid sphere of the second-order partial differential of each direction of the function representing the first solid sphere using a vector weighted by the partial differential vector If the effect of canceling the response waveform at the position is obtained, the third radius is larger than the second length.
  • the third radius and the second length be as equal as possible.
  • the difference between the third radius and the second length Is preferably 20% or less, and more preferably 10% or less.
  • the filtering unit calculates, for each function representing a solid sphere of a plurality of sizes, an evaluation matrix obtained by performing filtering using a second-order partial differential matrix of the function representing each solid sphere. It is preferable that
  • the filtering unit subjects the function representing the first solid sphere of a plurality of sizes to filtering by a second-order partial differential matrix of the function representing each first solid sphere. It is preferable that each of the evaluated matrices is calculated.
  • the image is a medical image and the structure is a blood vessel.
  • the filtering unit performs filtering using a second-order partial differential matrix of a function representing a solid sphere in Fourier space.
  • the evaluation unit discriminates at least one of a local point-like structure, a linear structure, and a planar structure of the structure.
  • the filtering unit that calculates the Hessian matrix by performing filtering using the second-order partial differentiation of the function representing the solid sphere at each pixel position in the image, and the calculated Hessian matrix as the eigenvalue
  • An evaluator that discriminates structures included in the image using eigenvalues and eigenvectors obtained by analysis, and the filtering unit has a radius equal to the radius of the solid sphere at each pixel position in the image.
  • the first-order partial differential vector is obtained by performing filtering using the first-order partial differential of the function representing the hollow sphere having the solid sphere by using the value obtained by projecting the obtained first-order partial differential vector in the direction of the eigenvector.
  • a correction unit that corrects one of the response waveforms appearing at two positions symmetrically spaced from the center of the solid sphere of the second-order partial differentiation in each direction of the function to be expressed. For this reason, only the position of the contour portion of the structure having a curvature larger than the curvature of the solid sphere represented by the function representing the solid sphere and the response waveform of the second-order partial differential in each direction of the function representing the solid sphere. It is possible to suppress misjudgment that occurs when the two match, and to improve the accuracy of the evaluation value. For this reason, it is possible to more accurately determine the structure included in the image based on the evaluation value.
  • the filtering unit that calculates the Hessian matrix by performing filtering using the second-order partial differentiation of the function representing the solid sphere at each pixel position in the image, and the calculated Hessian matrix as the eigenvalue
  • An evaluator that discriminates a structure included in the image using the eigenvalue and eigenvector obtained by the analysis, and the filtering unit has a first radius larger than the first radius that is the radius of the first solid sphere.
  • the partial differential vector is calculated, and the value obtained by projecting the first partial differential vector of the function representing the second solid sphere in the direction of the eigenvector and the first partial partial vector of the function representing the third solid sphere are the direction of the eigenvector.
  • a correction unit is provided that corrects the second-order partial differential of the function representing the solid sphere by canceling the response waveform at one of the response waveforms appearing at two positions symmetrically spaced from the center of the solid sphere. ing.
  • FIG. 1 is a schematic block diagram showing the configuration of an image processing apparatus according to a first embodiment of the present invention.
  • Diagram for explaining multi-resolution conversion Diagram for explaining eigenvalues of linear structure Diagram for explaining eigenvalues of planar structure
  • the figure for demonstrating the response of the 2nd-order partial differential of the x direction of the solid sphere model function used for the filtering of the 1st Embodiment of this invention The figure for demonstrating the response of the 1st-order partial differentiation of the x direction of the hollow sphere model function used for the correction process of the 1st Embodiment of this invention.
  • the figure for demonstrating the response of the 2nd-order partial differential of the x direction of the solid sphere corrected by the 1st-order partial differential vector of the hollow sphere model function of the first embodiment of the present invention The figure explaining the response which performed the correction process by the 1st Embodiment of this invention (when the size of a structure and the size of the solid sphere represented by the solid sphere model function correspond) The figure explaining the response which implemented the correction process by the 1st Embodiment of this invention (when the size of the solid sphere which a solid sphere model function represents with respect to the size of a structure is large) The flowchart which shows the process performed in the 1st Embodiment of this invention.
  • Example of applying image processing according to first embodiment of the present invention to blood vessel extraction (pseudo three-dimensional image)
  • Example (tomographic image) in which image processing according to the first embodiment of the present invention is applied to blood vessel extraction The figure for demonstrating the principle of the correction process by the 2nd Embodiment of this invention.
  • FIG. 1 is a schematic block diagram showing the configuration of an image processing apparatus according to an embodiment of the present invention.
  • the configuration of the image processing apparatus 1 as shown in FIG. 1 is realized by executing a program read into an auxiliary storage device (not shown) on a computer (for example, a personal computer). Further, this program is stored in an information storage medium such as a CD-ROM or distributed via a network such as the Internet and installed in a computer.
  • the image processing apparatus 1 generates a three-dimensional image M0 using a plurality of two-dimensional images captured by the X-ray CT apparatus 2, and the linear structure and the planar structure included in the three-dimensional image M0.
  • the image acquisition unit 10, the detection region setting unit 20, the determination unit 30, the division unit 40, the display unit 50, and the input unit 60 are provided.
  • the image acquisition unit 10 acquires, for example, a plurality of CT images (two-dimensional images) captured by the X-ray CT apparatus 2, and generates a three-dimensional image M0 from the plurality of two-dimensional images.
  • the image acquisition unit 10 may acquire not only CT images but also two-dimensional images such as so-called MRI images, RI images, PET images, X-ray images, and the like.
  • the discrimination unit 30 includes a filtering unit 32, a correction unit 33, and an evaluation unit 34.
  • the filtering unit 32 includes a function (solid sphere model function) representing a solid sphere and a Gaussian kernel to be described later on each of the three-dimensional multi-resolution images Msi.
  • the same filtering as that used in Non-Patent Document 2 is used. That is, the filter kernel of the same size is convolved with each of the three-dimensional multi-resolution images Msi having different resolutions.
  • the solid sphere model function is defined by the radius R of the solid sphere and ⁇ of the Gaussian kernel. These are set to appropriate values based on knowledge such as preliminary analysis.
  • the radius R of the solid sphere is set to a value that is at least larger than ⁇ of the Gaussian kernel.
  • R 2.0 (voxel)
  • ⁇ 0.5 (voxel)
  • a function representing a solid sphere of a plurality of sizes of radius R is used to calculate an evaluation matrix subjected to filtering using a second-order partial differential matrix of the function representing each solid sphere.
  • Hessian analysis using second-order partial differentiation of the solid sphere model function according to Non-Patent Document 2 will be described.
  • the Hessian matrix used for the Hessian analysis is a 3 ⁇ 3 matrix for a three-dimensional image as shown in the following equation (1).
  • the elements Ixx, Ixy, Iyy, Iyz, Izz, and Izz of the Hessian matrix filter the image data of the target image by filtering the second-order partial derivative of the hollow sphere model function f (r) and the Gaussian kernel function g (r). Calculation).
  • this filtering process is performed in Fourier space.
  • the image data of the target image is Fourier transformed.
  • the Fourier-transformed solid sphere model function, the Gaussian kernel function, and the second-order partial derivative are added to the Fourier-transformed target image (FT (image)).
  • FT Fourier-transformed target image
  • a filtering result is acquired by carrying out inverse Fourier transform of the integration result.
  • the obtained filtering result is the same as the filtering result obtained by performing the convolution operation in the real space.
  • the Fourier transform of the solid sphere model function is F ( ⁇ )
  • the Fourier transform of the Gaussian kernel function is G ( ⁇ )
  • x, y, and z are three-dimensional coordinates in the real space, and R is the radius of the solid sphere.
  • R is the radius of the solid sphere.
  • ⁇ ( ⁇ x 2 + ⁇ y 2 + ⁇ z 2 ) 1/2 as shown in the equation (4).
  • Equation (2) For each function used in Equation (2), the Fourier transformed Gaussian kernel function G ( ⁇ ) is shown in the above Equation (3), and the Fourier transformed solid sphere model function F ( ⁇ ) is given in Equation (4). Show.
  • the solid sphere model function F ( ⁇ ) shown in Expression (4) is obtained by performing Fourier transform on the function f (r) representing the solid sphere shown in Expression (5).
  • the Gaussian kernel function G ( ⁇ ) shown in the expression (3) is used to define the differential range of the pixel value of the target image at the time of filtering as in Non-Patent Documents 1 and 2.
  • the linear structure has two eigenvalues out of three, and one has a feature close to 0, as shown in FIG. It has been known.
  • the eigenvalue of the formula (1) has the relationship shown in the formula (6) with respect to the target tissue having a linear structure.
  • the eigenvalue is assumed to be
  • the eigenvalue of the formula (1) has a relationship like the formula (7) with respect to the target tissue having a planar structure.
  • the point-like structure has a characteristic that all three eigenvalues are large.
  • the eigenvalue of the formula (1) has a relationship as shown in the formula (8) with respect to the target tissue having a point-like structure.
  • the correction unit 33 in the present embodiment first eigenvalue-decomposes the Hessian matrix calculated by the filtering unit 32 to obtain three eigenvalues ⁇ 1 , ⁇ 2 , ⁇ 3 (where
  • FIGS. 5A, 5B, 5C, and 5D the principle of suppression of erroneous discrimination according to the present embodiment will be described using FIGS. 5A, 5B, 5C, and 5D, and then a specific correction process of the correction unit 33 based on this principle will be described.
  • FIG. 5A shows the response at each x position of the solid sphere model function (one-dot chain line) and the second-order partial differential (solid line) in the x direction of the solid sphere model function
  • FIG. 5B shows the hollow sphere model.
  • the response of each x position of the function (one-dot chain line) and the first-order partial derivative (dashed line) in the x direction of the hollow sphere model function is shown
  • (C) shows the first-order partial derivative of the solid sphere model function and the hollow sphere model function.
  • the response of each x position of the second-order partial differential (solid line) in the x direction of the solid sphere model function corrected by using the x-direction is shown.
  • the second-order partial differential in each direction of the solid sphere model function shows response waveforms at two spaced positions corresponding to the surface of the solid sphere. Then, when the line segment crossing the structure such as the diameter of the blood vessel and the diameter of the solid sphere represented by the solid sphere model function substantially coincide with each other, two opposing positions of the outline of the structure in the image with respect to each differential direction A large response (expected response) is obtained when both of the contour parts (both ends of the line segment crossing the structure) coincide with the two positions indicating the response waveform of the second-order partial differentiation of the solid sphere model function. Therefore, according to the Hessian analysis by the method of Non-Patent Document 2, structures such as blood vessels can be identified.
  • the second-order partial differential of each direction of the solid sphere model function has a characteristic of showing response waveforms at two positions symmetrical from the center of the solid sphere.
  • the response of the first-order partial differential in the x direction of the hollow sphere model function f 1 (r) having substantially the same size as the solid sphere is the solid sphere model function. It has a feature that it has a response waveform having substantially the same shape and the same sign as one response waveform of the solid sphere model function at the same position as one of the response waveforms of the second-order partial differential in the x direction.
  • the first-order partial differential value of the hollow sphere model function f 1 (r) x direction is obtained by using the characteristic of having a response waveform in the x direction of the solid sphere model function. Used to cancel the response waveform at one position of the second partial differential of direction.
  • the negative response waveforms in the x direction cancel each other because they have opposite signs at the same position, and the positive response waveforms in the x direction have the same position.
  • a response waveform is provided only at one position on the positive side in the x direction as shown in FIG. 5A (C).
  • FIG. 5B shows the solid sphere model function f (r), and FIG. 5B shows the response of the second-order partial differential in the x direction of the solid sphere model function f (r).
  • FIG. 5C shows the hollow sphere model function f 1 (r), and (B) shows the response of the first-order partial differential of the hollow sphere model function f 1 (r) in the x direction.
  • FIG. 5D shows the response of the second-order partial differential in the x direction of the solid sphere model function shown in FIG. 5B and the response of the first-order partial differential of the hollow sphere model function f 1 (r) shown in FIG. 5C. The result of addition is shown.
  • each of the above-mentioned responses is shown by a plurality of xy plan views in which z coordinate values are varied at equal intervals. These plurality of xy plan views are shown in the vertical direction from the one having a large z coordinate. They are shown in order from the top.
  • the correction unit 33 uses the first-order partial differential vector of the hollow sphere model function to determine the position of one position in the second-order partial differentiation in each direction of the solid sphere model function f (r).
  • the eigenvalue of the Hessian is corrected so as to cancel the response waveform.
  • the correction unit 33 first uses the second-order partial differential matrix of the solid sphere model function and the first-order partial differential vector of the hollow sphere model function in the above-described image for correction using Expression (9).
  • the first-order partial differential vector is calculated.
  • the Fourier transform-transformed Gaussian kernel function G ( ⁇ ) and the following Expression (11) are applied to the pixels to be processed of the Fourier-transformed three-dimensional multi-resolution image Msi.
  • First-order partial differential vector to be used for Hessian matrix eigenvalue correction is calculated by performing convolution of the first-order partial differential filter of the hollow sphere model function F 1 ( ⁇ ) subjected to Fourier transform and performing inverse Fourier transform on the filtered result.
  • the function F 1 ( ⁇ ) used in the equation (11) is obtained by substituting the hollow sphere model function f 1 (r) defined by the delta function ⁇ (r ⁇ R 4 ) represented by the equation (10) into the Fourier It is obtained by converting.
  • Equation (9) filtering is performed in Fourier space, but filtering may be performed in real space.
  • x, y, and z are coordinates in a three-dimensional space
  • r is a polar coordinate representation thereof
  • R 4 is a radius of the hollow sphere.
  • the radius R 4 of the hollow sphere represented by the hollow sphere model function is the same as the radius R of the solid sphere represented by the solid sphere model function.
  • the radius R 4 of the hollow sphere is within a range having an effect of canceling the response waveform at one position of the second partial differential of the solid sphere model function, and the radius R 4 of the hollow sphere is larger than the radius R of the solid sphere.
  • the response waveform of the second-order partial differential of the solid sphere model function and the response waveform of the first-order partial differential of the hollow sphere model function are more compatible. It is preferable to exactly match the radius R of the real sphere.
  • the first-order partial differential vectors ( ⁇ 1 , ⁇ 2 , ⁇ 3 ) in the X, Y, and Z directions calculated as described above are the eigenvectors e 1 , e 2 , e of the eigenvalues ⁇ 1 , ⁇ 2 , ⁇ 3. 3
  • the correction unit 33 converts the first-order partial differential vectors ⁇ 1 ′ , ⁇ 2 ′ , ⁇ 3 ′ corresponding to the directions of the eigenvectors e 1 , e 2 , e 3 to the following formula (12 ).
  • Equation (13) the eigenvalue of the Hessian matrix is corrected as shown in Equation (13).
  • the eigenvalues of the evaluation matrix are ⁇ 1 , ⁇ 2 , and ⁇ 3 and ⁇ is a predetermined coefficient.
  • ⁇ is a weight designed in advance so as to most cancel one of the response waveforms appearing at two equidistant positions from the origin of the second-order partial differentiation of the solid sphere model function (the center of the solid sphere). This weight is designed according to the radius R of the solid sphere and the kernel size s.
  • ⁇ 1 + ⁇ ′ is a correction that cancels one of the two response waveforms symmetrical to the center of the solid sphere model function in each differential direction (for example, the response waveform on the left side of FIG. 5A (A)). Indicates the response when it is made.
  • ⁇ 1 - ⁇ ′ represents a response when correction for canceling the other (for example, the response waveform on the right side of FIG. 5A (A)) is performed.
  • the eigenvalue ⁇ 1 ′ is corrected so as to be a smaller value of the two responses
  • a conventional method detects a structure by mistake, that is, if a response is obtained only in one of the two response waveforms symmetric about the center of the solid sphere model function, either of the solid sphere model functions If one of the response waveforms is canceled, at least one of the two responses
  • the structure is correctly detected, that is, both the two opposite contour portions (both ends of the line segment crossing the structure) of the contour of the structure in the image with respect to each differential direction are solid.
  • the evaluation unit 34 uses the corrected eigenvalues ⁇ 1 ′, ⁇ 2 ′, ⁇ 3 ′ instead of the eigenvalues ⁇ 1 , ⁇ 2 , ⁇ 3 in the following formulas (14) and (15), Using the calculated values R A , R B , and R C , an evaluation value L0 (Lineness) of the linear structure likelihood and an evaluation value P0 (Planeness) of the planar structure likelihood in each pixel of the three-dimensional multi-resolution image Msi are obtained. calculate.
  • the evaluation unit 34 determines the point structure or the linear structure based on the corrected eigenvalues ⁇ 1 ′, ⁇ 2 ′, ⁇ 3 ′ of the corrected Hessian matrix and the eigenvectors e 1 , e 2 , e 3. It is possible to discriminate the planar structure, but it is not always necessary to evaluate all of the dotted structure, the linear structure, and the planar structure, and depending on the required specifications, the dotted structure, the linear structure Only a part of the planar structure may be determined. In the present embodiment, since the purpose is to extract a linear structure and a planar structure from a medical image, evaluation values are calculated only for the linear structure and the planar structure.
  • R A , R B and R C are calculated by the following formulas (16) to (19).
  • S 2nd is the power of the second-order partial differential value, and is calculated by the following equation (19).
  • FIG. 6A is a diagram for explaining a response when the size of a structure (the length of a line segment crossing the structure) and the size of a solid sphere model function (diameter of the solid sphere) are substantially the same.
  • FIG. 6A (A) shows an xy plan view of the luminance value of the linear structure.
  • FIG. 6A (B) is a xy plan view showing a response obtained by performing filtering processing using the second-order partial differentiation in the x direction of the solid sphere model function on the image data of FIG. 6A (A).
  • FIG. 6A is a diagram for explaining a response when the size of a structure (the length of a line segment crossing the structure) and the size of a solid sphere model function (diameter of the solid sphere) are substantially the same.
  • FIG. 6A (A) shows an xy plan view of the luminance value of the linear structure.
  • FIG. 6A (B) is a xy plan view showing a response obtained by performing
  • 6A (C) shows a response obtained by performing a filtering process using the first-order partial differentiation of the hollow sphere model function in the x direction on the image data of FIG. 6A (A) in an xy plan view.
  • 6A (D) shows the response obtained by performing the filtering process using the second-order partial differential matrix in the x direction of the solid sphere model function on the image in FIG. 6A (A). What was corrected using the response which performed the filtering process of partial differentiation is represented with an xy top view.
  • the peak of the linear structure appears near the center of the linear structure, and a good response for discriminating the linear structure by the image processing according to the present embodiment is obtained. You can see that it is obtained.
  • FIG. 6B is a diagram for explaining a response when the size of the structure in the x direction (the length of a line segment crossing the structure) is larger than the diameter of the solid sphere model function.
  • FIG. 6B (A) shows an xy plan view of the luminance value of the linear structure.
  • FIG. 6B (B) is an xy plan view showing a response obtained by performing filtering processing using the second-order partial differentiation in the x direction of the solid sphere model function on the image data in FIG. 6B (A).
  • FIG. 6B (C) shows the response obtained by performing the filtering process using the first-order partial differential in the x direction of the hollow sphere model function in an xy plan view.
  • FIG. 6B (A) shows an xy plan view of the luminance value of the linear structure.
  • FIG. 6B (B) is an xy plan view showing a response obtained by performing filtering processing using the second-order partial differentiation in the x direction of the solid sphere model function on the image data in
  • FIG. 6B (D) shows a response obtained by performing filtering processing using the second-order partial differential matrix in the x direction of the solid sphere model function on the image data in FIG. 6B (A). What was corrected using the response which performed the filtering process of the partial partial differentiation is represented with an xy top view.
  • FIG. 6B (B) showing the evaluation values corresponding to the method of Non-Patent Document 2
  • certain responses that cause misclassification due to second-order partial differentiation of the solid sphere model function appear in two places (two negative values).
  • FIG. 6B (D) showing the response by the method of the present embodiment, the response becomes very weak, and it can be seen that the response causing the misclassification is suppressed.
  • FIG. 6A (D) shows the response when in place of the lambda 1 of the formula (13) was applied an element (Ixx) of Hessian matrix shown in Equation (1) . Since the process of calculating the eigenvalue is merely linearly transforming the space in the Hessian matrix and does not affect the response value, the eigenvalue ⁇ obtained by eigenvalue analysis of the Hessian matrix in the above equation (13). This is because even if the element of the Hessian matrix before the eigenvalue analysis is corrected instead of 1, the effect of canceling the response that causes the false detection can be observed in the corrected element.
  • the linear structure-like evaluation value L0 and the planar structure-like evaluation value P0 are calculated, respectively.
  • the calculated evaluation values L0 and P0 are the evaluation values of the corresponding pixel positions of the original three-dimensional image M0, but the evaluation values are calculated at the corresponding pixel positions of all the multi-resolution three-dimensional images Msi. Become.
  • the dividing unit 40 divides the blood vessel region and the region other than the blood vessel of the three-dimensional image M0 based on the linear structure-likeness evaluation value L0 and the planar structure-likeness evaluation value P0 calculated by the determination unit 30. Specifically, a blood vessel region is set as a target region, a region other than the blood vessel region is set as a background region, a discrimination region having a predetermined pixel size is set at all pixel positions in the three-dimensional image M0, and a Graph-Cut region dividing method is used. Thus, the discrimination area is divided into a target area and a background area. Note that the Graph ⁇ Cut region partitioning method is “Yuri Y.
  • the dividing unit 40 uses the same method as described in Japanese Patent Application Laid-Open No. 2011-206531 which is a past application filed by the present inventor, and the evaluation value L0 of the linear structure likelihood calculated by the determining unit 30 and Based on the evaluation value P0 of the planar structure likelihood, the blood vessel region and the region other than the blood vessel of the three-dimensional image M0 are divided into regions.
  • the display unit 50 is a monitor that displays a two-dimensional image or a three-dimensional image, a CRT screen, or the like.
  • the linear structure divided as the target region is displayed on the display unit 50 by volume rendering, so that the entire linear structure or planar structure can be overviewed and the continuity thereof can be visualized.
  • the input unit 60 is a keyboard, a mouse, or the like.
  • FIG. 7 is a flowchart showing processing performed in the present embodiment.
  • the image acquisition unit 10 generates a three-dimensional image M0 from the two-dimensional image captured by the X-ray CT apparatus 2 (step ST1).
  • the detection area setting unit 20 generates a plurality of three-dimensional multi-resolution images Msi having different resolutions by making the three-dimensional image M0 isotropic and performing multi-resolution conversion (step ST2).
  • the filtering unit 32 of the determination unit 30 performs filtering using the second-order partial differentiation of the solid sphere model function and the Gaussian kernel function for each three-dimensional multiresolution image Msi, and calculates the Hessian matrix at each pixel position.
  • the correcting unit 33 corrects the eigenvalue of the Hessian matrix using the first-order partial differential vector of the solid sphere model function (step ST4).
  • the evaluation unit 34 calculates the evaluation value L0 of the linear structure likelihood and the evaluation value P0 of the planar structure likelihood based on the corrected eigenvalue and eigenvector (step ST5).
  • the dividing unit 40 divides the three-dimensional image M0 into a target region (blood vessel region) and a background region using the Graph Cut region dividing method described above (step ST6). Then, the display area 50 displays the divided target area and background area by volume rendering (step ST7), and the process is terminated.
  • FIG. 8A and 8B show an example in which the image processing of the present embodiment is applied to blood vessel discrimination in an actual patient's chest CT image and an example in which the conventional method is applied to blood vessel discrimination in the same CT image.
  • FIG. 8A and 8B the left is an application example (comparative example) of image processing according to the conventional method, and the right is an application example of image processing of the present embodiment.
  • FIG. 8A shows a pseudo three-dimensional image in which the determined blood vessel region is represented by the volume rendering method
  • FIG. 8B shows an enlarged axial tomographic image of a part of FIG. 8A.
  • FIG. 8A in the conventional method, a large number of blood vessels due to misclassification are displayed. However, in the application example of this embodiment, misclassification is greatly suppressed and blood vessels are accurately extracted. I understand. Moreover, in the comparative example (left figure) of FIG. 8B, the portion mistakenly identified as a linear structure inside the blood vessel is shown in light gray, but in the application example (right figure) of this embodiment. It can be seen that there is no portion (light gray portion) determined to be a linear structure inside the blood vessel.
  • the response having the same shape as the one response waveform at the same position as the one response waveform of the second-order partial differentiation of the solid sphere model function and the sign of the positive / negative is reversed.
  • a first-order partial differential vector of a function having a waveform is calculated, and a value obtained by projecting the first-order partial differential vector in the direction of the eigenvector of the Hessian matrix is used to calculate the second-order partial differential of each direction of the solid sphere model function.
  • the response waveform appearing at one position is corrected so as to cancel, so that the contour of the structure and the direction of the solid sphere model function It is possible to suppress erroneous determination that occurs when only one position where the response waveform of the second-order partial differential appears matches, and to improve the accuracy of the evaluation value. For this reason, it is possible to more accurately determine the structure included in the image based on the evaluation value.
  • image filtering processing in Fourier space it is not easy to perform correction to change only some of the undesired values of the filter coefficients.
  • the second-order partial differentiation of the solid sphere model function By paying attention to the response characteristics and using the response characteristics of the first-order partial differentiation of the hollow sphere model function, the problem of misclassification by the technique of Non-Patent Document 1 can be solved satisfactorily.
  • the “first-order” partial differential vector of the hollow sphere model function is used for correction, once the partial differential values in the x, y, and z directions are calculated, the partial differential in any direction is performed. The value can be calculated. That is, by using the “first-order” partial differential vector of the hollow sphere model function in the direction of the eigenvector of the second-order partial differential, it is possible to correct the second-order partial differential in any direction. As a result, the increase in calculation processing for correction is relatively small, and computers having various processing capabilities can be used as the image processing apparatus of this embodiment.
  • the response waveform of one position of the first-order partial differential in each direction of the hollow sphere model function is opposite in sign to the response waveform of one position of the second-order partial differentiation of each direction of the solid sphere model function. Therefore, the response waveform at one position of the second-order partial differential of the solid sphere model function, which causes misclassification, can be preferably canceled out.
  • a Hessian matrix is calculated for each function representing a solid sphere having a plurality of sizes of radius R. Based on the plurality of Hessian matrices, a response that causes misclassification is obtained. The correction value which suppresses is performed, and the evaluation value of the linear structure-likeness and the planar structure-likeness at each pixel position is calculated. For this reason, the structure of the size corresponding to each solid sphere size can be evaluated suitably. In addition, based on this evaluation value, when the target area is divided as in the first embodiment, a structure having a plurality of sizes can be preferably divided as the target area.
  • a constant kernel size ⁇ and a constant radius R of a solid sphere model function are applied to an image having a plurality of resolutions as in this embodiment. Then, a method for performing Hessian analysis may be applied, and a method for performing Hessian analysis by varying the kernel size ⁇ and the radius R of the solid sphere model function may be applied to an image with a constant resolution.
  • the length of a line segment that crosses the detection target structure such as the length of the long axis or the short axis of the structure, is acquired in advance from data obtained experimentally, and the detection target structure is crossed. It is conceivable that the resolution of the three-dimensional multi-resolution image Msi or the kernel size ⁇ and the radius R of the solid sphere are set so that the length of the line segment and the size of the radius R match.
  • both the linear structure-like evaluation value L0 and the planar structure-like evaluation value P0 are calculated.
  • the linear structure-like evaluation value L0 and the planar structure-like evaluation value are calculated. Only one of the evaluation values P0 may be calculated.
  • a second embodiment of the present invention will be described below.
  • the correction processing by the correction unit 33 is different from that of the first embodiment, but other than that is the same as the first embodiment, and the functions of the respective functional blocks are also common.
  • the following description will focus on the differences from the first embodiment, and the description of the same points as in the first embodiment will be omitted.
  • the second-order partial differentiation x of the solid sphere model function is performed.
  • the response waveform located on the negative side of the direction has one positive peak from the negative side in the x direction and one negative peak adjacent thereto, and the response waveform located on the positive side in the x direction is the x direction Focusing on the feature of having one negative peak and one positive peak adjacent to the negative peak from the negative side, the one response waveform is located at the same position as this one response waveform.
  • a response waveform having an approximately identical shape with opposite signs is obtained from the first-order deviation of each direction of the solid sphere model function f 2 (r) (second solid sphere model function) having a slightly smaller size than the solid sphere.
  • a solid sphere model that is slightly larger than this solid sphere. (In the third real sphere model function) function f 3 (r) were found to be possible to create a combination of 1-order partial differential in each direction.
  • the solid sphere model function used by the filtering unit 32 for the filtering process is referred to as a first solid sphere model function.
  • FIG. 9A (A) shows the response of the first partial differential in the x direction of the second solid sphere model function f 2 (r) having the second radius R 2
  • FIG. 9A (B) shows the third FIG. 9A (C) shows the response of the first partial partial differentiation in the x direction of the third solid sphere model function f 3 (r) having the radius R 3 of FIG.
  • the response is obtained by adding the first-order partial differential in the x direction of the third solid-sphere model function f 3 (r) to the first-order partial differential in the x direction of the real sphere model function f 2 (r).
  • the response of the first partial differential of the second and third solid sphere model functions f 2 (r) and f 3 (r) is the center of the solid sphere.
  • FIG. 9A (C) has response waveforms at two symmetrical positions from the center of the solid sphere, and the response waveform located on the negative side in the x direction is one negative peak from the negative side in the x direction. It has one positive peak adjacent to this, and the response waveform located on the positive side in the x direction has one negative peak and one positive peak adjacent to it from the negative side in the x direction. It has characteristics.
  • the first-order partial differentiation in each direction of the second and third solid-sphere model functions f 2 (r) and f 3 (r) is used for filtering by the second-order partial differentiation matrix in each direction of the solid sphere model function.
  • the response waveforms on the positive side in each direction cancel each other because they have opposite signs at the same position, and the response waveforms on the negative side in the x direction are the same at the same position.
  • the second radius R 2 is the first solid sphere model function. Is equal to the distance from the center of the solid sphere to the positive peak of the response waveform located in the positive direction, and the third radius R 3 is the second floor of the first solid sphere model function. This is the same as the distance from the center of the solid sphere in the partial differentiation to the negative peak of the response waveform located in the positive direction.
  • the second radius R 2 (third radius R 3 ) is a range having an effect of canceling the response waveform at one position of the second-order partial differential in each direction of the solid sphere model function.
  • the “distance from the center to the positive (negative) peak of the response waveform” may be larger or smaller than R 2 (third radius R 3 ), but the response of the second partial differential of the solid model function In order to make the waveform more corresponding to the response waveform obtained by combining the first and second partial differentials of the second and third solid sphere model functions f 2 (r) and f 3 (r) in each direction
  • the second radius R 2 (third radius R 3 ) is exactly the same as the distance from the center to the positive (negative) peak of the response waveform.
  • FIG. 9B (A) shows the first-order partial differential response in the x direction of the second solid sphere model function f 2 (r)
  • FIG. 9B (B) shows the third solid sphere model function f 3
  • FIG. 9B (C) shows the response of the first-order partial differentiation in the x direction of the solid sphere model function f 2 (r) with the positive and negative signs reversed.
  • the response when the first-order partial derivative in the x direction of the sphere model function f 3 (r) is added is shown.
  • each of the above-mentioned responses is shown by a plurality of xy plan views with different z coordinate values at equal intervals, and the plurality of xy plan views are arranged in the vertical direction from the one having the largest z coordinate. ing.
  • the brighter (white) the greater the response in the positive direction and the darker (black) the greater the response in the negative direction.
  • the negative response waveform shown in FIG. 9B (C) is a negative peak located on the negative side from the center of the first partial partial differentiation in the x direction of the second solid sphere model function shown in FIG. 9B (A).
  • the negative peak located on the negative side is shifted from the center of the first-order partial differential in the x direction of the third solid sphere model function shown in FIG. 9B (B) and overlaps, and the negative peak and the positive peak from the negative direction overlap. It can be seen that are adjacent waveforms.
  • the positive response waveform shown in FIG. 9B (C) is a positive peak located on the positive side from the center of the first partial differential of the second solid sphere model function shown in FIG.
  • the first-order partial differentials of the second and third solid sphere model functions are combined so as to show the response of FIG. 9B (C), and the first solid sphere model function shown in FIG.
  • the response waveform on the negative side from the center of the solid sphere can be canceled as shown in FIG. 5D.
  • the response waveform on the negative side from the center of the solid sphere can be canceled for each direction.
  • the correction unit 33 uses the first partial differential of the second and third solid sphere model functions to calculate the first solid sphere model function f (r The Hessian is corrected so as to cancel only the response waveform at one position.
  • the Hessian is corrected so as to cancel only the response waveform at one position.
  • the correction unit 33 uses the second radius R 2 instead of the radius R in the equation (2) for the second radius R 2 solid sphere model function to obtain the second solid sphere model function.
  • a first-order partial differential vector is calculated.
  • Each of the values calculated for the second solid sphere model function f 2 (r) by substituting the coefficients l, m, and n corresponding to the respective differential directions into the equation (9) so that ⁇ n).
  • the partial partial differential vectors (s x , s y , s z ) can be calculated respectively.
  • the correction unit 33 uses the third radius R 3 instead of the radius R in the equation (2) for the third solid sphere model function f 3 (r) to obtain the first-order partial differential.
  • a vector (t x , t y , t z ) is calculated.
  • the first-order partial differential vectors ( ⁇ 1 , ⁇ 2 , ⁇ 3 ) in the X, Y, and Z directions calculated as described above are the eigenvectors e 1 , e 2 , e of the eigenvalues ⁇ 1 , ⁇ 2 , ⁇ 3. 3
  • the correction unit 33 converts the first-order partial differential vectors ⁇ 1 ′ , ⁇ 2 ′ , ⁇ 3 ′ corresponding to the directions of the eigenvectors e 1 , e 2 , e 3 to the following formula (12 ).
  • the eigenvalue of the Hessian matrix is corrected based on Expression (13).
  • eigenvalues of the evaluation matrix are ⁇ 1 , ⁇ 2 , and ⁇ 3, and ⁇ and ⁇ are predetermined coefficients.
  • ⁇ and ⁇ are weights designed in advance so as to most cancel one of the response waveforms appearing at two equidistant positions from the origin of the second-order partial differential of the solid sphere model function (the center of the solid sphere). . This weight is designed according to the radius R of the solid sphere and the kernel size s.
  • the evaluation unit 34 replaces the eigenvalues ⁇ 1 , ⁇ 2 , ⁇ 3 with the eigenvalues ⁇ 1 ′, ⁇ 2 ′, ⁇ 3 ′ in the equations (14) and (15). , And using the calculated values R A , R B , and R C , the evaluation value L 0 (Lineness) of the linear structure likelihood and the evaluation value P 0 of the planar structure likelihood in each pixel of the three-dimensional multi-resolution image Msi (Planeness) is calculated.
  • the single response waveform has substantially the same shape and a positive / negative sign at the same position as one response waveform of the second-order partial differentiation of the solid sphere model function.
  • the first-order partial differential vector of a function having a response waveform in which is reversed (a function combining the first-order partial derivatives of solid sphere model functions of two sizes) is calculated, and this first-order partial differential vector is calculated as the Hessian matrix Using the value projected in the direction of the eigenvector, the response waveform appearing at two positions symmetrically spaced from the center of the solid sphere of the second-order partial differentiation of each direction of the solid sphere model function. The response waveform that appears is corrected to cancel.
  • the first-order partial differential of the hollow sphere model function of one size is corrected by using the first-order partial differential of the solid sphere model function of two types for correction.
  • the calculation amount of the correction process by the correction unit 33 is slightly increased as compared with the method of the first embodiment to be used, the second floor in either direction of the contour of the structure and the solid sphere model function as in the first embodiment. It is possible to suppress erroneous discrimination that occurs when only one position where a partial differential response waveform appears, and to improve the accuracy of the evaluation value.
  • the second radius R 2 out from the center of the medium from the center of the actual ball to a positive peak length and the solid sphere length up negative peak, longer lengths
  • the third radius R 3 is equal to the shorter one of the length from the center of the solid sphere to the positive peak and the length from the center of the solid sphere to the negative peak. Therefore, the effect of canceling the response waveform corresponds to one of the response waveforms in which the positive peak and the negative peak of the first solid sphere are adjacent to each other. Can be suppressed.
  • a Hessian matrix is calculated for each function representing the first solid sphere having a plurality of sizes of radius R. Based on the plurality of Hessian matrices, It is preferable to calculate an evaluation value of the linear structure-likeness and the planar structure-likeness at each pixel position by performing correction for suppressing the response. In this case, a structure having a size corresponding to each kernel size can be suitably evaluated. For this reason, the image processing of the second embodiment may be performed on a plurality of resolution images using the first solid sphere model function having a certain radius R, and the image of one kind of resolution may be obtained. On the other hand, the first radius R may be different, and the image processing of the second embodiment may be performed using the first solid sphere model function of the first radius R having a plurality of sizes.
  • the image processing apparatus for direction determination does not include the division unit 40, the evaluation unit 34, and the display unit 50 in the first embodiment, but can include a direction determination unit instead.
  • a direction determination unit instead of such an image processing apparatus, an image acquisition process by the image acquisition unit 10, a detection area setting process by the detection area setting unit 20, and a filtering process by the filtering unit 32 are the same as those in the first embodiment.
  • the direction determination unit can be configured as follows.
  • the direction determination unit performs eigenvalue analysis on the evaluation matrix calculated by the filtering unit 32 to calculate eigenvalues and eigenvectors e 1 , e 2 , e 3 , and based on the calculated eigenvector direction, The normal direction of the planar structure is determined. Specifically, the eigenvalues ⁇ 1 , ⁇ 2 , ⁇ 3 and the eigenvectors e 1 , e 2 , e 3 of the Hessian matrix obtained by the filtering unit 32 are acquired, and the eigenvalues are satisfied if the relationship shown in Expression (6) is satisfied.
  • the direction pointed to by the eigenvector e 1 corresponding to ⁇ 1 is discriminated as the axial direction of the linear structure, and the direction pointed to by the eigen vector e 3 corresponding to the eigenvalue ⁇ 3 is planar if the relationship shown in equation (7) is satisfied. Discriminated from the normal direction of the structure. Note that the direction discriminating unit only selects either the axial direction of the linear structure based on the equation (6) or the normal direction of the planar structure based on the equation (7) based on the eigenvalue and eigenvector of the Hessian matrix. May be determined.
  • the axial direction of the linear structure and the normal direction of the planar structure obtained by the image processing apparatus for determining the direction as described above are, for example, CPR (Curved Planer Reformation / Reconstruction) method of a tubular structure such as a blood vessel or a large intestine.
  • CPR Chemical Planer Reformation / Reconstruction
  • the axial direction of the linear structure obtained by the present invention and the normal direction of the planar structure are, for example, ⁇ Armin Kanitsar, Dominik Fleischmann, Rainer Wegenkittl, Petr Felkel, Eduard Groller: CPR-Curved Planar Reformation. IEEE Visualization. It can be used for various known CPR image generation methods such as “2002.”.
  • the image is a medical image and the structure is a blood vessel
  • the accuracy of the structure determination of the medical image in which many misclassifications have occurred due to the inclusion of blood vessels of various diameters is remarkable. Can be improved.
  • the filtering unit performs filtering using a second-order partial differential matrix of a function representing a solid sphere, filtering using a first-order partial differential matrix of a function representing a solid sphere, and hollow spheres. Since the filtering using the first-order partial differential matrix of the function to be expressed is performed in the Fourier space, the processing cost and the processing speed can be reduced as compared with the case where the filtering processing is performed in the real space. However, the filtering unit performs filtering using a second-order partial differentiation matrix of a function representing a solid sphere, filtering using a first-order partial differentiation matrix of a function representing a solid sphere, and first-order partial differentiation of a function representing a hollow sphere. Filtering using a matrix may be performed in real space.
  • the evaluation unit preferably uses at least one of a local point-like structure, a linear structure, and a planar structure of the structure based on the eigenvector and the corrected eigenvalue. Can be determined.
  • these eigenvectors and eigenvalues may be used to discriminate not only point-like structures, linear structures, and planar structures, but also structures having any shape.
  • the discrimination of blood vessels is taken as an example of a linear structure, but the present invention can also be applied to discrimination of other linear structures such as bronchi. It can also be used to discriminate not only bones but also planar structures such as skin and interlobar pleura.
  • each evaluation process or direction determination process is performed on the linear structure and the planar structure included in the three-dimensional image M0.
  • the linear structure included in the two-dimensional image is used.
  • Each evaluation process or direction determination process may be performed on an object and a planar structure.
  • the linear structure and the planar structure included in the three-dimensional image M0 are divided using the Graph-Cut area dividing method, but other area dividing methods such as the Watershed algorithm are used.
  • the Watershed algorithm is a method of dividing an image so that a boundary is formed between water accumulated in different depressions when water is filled in the terrain where the pixel value information of the image is regarded as altitude. For this reason, it is possible to divide the linear structure and the planar structure by executing the Watershed algorithm after appropriately smoothing the evaluation values L0 and P0 on the three-dimensional image M0. is there.

Landscapes

  • Engineering & Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Health & Medical Sciences (AREA)
  • General Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Radiology & Medical Imaging (AREA)
  • Quality & Reliability (AREA)
  • Software Systems (AREA)
  • Image Processing (AREA)
  • Image Analysis (AREA)
  • Apparatus For Radiation Diagnosis (AREA)

Abstract

【課題】画像に含まれる構造物の誤判別を抑制する。 【解決手段】画像に中実球関数の２階偏微分を用いたフィルタリングを施してヘッセ行列を算出するフィルタリング部（３２）と、ヘッセ行列の固有値および固有ベクトルを用いて画像中に含まれる構造物を判別する評価部（３０）とを備え、フィルタリング部（３２）が、画像に中実球の半径と同じ半径を有する中空球関数の１階偏微分を用いたフィルタリングを施して１階偏微分ベクトルを算出し、１階偏微分ベクトルを固有ベクトルの方向に射影した値を用いて、中実球を表す関数の各方向の２階偏微分の、中実球の中心から対称に離間した２つの位置に表れる応答波形のうち一方の応答波形を打ち消す補正を行う補正部（３３）を備える。

Description

画像処理装置および方法並びにプログラム
　本発明は、画像中の構造物の構造を判別するための画像処理装置および方法並びに画像処理方法をコンピュータに実行させるためのプログラムに関するものである。
　近年、医療機器（例えば多検出器型ＣＴ等）の進歩により質の高い３次元画像が画像診断に用いられる。３次元画像は多数の２次元画像から構成され情報量が多いため、医師が所望の観察部位を見つけ診断することに時間を要する場合がある。そこで、注目する臓器を抽出しＭＩＰ、ＶＲ、ＣＰＲ等の表示を行うことにより、臓器全体や病変の視認性を高め診断の効率化を図ることが行われている。
　一方、医用画像中の血管および骨を抽出する手法として、ヘッセ行列を用いたヘッセ解析が提案されている（非特許文献１参照）。ヘッセ解析は、ガウシアンカーネル等の所定のフィルタについての二次微分カーネルを用いることにより算出した、２階の偏微分係数を要素とするヘッセ行列の固有値を解析することにより、画像中の局所構造が点、線および面のいずれであるかを判別するものである。ヘッセ解析を用いることにより血管は線状構造物として、骨は面状構造物として判別することが可能である。
　しかしながら、非特許文献１の手法は、線状構造物の周辺に他の構造物（周辺構造物）が位置している場合に、この周辺構造物の一部を誤って線状構造物として判別する場合があった。非特許文献２の方法は、非特許文献１に提案されたフィルタを、球（内側を１とし、外側を０とする中実球）を表す関数（中実球モデル関数）の２階偏微分をたたみ込んだものに改良することにより、各画素のフィルタリング範囲を画素を中心とする球表面に限定し、周辺構造物のフィルタリング結果に及ぼす影響を小さくすることができる。
A. F. Frangi et. al., "Multiscale vessel enhancement filtering", Proceedings of MICCAI, pp130-137, 1998. M. Law et. al., "Three Dimensional Curvilinear Structure Detection Using Optimally Oriented Flux",Proceedings of ECCV, pp368-382, 2008.
　しかしながら、非特許文献２の手法を用いて、様々な太さの血管を含む医用画像から線状構造である血管を判別すると、血管の内側に実際の血管よりも細い血管が誤って認識されてしまう場合がある。言い換えると、非特許文献２の手法で線状構造物を判別する場合に、構造物の輪郭部分のうち、中実球モデル関数の中実球の曲率半径より大きい曲率半径を有する輪郭部分で、中実球モデル関数の表す中実球の直径と略同一の直径を有する線状の構造物が存在すると誤って判別されてしまうという問題がある。
　上記問題は、中実球モデル関数の中実球よりも曲率半径が大きい構造物の輪郭部分であればどこでも起こりうるため、本発明は、この問題に鑑みて、ヘッセ解析において中実球モデル関数を用いたフィルタリングを行う画像処理方法において、中実球モデル関数の表す中実球よりも曲率半径が大きい構造物の輪郭部分で生じる構造物の誤判別を防止することを目的とする。
　本第１の発明による画像処理装置は、画像中の各画素位置において、中実球を表す関数の２階偏微分を用いたフィルタリングを施してヘッセ行列を算出するフィルタリング部と、算出されたヘッセ行列を固有値解析して得られた固有値および固有ベクトルを用いて、画像中に含まれる構造物を判別する評価部とを備え、フィルタリング部が、画像中の各画素位置において、中実球の半径と同じ半径を有する中空球を表す関数の１階偏微分を用いたフィルタリングを施して１階偏微分ベクトルを取得し、取得した１階偏微分ベクトルを固有ベクトルの方向に射影した値を用いて、中実球を表す関数の各方向の２階偏微分の、中実球の中心から対称に離間した２つの位置に表れる応答波形のうち一方の応答波形を打ち消す補正を行う補正部を備えたものであることを特徴とするものである。
　本第１の発明による画像処理方法は、画像中の各画素位置において、中実球を表す関数の２階偏微分を用いたフィルタリングを施してヘッセ行列を算出するフィルタリング工程と、算出されたヘッセ行列を固有値解析して得られた固有値および固有ベクトルを用いて、画像中に含まれる構造物を判別する評価工程とを備えた画像処理方法であって、フィルタリング工程が、画像中の各画素位置において、中実球の半径と同じ半径を有する中空球を表す関数の１階偏微分を用いたフィルタリングを施して１階偏微分ベクトルを取得し、取得した１階偏微分ベクトルを固有ベクトルの方向に射影した値を用いて、中実球を表す関数の各方向の２階偏微分の、中実球の中心から対称に離間した２つの位置に表れる応答波形のうち一方の応答波形を打ち消す補正を行うものである。
　なお、本第１の発明による画像処理方法を、コンピュータに実行させるためのプログラムとして提供してもよい。
　上記、「中実球の半径と同じ半径を有する中空球」とは、中実球と中空球の半径が厳密に一致する場合だけでなく、中空球を表す関数の１階偏微分ベクトルを用いて、中実球を表す関数の各方向の２階偏微分の中実球の中心から対称に離間した２つの位置に表れる応答波形のうち一方の位置の応答波形を打ち消す効果が得られる範囲であれば、中空球の半径が中実球の半径より大きい場合や小さい場合も含む。上記「一方の位置の応答波形を打ち消す効果」を良好に得るために、中空球の半径と中実球の半径ができるだけ等しいことが好ましく、例えば、中空球の半径と中実球の半径の差が２割以下であることが好ましく、中空球の半径と中実球の半径の差が１割以下であることが好ましい。
　本第１の発明による画像処理装置において、評価行列の固有値をλ1、λ２、λ３とし、補正後の固有値をλ1 ’、λ２ ’、λ３ ’ とし、固有ベクトルをｅ１＝（ｘ１,ｙ１、ｚ１），ｅ２＝（ｘ２,ｙ２、ｚ２），ｅ３＝（ｘ３,ｙ３、ｚ３）とし、第１の半径を有する中空球を表す関数の１階偏微分ベクトルを（ρ１,ρ２,ρ３）とすると、補正部は、下記式（１２）で算出されるρ１ ’,ρ２ ’,ρ３ ’および所定の係数αを用いて、下記式（１３）に示すように、固有値を補正することにより、一方の応答波形を打ち消す補正を行うものであることが好ましい。
Figure JPOXMLDOC01-appb-M000001
Figure JPOXMLDOC01-appb-M000002
　本第１の発明による画像処理装置において、中空球を表す関数は、下記式（１０）により表されるものであることが好ましい。
Figure JPOXMLDOC01-appb-M000003
 ここで、ｘ、ｙ、ｚは３次元空間の座標、ｒはその極座標表現であり、Ｒ４は中空球の半径とする。
　本第２の発明による画像処理装置は、画像中の各画素位置において、中実球を表す関数の２階偏微分を用いたフィルタリングを施してヘッセ行列を算出するフィルタリング部と、算出されたヘッセ行列を固有値解析して得られた固有値および固有ベクトルを用いて、画像中に含まれる構造物を判別する評価部とを備え、フィルタリング部が、画像中の各画素位置において、第１の中実球の半径である第１の半径より大きい第２の半径を有する第２の中実球を表す関数の１階偏微分を用いたフィルタリングを施して１階偏微分ベクトルを算出し、かつ、画像中の各画素位置において、第１の半径より小さい第３の半径を有する第３の中実球を表す関数の１階偏微分を用いたフィルタリングを施して１階偏微分ベクトルを算出し、第２の中実球を表す関数の１階偏微分ベクトルを固有ベクトルの方向に射影した値および第３の中実球を表す関数の１階偏微分ベクトルを固有ベクトルの方向に射影した値を用いて、中実球を表す関数の各方向の２階偏微分の、中実球の中心から対称に離間した２つの位置に表れる応答波形のうち一方の位置の応答波形を打ち消す補正を行う補正部を備えたものである。
　本第２の発明による画像処理装置において、一方の位置の応答波形は、１つの正のピークと１つ負のピークが隣接する波形であり、第２の半径は、中実球の中心から正のピークまでの長さと中実球の中心から負のピークまでの長さのうち、長い方の長さと一致するものであり、第３の半径は、中実球の中心から正のピークまでの長さと中実球の中心から負のピークまでの長さのうち、短い方の長さと一致するものであることが好ましい。
　上記、「第２の半径は、第１の中実球の中心から前記正のピークまでの長さと第１の中実球の中心から前記負のピークまでの長さのうち、長い方の長さと一致する」とは、第２の半径が中実球の中心から正のピークまでの長さと中実球の中心から負のピークまでの長さのうち、長い方の長さ（以下、第１の長さ）と厳密に一致する場合だけでなく、第２の中実球を表す関数の１階偏微分ベクトルを重み付けしたベクトルおよび第３の中実球を表す関数の１階偏微分ベクトルを重み付けしたベクトルを用いて、第１の中実球を表す関数の各方向の２階偏微分の、中実球の中心から対称に離間した２つの位置に表れる応答波形のうち一方の位置の応答波形を打ち消す効果が得られる範囲であれば、第２の半径が第１の長さより大きい場合や小さい場合も含む。上記「一方の位置の応答波形を打ち消す効果」を良好に得るために、第２の半径と第１の長さができるだけ等しいことが好ましく、例えば、第２の半径と第１の長さとの差が２割以下であることが好ましく、１割以下であることがさらに好ましい。
　同様に、上記、「第３の半径は、第１の中実球の中心から前記正のピークまでの長さと第１の中実球の中心から前記負のピークまでの長さのうち、長い方の長さと一致する」とは、第３の半径が中実球の中心から正のピークまでの長さと中実球の中心から負のピークまでの長さのうち、短い方の長さ（以下、第２の長さ）と厳密に一致する場合だけでなく、第３の中実球を表す関数の１階偏微分ベクトルを重み付けしたベクトルおよび第３の中実球を表す関数の１階偏微分ベクトルを重み付けしたベクトルを用いて、第１の中実球を表す関数の各方向の２階偏微分の、中実球の中心から対称に離間した２つの位置に表れる応答波形のうち一方の位置の応答波形を打ち消す効果が得られる範囲であれば、第３の半径が第２の長さより大きい場合や小さい場合も含む。上記「一方の位置の応答波形を打ち消す効果」を良好に得るために、第３の半径と第２の長さができるだけ等しいことが好ましく、例えば、第３の半径と第２の長さとの差が２割以下であることが好ましく、１割以下であることがさらに好ましい。
　本第１の発明による画像処理装置において、フィルタリング部は、複数のサイズの中実球を表す関数について、各中実球を表す関数の２階偏微分行列によるフィルタリングを施した評価行列をそれぞれ算出するものであることが好ましい。
　本第２の発明による画像処理装置において、フィルタリング部は、複数のサイズの第１の中実球を表す関数について、各第１の中実球を表す関数の２階偏微分行列によるフィルタリングを施した評価行列をそれぞれ算出するものであることが好ましい。
　本第１および第２の発明による画像処理装置において、画像が医用画像であり、構造物が血管であることが好ましい。
　本第１および第２の発明による画像処理装置において、フィルタリング部は、中実球を表す関数の２階偏微分行列を用いたフィルタリングを、フーリエ空間において実施するものであることが好ましい。
　本第１および第２の発明による画像処理装置において、評価部は、構造物の局所的な点状構造、線状構造、面状構造の少なくとも１つを判別するものとすることが好ましい。
　本第１発明によれば、画像中の各画素位置において、中実球を表す関数の２階偏微分を用いたフィルタリングを施してヘッセ行列を算出するフィルタリング部と、算出されたヘッセ行列を固有値解析して得られた固有値および固有ベクトルを用いて、画像中に含まれる構造物を判別する評価部とを備え、フィルタリング部が、画像中の各画素位置において、中実球の半径と同じ半径を有する中空球を表す関数の１階偏微分を用いたフィルタリングを施して１階偏微分ベクトルを取得し、取得した１階偏微分ベクトルを固有ベクトルの方向に射影した値を用いて、中実球を表す関数の各方向の２階偏微分の、中実球の中心から対称に離間した２つの位置に表れる応答波形のうち一方の応答波形を打ち消す補正を行う補正部を備えている。このため、中実球を表す関数が表す中実球の曲率より大きい曲率を有する構造物の輪郭部分と、中実球を表す関数の各方向の２階偏微分の応答波形の１つの位置のみが一致した場合に生じる誤判別を抑制でき、評価値の精度を向上することができる。このため、評価値に基づいて画像に含まれる構造物をより正確に判別することができる。
　本第２発明によれば、画像中の各画素位置において、中実球を表す関数の２階偏微分を用いたフィルタリングを施してヘッセ行列を算出するフィルタリング部と、算出されたヘッセ行列を固有値解析して得られた固有値および固有ベクトルを用いて、画像中に含まれる構造物を判別する評価部とを備え、フィルタリング部が、第１の中実球の半径である第１の半径より大きい第２の半径を有する第２の中実球を表す関数の１階偏微分ベクトルを算出し、かつ、第１の半径より小さい第３の半径を有する第３の中実球を表す関数の１階偏微分ベクトルを算出し、第２の中実球を表す関数の１階偏微分ベクトルを固有ベクトルの方向に射影した値および第３の中実球を表す関数の１階偏微分ベクトルを固有ベクトルの方向に射影した値を用いて、中実球を表す関数の各方向の２階偏微分の、中実球の中心から対称に離間した２つの位置に表れる応答波形のうち一方の位置の応答波形を打ち消す補正を行う補正部を備えている。このため、中実球を表す関数が表す中実球の曲率より大きい曲率を有する構造物の輪郭部分と、中実球を表す関数の各方向の２階偏微分の応答波形の１つの位置のみが一致した場合に生じる誤判別を抑制でき、評価値の精度を向上することができる。このため、評価値に基づいて画像に含まれる構造物をより正確に判別することができる。
本発明の第１の実施形態による画像処理装置の構成を示す概略ブロック図 多重解像度変換を説明するための図 線状構造の固有値を説明するための図 面状構造の固有値を説明するための図 本発明の第１の実施形態による補正処理の原理を説明するための図 本発明の第１の実施形態のフィルタリングに用いる中実球モデル関数のｘ方向の２階偏微分のレスポンスを説明するための図 本発明の第１の実施形態の補正処理に用いる中空球モデル関数のｘ方向の１階偏微分のレスポンスを説明するための図 本発明の第１の実施形態の中空球モデル関数の１階偏微分ベクトルにより補正された中実球のｘ方向の２階偏微分のレスポンスを説明するための図 本発明の第１の実施形態による補正処理を実施したレスポンスを説明する図（構造物のサイズと中実球モデル関数の表す中実球のサイズが一致する場合） 本発明の第１の実施形態による補正処理を実施したレスポンスを説明する図（構造物のサイズに対して中実球モデル関数の表す中実球のサイズが大きい場合） 本発明の第１の実施形態において行われる処理を示すフローチャート 本発明の第１の実施形態による画像処理を血管抽出に適用した例（擬似３次元画像） 本発明の第１の実施形態による画像処理を血管抽出に適用した例（断層像） 本発明の第２の実施形態による補正処理の原理を説明するための図 本発明の第２の実施形態の第２および第３の中実球モデル関数の１階偏微分ベクトルによる補正の原理を説明するための図
　以下、図面を参照して本発明の実施形態について説明する。図１は本発明の実施形態による画像処理装置の構成を示す概略ブロック図である。なお、図１のような画像処理装置１の構成は、補助記憶装置（不図示）に読み込まれたプログラムをコンピュータ（例えばパーソナルコンピュータ等）上で実行することにより実現される。また、このプログラムは、ＣＤ－ＲＯＭ等の情報記憶媒体に記憶され、もしくはインターネット等のネットワークを介して配布され、コンピュータにインストールされることになる。
　画像処理装置１は、例えばＸ線ＣＴ装置２により撮像された複数の２次元画像を用いて、３次元画像Ｍ０を生成し、この３次元画像Ｍ０に含まれる線状構造物と面状構造物とを自動的に分割するものであって、画像取得部１０、検出領域設定部２０、判別部３０、分割部４０、表示部５０および入力部６０を備える。
　画像取得部１０は、例えばＸ線ＣＴ装置２により撮像された複数のＣＴ画像（２次元画像）を取得し、複数の２次元画像から３次元画像Ｍ０を生成する。なお、画像取得部１０は、ＣＴ画像のみならず、いわゆるＭＲＩ画像、ＲＩ画像、ＰＥＴ画像、Ｘ線画像等の２次元画像を取得するものであってもよい。
　検出領域設定部２０は、まず３次元画像Ｍ０のボクセルサイズを等方化する。例えば、３次元画像Ｍ０のボクセルサイズが、３次元画像Ｍ０におけるＸ，Ｙ，Ｚ方向のそれぞれにおいて、０．３ｍｍ、０．３ｍｍ、０．６ｍｍであった場合、（Ｘ，Ｙ，Ｚ）＝（０．５，０．５，０．５）（ｍｍ）に等方化する。
　検出領域設定部２０は、等方化後に３次元画像Ｍ０を多重解像度変換して、図２に示すように、解像度が異なる複数の３次元多重解像度画像Ｍｓｉ（ｉ＝０～ｎ）（ガウシアンピラミッド）を生成する。なお、ｉ＝０は３次元画像Ｍ０と同一解像度、ｉ＝ｎは最低解像度を表す。√２倍刻みで画像を縮小し、３次元多重解像度画像Ｍｓｉのボクセルサイズは、解像度が高い順に、（Ｘ，Ｙ，Ｚ）＝（０．５，０．５，０．５）、（０．７，０．７，０．７）、（１．０，１．０，１．０）…となる。
　判別部３０は、フィルタリング部３２および補正部３３および評価部３４を備える。フィルタリング部３２は、ヘッセ行列（評価行列）を用いたヘッセ解析を行うために、３次元多重解像度画像Ｍｓｉのそれぞれに、後述の中実球を表す関数（中実球モデル関数）とガウシアンカーネルを用いた、非特許文献２による手法と同様のフィルタリングを行う。すなわち、解像度が異なる３次元多重解像度画像Ｍｓｉのそれぞれに対して、同一サイズのフィルタカーネルを畳み込む。なお、中実球モデル関数は中実球の半径Ｒとガウスカーネルのσで規定される。これらは、予備的な解析等の知見に基づいて応じた適切な値が設定される。なお、中実球の半径Ｒはガウスカーネルのσよりも少なくとも大きい値が設定される。
　解像度が異なる３次元多重解像度画像Ｍｓｉのそれぞれに対して、同一サイズのフィルタカーネル（例えばR=2.0(voxel)、σ=0.5(voxel)とする）を畳み込むことにより、３次元画像Ｍ０に対して、実質的に異なるサイズのフィルタカーネルを適用することとなるため、異なるサイズの点状構造物、線状構造物（例えば、血管）および面状構造物（例えば皮質骨等の骨）を検出することができる。言い換えると、複数サイズの半径Ｒの中実球を表す関数を用いて、該各中実球を表す関数の２階偏微分行列を用いたフィルタリングを施した評価行列をそれぞれ算出している。
　以下、非特許文献２による中実球モデル関数の２階偏微分を用いたヘッセ解析について説明する。ヘッセ解析に使用されるヘッセ行列は３次元画像に対しては下記の式（１）に示すように３×３の行列となる。
Figure JPOXMLDOC01-appb-M000004
　上記ヘッセ行列の各要素Ｉｘｘ、Ｉｘｙ、Ｉｙｙ、Ｉｙｚ、Ｉｚｚ、Ｉｚｘは、対象画像の画像データを中空球モデル関数ｆ（ｒ）およびガウシアンカーネル関数ｇ（ｒ）の２階偏微分をフィルタリング（コンボリューション演算）して算出する。
　非特許文献２ではこのフィルタリング処理をフーリエ空間で行う。まず対象画像の画像データをフーリエ変換する。次いで、このフーリエ変換された対象画像（FT(image)）に、それぞれフーリエ変換された中実球モデル関数、ガウシアンカーネル関数、および２階偏微分を積算する。そして、積算結果を逆フーリエ変換することによってフィルタリング結果を取得する。なお、この取得されたフィルタリング結果は、実空間でコンボリューション演算をして得られたフィルタリング結果と同じものとなる。中実球モデル関数のフーリエ変換をF（ν）、ガウシアンカーネル関数のフーリエ変換をＧ（ν）としたとき、以上の関係は式（２）のように表される。
Figure JPOXMLDOC01-appb-M000005
Figure JPOXMLDOC01-appb-M000006
　ここでｘ、ｙ、ｚは実空間における３次元座標であり、Ｒは中実球の半径である。また３次元フーリエ空間の変数（周波数）の極座標表現を式（４）に示したように、ν＝（νｘ ２＋νｙ ２＋νｚ ２）１／２で表す。
　式（２）に用いられている各関数について、フーリエ変換したガウシアンカーネル関数Ｇ（ν）を上記式（３）に示し、フーリエ変換した中実球モデル関数Ｆ（ν）を式（４）に示す。式（４）に示す中実球モデル関数Ｆ（ν）は、式（５）に示す中実球を表す関数ｆ（ｒ）をフーリエ変換することにより求められる。また、式（３）に示すガウシアンカーネル関数Ｇ（ν）は、非特許文献１および２と同様にフィルタリングの際に対象画像の画素値の微分範囲を規定するために用いられている。
　また、式（２）において、（２πνｘ）ｌ×（２πνｙ）ｍ×（２πνｚ）ｎの部分はフーリエ空間における微分処理に相当し、ｌ＋ｍ＋ｎ＝２（０＜ｌ、０＜ｍ、０＜ｎ）となるように、各微分方向に対応する係数ｌ、ｍ、ｎをそれぞれ式（４）に代入することにより、各要素Ｉｘｘ、Ｉｘｙ、Ｉｙｙ、Ｉｙｚ、Ｉｚｚ、Ｉｚｘに対応する値をそれぞれ算出できる。例えば、３次元多重解像度画像Ｍｓｉの処理の対象となる画素において、ヘッセ行列におけるＸ方向における２階偏微分である要素Ｉｘｘは、式（２）にＸ方向にｌ＝２を代入し、Ｙ方向およびＺ方向にそれぞれｍ＝ｎ＝０を代入することにより算出できる。また、３次元多重解像度画像Ｍｓｉの処理の対象となる画素において、ヘッセ行列における要素Ｉｘｙは、式（２）にＸ方向およびＹ方向のそれぞれにｌ＝ｍ＝１を代入し、Ｚ方向にｎ＝０を代入することにより算出できる。
　このようにして算出したヘッセ行列を固有値分解して固有値を得たとき、線状構造は図３に示すように、３つのうち２つの固有値の値が大きく、１つが０に近い特徴を持つことが知られている。例えば、式（１）の固有値は、線状構造からなる対象組織に対して、式（６）に示す関係を有する。なお、固有値は｜λ１｜＜｜λ２｜＜｜λ３｜とする。
Figure JPOXMLDOC01-appb-M000007
　また、面状構造は図４に示すように、３つのうち１つの固有値の値が大きく、２つが０に近い特徴を持つことが知られている。例えば、式（１）の固有値は、面状構造からなる対象組織に対して、式（７）のような関係を有する。
Figure JPOXMLDOC01-appb-M000008
　また、点状構造は、３つの固有値の全てが大きい特徴を持つことが知られている。例えば、式（１）の固有値は、点状構造からなる対象組織に対して、式（８）のような関係を有する。
Figure JPOXMLDOC01-appb-M000009
　したがって、固有値から線状構造らしさ、面状構造および点状構造らしさを判別することができ、判別結果を用いて、３次元画像Ｍ０において、線状構造物である血管領域、および面状構造物である骨領域を分割することができる。
　本実施形態における補正部３３は、まず、フィルタリング部３２が算出したヘッセ行列を固有値分解し、３つの固有値λ１，λ２，λ３（ただし、｜λ１｜≦｜λ２｜≦｜λ３｜とする）を算出する。
　ここで、図５Ａ、５Ｂ、５Ｃ、５Ｄを用いて本実施形態による誤判別の抑制の原理を説明し、その後、この原理に基づく具体的な補正部３３の補正処理を説明する。
　図５Ａは、（Ａ）に中実球モデル関数（一点鎖線）と中実球モデル関数のｘ方向の２階偏微分（実線）の各ｘ位置のレスポンスを示し、（Ｂ）に中空球モデル関数（一点鎖線）と中空球モデル関数のｘ方向の１階偏微分（破線）の各ｘ位置のレスポンスを示し、（Ｃ）に中実球モデル関数と中空球モデル関数の１階偏微分を用いて補正された中実球モデル関数のｘ方向の２階偏微分（実線）の各ｘ位置のレスポンスを示す。
　図５Ａ（Ａ）に示すように、中実球モデル関数の各方向の２階偏微分は中実球の表面に相当する離間した２つの位置に応答波形を示す。そして、血管の直径などの構造物を横断する線分と中実球モデル関数の表す中実球の直径が略一致した場合、各微分方向について画像中の構造物の輪郭のうち対向する２箇所の輪郭部分（構造物を横断する線分の両端）の両方が中実球モデル関数の２階偏微分の上記応答波形を示す２つの位置とそれぞれ一致して大きいレスポンス（期待したレスポンス）が得られるため、非特許文献２の手法によるヘッセ解析によれば血管などの構造物が判別できる。
　一方、中実球より曲率の大きい構造物の輪郭部分において、画像中の構造物の輪郭部分に中実球モデル関数の各微分方向について２階偏微分による応答波形の１つの位置のみが一致した場合にも、この一つの位置の応答波形に対応する所定の大きさのレスポンスが得られる。本発明者らは、この所定の大きさのレスポンスと期待したレスポンスとの区別がつかないことが中実球の直径と略同じ径の構造物が誤判別されてしまう原因となっていると推測した。そして、この中実球モデル関数の各方向の２階偏微分の１つの位置における応答波形を打ち消すことにより、誤判別が解消可能であることを見出した。そして、図５Ａ（Ａ）に示すように、中実球モデル関数の各方向２階偏微分は中実球の中心から対称な２つの位置に応答波形を示すという特徴を有するため、中実球モデル関数の２階偏微分の１つの応答波形と同じ位置にこの１つの応答波形と略同形状であるとともに正負の符号を逆にした応答波形を有する関数を利用することにより、中実球モデル関数の２階偏微分行列の１つの位置における応答波形を打ち消すことに注目した。
　本実施形態では、図５Ａ（Ｂ）に示すように、中実球と略同じサイズの中空球モデル関数ｆ１（ｒ）のｘ方向の１階偏微分のレスポンスが、中実球モデル関数のｘ方向の２階偏微分の一方の応答波形と同じ位置に、中実球モデル関数の１つの応答波形と略同形状かつ同符号の応答波形を有するという特徴を有し、さらに、中実球モデル関数のｘ方向の２階偏微分の他方の応答波形と同じ位置に、中実球モデル関数のｘ方向の２階偏微分の他方の応答波形と略同形状であるとともに正負の符号を逆にした応答波形を有するという特徴を有することを利用して、中実球と略同じサイズの中空球モデル関数ｆ１（ｒ）ｘ方向の１階偏微分値を、中実球モデル関数のｘ方向の２階偏微分の１つの位置における応答波形を打ち消すために用いる。
　つまり、中実球モデル関数のｘ方向の２階偏微分を用いてフィルタリングを行ったレスポンスと、中実球モデル関数の表す中実球と略同じ半径の中空球を表す中空球モデル関数の１階偏微分を用いてフィルタリングを行ったレスポンスとを加算すると、ｘ方向の負側の応答波形は、同じ位置において逆の符号を有するため打ち消し合い、ｘ方向の正側の応答波形は、同じ位置において同じ符号を有するため強め合うため、図５Ａ（Ｃ）に示すようにｘ方向において正側の１つの位置にのみ応答波形を有するものとなる。
　図５Ｂ（Ａ）は、中実球モデル関数ｆ（ｒ）を示し、（Ｂ）は中実球モデル関数ｆ（ｒ）のｘ方向の２階偏微分のレスポンスを示す。図５Ｃ（Ａ）は、中空球モデル関数ｆ１（ｒ）を示し、（Ｂ）は中空球モデル関数ｆ１（ｒ）のｘ方向の１階偏微分のレスポンスを示す。図５Ｄは、図５Ｂの（Ｂ）に示す中実球モデル関数のｘ方向の２階偏微分のレスポンスと図５Ｃに示す中空球モデル関数ｆ１（ｒ）の１階偏微分のレスポンスとを加算した結果を示す。なお、図５Ｂ～Ｄにおいて、各上述のレスポンスを等間隔にｚ座標値を異ならせた複数のｘｙ平面図により示しており、これらの複数のｘｙ平面図をｚ座標の大きいものから垂直方向に上から順に並べて示している。
　図５Ｂ～Ｄの各レスポンスを表すｘｙ平面図において、明るい（白い）ほどレスポンスが正の方向に大きくなり、暗い（黒い）ほどレスポンスが負の方向に大きくなる。図５Ｂ（Ｂ）に示す中実球の中心から負側に位置する正のピークと負のピークの隣接した応答波形と図５Ｃ（Ｂ）に示す中実球の中心から負側に位置する負のピークと正のピークの隣接した応答波形が相互に打ち消し合い、図５Ｄでは、中実球の中心から負側には応答波形が見られなくなっているのが分かる。同様に各方向について中実球の中心から負側の応答波形を打ち消すことができる。
　本実施形態では、上記原理に基づいて、補正部３３が中空球モデル関数の１階偏微分ベクトルを用いて中実球モデル関数ｆ（ｒ）の各方向の２階偏微分における１つの位置の応答波形を打ち消すようにヘッセ行列の固有値を補正する。具体的な補正方法を以下に説明する。
　先述の画像に中実球モデル関数の２階偏微分行列と中空球モデル関数の１階偏微分ベクトルを、本実施形態において補正部３３は、まず、式（９）を用いて補正に用いるための１階偏微分ベクトルを算出する。詳細には、式（９）に示すように、フーリエ変換した３次元多重解像度画像Ｍｓｉの処理の対象となる画素に、フーリエ変換したガウシアンカーネル関数Ｇ（ν）と下記の式（１１）に表すフーリエ変換した中空球モデル関数Ｆ１（ν）の１階偏微分フィルタの畳み込みを行い、そのフィルタリング結果を逆フーリエ変換することにより、ヘッセ行列の固有値補正に用いるための１階偏微分ベクトルを算出する。なお、式（１１）に用いた関数Ｆ１（ν）は、式（１０）で表されるデルタ関数δ（ｒ－Ｒ４）により定義される中空球モデル関数ｆ１（ｒ）を、フーリエ変換することにより得られる。
Figure JPOXMLDOC01-appb-M000010
Figure JPOXMLDOC01-appb-M000011
Figure JPOXMLDOC01-appb-M000012
　また、式（９）において、（２πνｘ）ｌ×（２πνｙ）ｍ×（２πνｚ）ｎの部分はフーリエ空間における微分処理に相当し、ｌ＋ｍ＋ｎ＝１（０＜ｌ、０＜ｍ、０＜ｎ）となるように、各微分方向に対応する係数ｌ、ｍ、ｎをそれぞれ式（９）に代入することにより、1階偏微分ベクトルの各要素ρ１、ρ２、ρ３に対応する値をそれぞれ算出できる。
　なお、式（９）においては、フーリエ空間においてフィルタリングを行っているが、実空間においてフィルタリングを行ってもよい。なお、式（１０）において、ｘ、ｙ、ｚは３次元空間の座標、ｒはその極座標表現であり、Ｒ４は中空球の半径とする。また、本実施形態において中空球モデル関数の表す中空球の半径Ｒ４は、中実球モデル関数の表す中実球の半径Ｒと同じである。なお、中空球の半径Ｒ４は、上記中実球モデル関数の２階偏微分の１つの位置における応答波形を打ち消す効果を有する範囲で、中実球の半径Ｒよりも中空球の半径Ｒ４が大きくてもよく、小さくてもよいが、中実球モデル関数の２階偏微分の応答波形と、中空球モデル関数の１階偏微分の応答波形をより対応したものとするために、中実球の半径Ｒと厳密に一致することが好ましい。
　上述したように算出したＸ方向、Ｙ方向およびＺ方向の１階偏微分ベクトル（ρ１,ρ２,ρ３）は、固有値λ１，λ２，λ３の固有ベクトルｅ１，ｅ２，ｅ３の方向とずれているため、補正部３３は、固有ベクトルｅ１，ｅ２，ｅ３の方向に対応した１階偏微分ベクトルρ１ ’,ρ２ ’,ρ３ ’を下記の式（１２）により算出する。なお、評価行列の固有ベクトルをｅ１＝（ｘ１,ｙ１、ｚ１），ｅ２＝（ｘ２,ｙ２、ｚ２），ｅ３＝（ｘ３,ｙ３、ｚ３）とし、中空球を表す関数の１階偏微分ベクトルを（ρ１,ρ２,ρ３）とする。
Figure JPOXMLDOC01-appb-M000013
　そして、ヘッセ行列の固有値を式（１３）に示すように補正する。なお、評価行列の固有値をλ1、λ２、λ３とし、αを所定の係数とする。ここでαは中実球モデル関数の２階偏微分の原点（中実球の中心）から等間隔な２つの位置に表れる応答波形の一方を最も打ち消すように予め設計された重みである。この重みは中実球の半径Ｒとカーネルサイズsに合わせて設計しておく。
Figure JPOXMLDOC01-appb-M000014
　式（１３）において、λ１＋αρ’は各微分方向において中実球モデル関数の中心に対称な２つの応答波形のうち、一方（例えば図５Ａ（Ａ）の左側の応答波形）を打ち消す補正を行ったときのレスポンスを表す。λ１－αρ’はもう一方（例えば図５Ａ（Ａ）の右側の応答波形）を打ち消す補正を行ったときのレスポンスを表す。式（１３）に示すように、固有値λ１’は、２つのレスポンス｜λ１＋αρ’｜、｜λ１－αρ’｜の小さいほうの値となるように補正されている。この結果、２つのレスポンス｜λ１＋αρ’｜、｜λ１－αρ’｜のいずれかの値が小さい場合には、各固有値λ１’の値が小さくなり、２つのレスポンス｜λ１＋αρ’｜、｜λ１－αρ’｜がいずれも大きい値となる場合にのみ各固有値λ１’の値が大きくなる。λ２’、λ３’ついても同様である。
　従来の方法では誤って構造物を検出していた場合、つまり、中実球モデル関数の中心に対称な２つの応答波形の一方にのみレスポンスが得られる場合には、中実球モデル関数のいずれか一方の応答波形を打ち消すと、２つのレスポンス｜λ１＋αρ’｜、｜λ１－αρ’｜の少なくとも一方の値が小さくなるため、ｍｉｎ（｜λ１＋αρ’｜、｜λ１－αρ’｜）が小さい値となる。一方、構造物を正しく検出している場合、すなわち、各微分方向について画像中の構造物の輪郭のうち対向する２箇所の輪郭部分（構造物を横断する線分の両端）の両方が中実球モデル関数の２階偏微分の応答波形を示す２つの位置とそれぞれ一致している場合には、中実球モデル関数の中心に対称な２つの応答波形のうち、いずれか一方の応答波形を打ち消しても、もう一方の応答波形に基づいて大きいレスポンスが得られるため、２つのレスポンス｜λ１＋αρ’｜、｜λ１－αρ’｜の最小値が大きくなる。従って、誤って構造物を検出していた場合、つまり、２つの応答波形の一方にのみレスポンスが得られる場合を区別することができるため、正確に構造物を判別することができる。なお、式（１３）はレスポンスの絶対値を用いて補正をする例を示しているが、絶対値を用いないで補正を行ってもよい。また、ｍｉｎ（｜λ１＋αρ’｜、｜λ１－αρ’｜）を用いずにλ１＋αρ’とλ１－αρ’を別々に用いて対象構造の判別を行ってもよい。
　そして、評価部３４は、下記の式（１４）、（１５）において、固有値λ１，λ２，λ３に代えて補正後の固有値λ１’，λ２’，λ３’を使用し、これにより算出した値ＲA，ＲB，ＲCを用いて３次元多重解像度画像Ｍｓｉの各画素における線状構造らしさの評価値Ｌ０（Lineness）および面状構造らしさの評価値Ｐ０（Planeness）を算出する。なお、先述の通り、評価部３４は、補正されたヘッセ行列の固有値λ１’，λ２’，λ３’と固有ベクトルｅ１，ｅ２，ｅ３に基づいて、点状構造、線状構造、面状構造を判別することが可能であるが、必ずしも点状構造、線状構造、面状構造の全てを評価する必要はなく、要求される仕様に応じて、点状構造、線状構造、面状構造の一部のみを判別してもよい。本実施形態では、医用画像から線状構造および面状構造を抽出することを目的としているため、線状構造および面状構造についてのみ評価値を算出している。
Figure JPOXMLDOC01-appb-M000015
　なお、式（１４）、（１５）におけるａ～ｈは定数である。また、ＲA，ＲB，ＲCは下記の式（１６）～（１９）により算出する。また、Ｓ2ndは２階偏微分値のパワーであり、下記の式（１９）により算出する。
Figure JPOXMLDOC01-appb-M000016
　図６Ａ、６Ｂを用いて、線状構造物を表すサンプル画像に、フィルタリング部３２によるフィルタリング処理の後に補正部３３による補正処理を施したレスポンスの例を以下に示す。図６Ａは、構造物のサイズ（構造物を横断する線分の長さ）と中実球モデル関数のサイズ（中実球の直径）がほぼ一致する場合のレスポンスを説明する図である。図６Ａ（Ａ）は、線状構造物の輝度値のｘｙ平面図を示す。図６Ａ（Ｂ）は、図６Ａ（Ａ）の画像データに中実球モデル関数のｘ方向の２階偏微分を用いたフィルタリング処理を行って得られたレスポンスをｘｙ平面図で示す。図６Ａ（Ｃ）は図６Ａ（Ａ）の画像データに中空球モデル関数のｘ方向の１階偏微分を用いたフィルタリング処理を行って得られたレスポンスをｘｙ平面図で示す。図６Ａ（Ｄ）は、図６Ａ（Ａ）の画像に中実球モデル関数のｘ方向の２階偏微分行列を用いたフィルタリング処理を行って得られたレスポンスを中空球のｘ方向の１階偏微分のフィルタリング処理を行ったレスポンスを用いて補正したものをｘｙ平面図で表す。
　図６Ａ（Ｄ）に示されるように、線状構造物の中心付近に線状構造らしさのピークが表れており、本実施形態による画像処理によって線状構造物を判別するための良好なレスポンスが得られていることが分かる。
　図６Ｂは、ｘ方向の構造物のサイズ（構造物を横断する線分の長さ）が中実球モデル関数の中実球の直径よりも大きい場合のレスポンスを説明する図である。図６Ｂ（Ａ）は、線状構造物の輝度値のｘｙ平面図を示す。図６Ｂ（Ｂ）は、図６Ｂ（Ａ）画像データに中実球モデル関数のｘ方向の２階偏微分を用いたフィルタリング処理を行って得られたレスポンスをｘｙ平面図で示す。図６Ｂ（Ｃ）は、中空球モデル関数のｘ方向の１階偏微分を用いたフィルタリング処理を行って得られたレスポンスをｘｙ平面図で示す。図６Ｂ（Ｄ）は、図６Ｂ（Ａ）の画像データに中実球モデル関数のｘ方向の２階偏微分行列を用いたフィルタリング処理を行って得られたレスポンスを中空球のｘ方向の１階偏微分のフィルタリング処理を行ったレスポンスを用いて補正したものをｘｙ平面図で表す。
　非特許文献２の手法に対応する評価値を示す図６Ｂ（Ｂ）では中実球モデル関数の２階偏微分による誤判別の原因となる一定のレスポンスが２箇所に表れている（２つの負のピーク）のに対し、本実施形態の手法によるレスポンスを示す図６Ｂ（Ｄ）では、レスポンスが非常に弱くなり、誤判別の原因となるレスポンスが抑制されていることが分かる。
　なお、上記図６Ａ（Ｄ）、図６Ｂ（Ｄ）は、式（１３）のλ１に代えて、式（１）に示すヘッセ行列の要素（Ｉｘｘ）を適用した場合のレスポンスを示している。固有値を算出する処理は、ヘッセ行列における空間を線形変換しているだけでレスポンスの値に影響を与えるものではないため、上記式（１３）において、ヘッセ行列を固有値解析して得られた固有値λ１に替えて、固有値解析をする前のヘッセ行列の要素を補正しても、補正された要素において誤検出の原因となるレスポンスが打ち消されている効果が観察できるからである。
　なお、本実施形態においては、解像度が異なる多重解像度３次元画像Ｍｓｉにおいて、それぞれ線状構造らしさの評価値Ｌ０および面状構造らしさの評価値Ｐ０が算出される。算出された評価値Ｌ０，Ｐ０は、元の３次元画像Ｍ０の対応する画素位置の評価値となるが、すべての多重解像度３次元画像Ｍｓｉの対応する画素位置において評価値が算出されることとなる。
　分割部４０は、判別部３０が算出した線状構造らしさの評価値Ｌ０および面状構造らしさの評価値Ｐ０に基づいて、３次元画像Ｍ０の血管領域および血管以外の領域を領域分割する。具体的には、血管領域を対象領域、血管領域以外の領域を背景領域に設定し、３次元画像Ｍ０内の全画素位置において所定画素サイズの判別領域を設定し、Graph Cut領域分割方法を用いて、判別領域を対象領域と背景領域とに分割する。なお、Graph Cut領域分割方法は、「Yuri Y. Boykov, Marie-Pierre Jolly, “Interactive Graph Cuts for Optimal Boundary and Region Segmentation of Objects in N-D images”, Proceedings of “International Conference on Computer Vision”, Vancouver, Canada, July 2001 vol.I, p.105-112.」に記載されている。本実施形態では、分割部４０は、本発明者による過去出願である特開２０１１－２０６５３１号公報に記載されたものと同じ方法により、判別部３０が算出した線状構造らしさの評価値Ｌ０および面状構造らしさの評価値Ｐ０に基づいて、３次元画像Ｍ０の血管領域および血管以外の領域を領域分割する。
　表示部５０は、２次元画像または３次元画像等を表示するモニタ、ＣＲＴ画面等である。本実施形態においては、表示部５０に、対象領域として分割された線状構造をボリュームレンダリング表示することにより線状構造や面状構造の全体を概観し、その連続性を可視化することができる。
　入力部６０は、キーボード、マウス等である。
　次いで、本実施形態において行われる処理について説明する。図７は本実施形態において行われる処理を示すフローチャートである。まず、画像取得部１０が、Ｘ線ＣＴ装置２により撮像された２次元画像から３次元画像Ｍ０を生成する（ステップＳＴ１）。次に、検出領域設定部２０が、３次元画像Ｍ０を等方化するとともに多重解像度変換して解像度が異なる複数の３次元多重解像度画像Ｍｓｉを生成する（ステップＳＴ２）。
　次いで、判別部３０のフィルタリング部３２が、３次元多重解像度画像Ｍｓｉ毎に、中実球モデル関数の２階偏微分およびガウシアンカーネル関数を用いたフィルタリングを行い、各画素位置におけるヘッセ行列を算出する（ステップＳＴ３）。次いで補正部３３が、上記ヘッセ行列の固有値を中実球モデル関数の１階偏微分ベクトルを用いて補正する（ステップＳＴ４）。次いで、評価部３４が、補正された固有値および固有ベクトルに基づいて線状構造らしさの評価値Ｌ０および面状構造らしさの評価値Ｐ０を算出する（ステップＳＴ５）。
　そして、分割部４０が上述したGraph Cut領域分割方法を用いて、３次元画像Ｍ０を対象領域（血管領域）と背景領域とに分割する（ステップＳＴ６）。そして、表示部５０が分割された対象領域および背景領域をボリュームレンダリング表示し（ステップＳＴ７）、処理を終了する。
　本実施形態の画像処理を実際の患者の胸部ＣＴ画像の血管判別に適用した例と従来の方法を同一のＣＴ画像の血管判別にそれぞれ適用した例を図８Ａ、８Ｂに示す。図８Ａ、８Ｂにおいて、左が従来手法による画像処理の適用例（比較例）であり、右が本実施形態の画像処理の適用例である。図８Ａは、判別された血管領域をボリュームレンダリング法により表した擬似３次元画像を示し、図８Ｂは、図８Ａの一部を拡大したアキシャル断層像で表したものである。
　図８Ａに示すように従来の手法では、多数の誤判別による血管が表示されているところ、本実施形態の適用例においては大幅に誤判別が抑制されて正確に血管が抽出されていることが分かる。また、図８Ｂの比較例（左図）では、血管の内側に線状構造物らしいと誤って判別された部分が薄い灰色で示されているが、本実施形態の適用例（右図）では、血管の内側に線状構造物らしいと判断された部分（薄い灰色の部分）がないことが分かる。
　このように、本実施形態によれば、中実球モデル関数の２階偏微分の１つの応答波形と同じ位置にこの１つの応答波形と略同形状であるとともに正負の符号を逆にした応答波形を有する関数の１階偏微分ベクトルを算出し、この１階偏微分ベクトルをヘッセ行列の固有ベクトルの方向に射影した値を用いて、中実球モデル関数の各方向の２階偏微分の、中実球の中心から対称に離間した２つの位置に表れる応答波形のうち、１つの位置に表れる応答波形を打ち消すように補正したため、構造物の輪郭と中実球モデル関数のいずれかの方向の２階偏微分の応答波形の表れる１つの位置のみが一致した場合に生じる誤判別を抑制でき、評価値の精度を向上することができる。このため、評価値に基づいて画像に含まれる構造物をより正確に判別することができる。一般にフーリエ空間における画像のフィルタリング処理において、フィルタ係数の一部の好ましくない値のみを変更する補正を行うことは容易でないが、上記本実施形態のように中実球モデル関数の２階偏微分のレスポンス特性に着目して中空球モデル関数の１階偏微分のレスポンス特性を利用することにより、良好に非特許文献１の手法による誤判別の問題を解消することができる。
　また、第１の実施形態においては中空球モデル関数の「１階」偏微分ベクトルを補正に用いているため、一旦ｘ、ｙ、ｚ方向の偏微分値を計算すると、任意の方向の偏微分値を計算することができる。すなわち、中空球モデル関数の「１階」偏微分ベクトルを２階偏微分の固有ベクトルの方向に射影して用いることにより、任意の方向の２階偏微分の補正が可能である。結果として、補正のための計算処理の増加が比較的小さく、様々な処理能力の計算機を本実施形態の画像処理装置として使用可能である。
　また、中空球モデル関数の各方向の１階偏微分の１つの位置の応答波形が、中実球モデル関数の各方向の２階偏微分の１つの位置の応答波形と、正負の符号が逆であるとともに非常に近似した形状であるため、誤判別の原因となる中実球モデル関数の２階偏微分の１つの位置の応答波形を好適に打ち消すことができる。
　また、上記第１の実施形態においては、複数サイズの半径Ｒの中実球を表す関数について、ヘッセ行列をそれぞれ算出し、この複数のヘッセ行列に基づいて、それぞれ誤判別の原因となるレスポンスを抑制する補正を行って、各画素位置における線状構造らしさおよび面状構造らしさの評価値を算出している。このため、各中実球サイズに対応するサイズの構造物を好適に評価することができる。また、この評価値に基づいて、上記の第１の実施形態のように対象領域を分割した場合には、好適に複数サイズの構造物を対象領域として分割することができる。
　なお、異なるサイズの構造物に対してヘッセ解析を行う手法として、本実施形態のように複数の解像度の画像に対して、一定のカーネルサイズσおよび一定の中実球モデル関数の半径Ｒを適用してヘッセ解析を行う手法を適用してもよく、一定の解像度の画像に対して、カーネルサイズσおよび中実球モデル関数の半径Ｒを異ならせてヘッセ解析を行う手法を適用してもよい。例えば、構造物の長軸や短軸の長さなど、検出対象構造物を横断する線分の長さを試験的に得られたデータなどから予め取得しておき、検出対象構造物を横断する線分の長さと半径Ｒの大きさが一致するように、３次元多重解像度画像Ｍｓｉの解像度、あるいは、カーネルサイズσおよび中実球の半径Ｒを設定することが考えられる。
　また、上記第１の実施形態においては、線状構造らしさの評価値Ｌ０および面状構造らしさの評価値Ｐ０の双方を算出しているが、線状構造らしさの評価値Ｌ０および面状構造らしさの評価値Ｐ０のいずれか一方のみを算出するようにしてもよい。
　本発明の第２の実施形態を以下に説明する。第２の実施形態においては、補正部３３による補正処理が第１の実施形態と異なるが、それ以外は第１の実施形態と同じであり、各機能ブロックの機能も共通している。以下、第１の実施形態と異なる点を中心に説明し、第１の実施形態と同じ点については説明を省略する。
　先述の中実球モデル関数の２階偏微分による２つの位置における応答波形の１つの位置における応答波形を打ち消すために、第２の実施形態では、中実球モデル関数の２階偏微分のｘ方向の負側に位置する応答波形がｘ方向の負側から１つの正のピークとこれに隣接する１つの負のピークを有するものであり、ｘ方向の正側に位置する応答波形がｘ方向の負側から１つの負のピークとこれに隣接する１つの正のピークを有するものであるという特徴に着目し、この一方の応答波形と同じ位置に位置し、かつ、この一方の応答波形と正負の符号が逆で略同形状の応答波形を、この中実球よりやや小さいサイズの中実球モデル関数ｆ２（ｒ）（第２の中実球モデル関数）の各方向の１階偏微分と、この中実球よりやや大きいサイズの中実球モデル関数ｆ３（ｒ）（第３の中実球モデル関数）の各方向の１階偏微分を組み合わせて作成可能であることを見出した。なお、ここでは、区別のために、フィルタリング部３２がフィルタリング処理に用いる中実球モデル関数を第１の中実球モデル関数と称する。
　図９Ａ（Ａ）は、第２の半径Ｒ２を有する第２の中実球モデル関数ｆ２（ｒ）のｘ方向の１階偏微分のレスポンスを示し、図９Ａ（Ｂ）は、第３の半径Ｒ３を有する第３の中実球モデル関数ｆ３（ｒ）のｘ方向の１階偏微分のレスポンスを示し、図９Ａ（Ｃ）は、正負の符号を反転させた第２の中実球モデル関数ｆ２（ｒ）のｘ方向の１階偏微分に第３の中実球モデル関数ｆ３（ｒ）のｘ方向の１階偏微分を加算したレスポンスを示す。
　図９Ａ（Ａ）、（Ｂ）に示すように、第２および第３の中実球モデル関数ｆ２（ｒ）、ｆ３（ｒ）の１階偏微分のレスポンスは、中実球の中心から負側の位置に１つの正のピークを有する応答波形を有し、中実球の中心から正側の位置に１つの負のピークを有する応答波形を有し、各応答波形は中心から対称に表れるという特徴を有する。
　そして、図９Ａ（Ａ）に示す中実球モデル関数ｆ２（ｒ）のｘ方向の１階偏微分の正負の符号を反転させたものと、図９Ａ（Ｂ）に示す中実球モデル関数ｆ３（ｒ）のｘ方向の１階偏微分を組み合わせて用いることにより、図９Ａ（Ｃ）に示すようなｘ方向のレスポンスを実現できる。図９Ａ（Ｃ）は、中実球の中心から対称な２つの位置に応答波形を有するものであり、ｘ方向の負側に位置する応答波形はｘ方向の負側から１つの負のピークとこれに隣接する１つの正のピークを有するものであり、ｘ方向の正側に位置する応答波形はｘ方向の負側から１つの負のピークとこれに隣接する１つの正のピークを有するという特徴を有するものである。
　つまり、中実球モデル関数の各方向の２階偏微分行列によるフィルタリングに、第２および第３の中実球モデル関数ｆ２（ｒ）、ｆ３（ｒ）の各方向の１階偏微分を図９Ａ（Ｃ）のようなレスポンスを組み合わせると、各方向の正側の応答波形は、同じ位置において逆の符号を有するため打ち消し合い、ｘ方向の負側の応答波形は、同じ位置において同じ符号を有するため強め合うため、図５Ａ（Ｃ）に示すようなｘ方向の正側の１つの位置にのみ応答波形を有するものとなる。
　なお、本実施形態では、第１の中実球モデル関数の２階偏微分の１つの位置の応答波形の打ち消しを行うために、第２の半径Ｒ２は、第１の中実球モデル関数の２階偏微分における中実球の中心から正方向に位置する応答波形の正のピークまでの距離と同じであり、第３の半径Ｒ３は、第１の中実球モデル関数の２階偏微分における中実球の中心から正方向に位置する応答波形の負のピークまでの距離と同じである。この第２の半径Ｒ２（第３の半径Ｒ３）は、上記中実球モデル関数の各方向の２階偏微分の１つの位置における応答波形を打ち消す効果を有する範囲で、第２の半径Ｒ２（第３の半径Ｒ３）より「中心から応答波形の正（負）のピークまでの距離」が大きくてもよく、小さくてもよいが、中実モデル関数の２階偏微分の応答波形と、第２および第３の中実球モデル関数ｆ２（ｒ）、ｆ３（ｒ）の各方向の１階偏微分を組み合わせた応答波形をより対応したものとするために、本実施形態のように第２の半径Ｒ２（第３の半径Ｒ３）が「中心から応答波形の正（負）のピークまでの距離と厳密に一致することが好ましい。
　ｘ方向における第２および第３の中実球モデル関数の１階偏微分のレスポンスを図９Ｂを用いて説明する。図９Ｂ（Ａ）は、第２の中実球モデル関数ｆ２（ｒ）のｘ方向の１階偏微分のレスポンスを示し、図９Ｂ（Ｂ）は第３の中実球モデル関数ｆ３（ｒ）のｘ方向の１階偏微分のレスポンスを示し、図９Ｂ（Ｃ）は、正負の符号を反転させた中実球モデル関数ｆ２（ｒ）のｘ方向の１階偏微分に中実球モデル関数ｆ３（ｒ）のｘ方向の１階偏微分を加算した場合のレスポンスを示す。なお、図９Ｂにおいて、各上述のレスポンスを等間隔にｚ座標値を異ならせた複数のｘｙ平面図により示しており、これらの複数のｘｙ平面図をｚ座標の大きいものから垂直方向に並べて示している。図９Ｂにおける各ｘｙ平面図において、明るい（白い）ほどレスポンスが正の方向に大きくなり、暗い（黒い）ほどレスポンスが負の方向に大きくなる。
　図９Ｂ（Ｃ）に示す負側の応答波形は、図９Ｂ（Ａ）に示す第２の中実球モデル関数のｘ方向の１階偏微分の中心から負側に位置する負のピークと、図９Ｂ（Ｂ）に示す第３の中実球モデル関数のｘ方向の１階偏微分の中心から負側に位置する負のピークがずれて重なり合い、負の方向から負のピークと正のピークが隣接する波形となっていることが分かる。また、図９Ｂ（Ｃ）に示す正側の応答波形は、図９Ｂ（Ａ）に示す第２の中実球モデル関数の１階偏微分の中心から正側に位置する正のピークと、図９Ｂ（Ｃ）に示す第３の中実球モデル関数の１階偏微分の中心から正側に位置する正のピークがずれて重なり合い、負の方向から負のピークと正のピークが隣接する波形となっていることが分かる。第２および第３の中実球モデル関数の１階偏微分を図９Ｂ（Ｃ）のレスポンスを示すように組み合わせて、図５Ｃ（Ｂ）に示す第１の中実球モデル関数のｘ方向の２階偏微分を補正すると、図５Ｄに示すように中実球の中心から負側の応答波形を打ち消すことができる。同様に各方向について中実球の中心から負側の応答波形を打ち消すことができる。
　本第２の実施形態では、上記原理に基づいて、補正部３３は、第２および第３の中実球モデル関数の各１階偏微分を用いて第１の中実球モデル関数ｆ（ｒ）の１つの位置の応答波形のみを打ち消すようにヘッセ行列を補正する。具体的な補正方法を以下に説明する。
　まず、補正部３３は、第２の半径Ｒ２の中実球モデル関数について、式（２）において半径Ｒに替えて第２の半径Ｒ２を用いることにより第２の中実球モデル関数の１階偏微分ベクトルを算出する。なお、式（２）では、（２πνｘ）ｌ×（２πνｙ）ｍ×（２πνｘ）ｎの部分はフーリエ空間における微分処理に相当し、ｌ＋ｍ＋ｎ＝１（０＜ｌ、０＜ｍ、０＜ｎ）となるように、各微分方向に対応する係数ｌ、ｍ、ｎをそれぞれ式（９）に代入することにより、第２の中実球モデル関数ｆ２（ｒ）について算出した各１階偏微分ベクトル（ｓｘ、ｓｙ、ｓｚ）をそれぞれ算出できる。また、同様に、補正部３３は、第３の中実球モデル関数ｆ３（ｒ）についても、式（２）において半径Ｒに替えて第３の半径Ｒ３を用いることにより１階偏微分ベクトル（ｔｘ、ｔｙ、ｔｚ）を算出する。そして、第２の中実球モデル関数ｆ２（ｒ）について算出した各１階偏微分ベクトル（ｓｘ、ｓｙ、ｓｚ）、および第３の中実球モデル関数ｆ３（ｒ）について算出した各１階偏微分ベクトル（ｔｘ、ｔｙ、ｔｚ）を重み付け加算して、（ρ１、ρ２、ρ３）＝（ｓｘ－βｔｘ、ｓｙ－βｔｙ、ｓｚ－βｔｚ）を算出する。
　なお、式（２）においては、フーリエ空間においてフィルタリングを行っているが、実空間においてフィルタリングを行ってもよい。
　上述したように算出したＸ方向、Ｙ方向およびＺ方向の１階偏微分ベクトル（ρ１,ρ２,ρ３）は、固有値λ１，λ２，λ３の固有ベクトルｅ１，ｅ２，ｅ３の方向とずれているため、補正部３３は、固有ベクトルｅ１，ｅ２，ｅ３の方向に対応した１階偏微分ベクトルρ１ ’,ρ２ ’,ρ３ ’を下記の式（１２）により算出する。
　そして、第１の実施形態同様に、式（１３）に基づいてヘッセ行列の固有値を補正する。なお、評価行列の固有値をλ1、λ２、λ３とし、α、βを所定の係数とする。ここでα、βは中実球モデル関数の２階偏微分の原点（中実球の中心）から等間隔な２つの位置に表れる応答波形の一方を最も打ち消すように予め設計された重みである。この重みは中実球の半径Ｒとカーネルサイズsに合わせて設計しておく。
　そして、評価部３４は、第１の実施形態と同様に、式（１４）、（１５）において、固有値λ１，λ２，λ３に代えて固有値λ１’，λ２’，λ３’を使用し、これにより算出した値ＲA，ＲB，ＲCを用いて３次元多重解像度画像Ｍｓｉの各画素における線状構造らしさの評価値Ｌ０（Lineness）および面状構造らしさの評価値Ｐ０（Planeness）を算出する。
　第２の実施形態においても、第１の実施形態同様に、中実球モデル関数の２階偏微分の１つの応答波形と同じ位置にこの１つの応答波形と略同形状であるとともに正負の符号を逆にした応答波形を有する関数（２種類のサイズの中実球モデル関数の１階偏微分を組み合わせた関数）の１階偏微分ベクトルを算出し、この１階偏微分ベクトルをヘッセ行列の固有ベクトルの方向に射影した値を用いて、中実球モデル関数の各方向の２階偏微分の、中実球の中心から対称に離間した２つの位置に表れる応答波形のうち、１つの位置に表れる応答波形を打ち消すように補正している。このため、第２の実施形態では、２種類のサイズの中実球モデル関数のそれぞれの１階偏微分を補正に用いる分、１種類のサイズの中空球モデル関数の１階偏微分を補正に用いる第１の実施形態の手法よりも補正部３３による補正処理の計算量は若干増加するものの、第１の実施形態同様に構造物の輪郭と中実球モデル関数のいずれかの方向の２階偏微分の応答波形の表れる１つの位置のみが一致した場合に生じる誤判別を抑制でき、評価値の精度を向上することができる。
　また、第２の実施形態において、第２の半径Ｒ２は、中実球の中心から正のピークまでの長さと中実球の中心から負のピークまでの長さのうち、長い方の長さと一致するものであり、第３の半径Ｒ３は、中実球の中心から正のピークまでの長さと中実球の中心から負のピークまでの長さのうち、短い方の長さと一致するものであるため、応答波形を打ち消す効果が、第１の中実球の正のピークと負のピークが隣接した一方の応答波形により対応したものとなり、好適に誤判別の原因となる応答波形の抑制を実現できる。
　また、上記第２の実施形態においても、複数サイズの半径Ｒの第１の中実球を表す関数について、ヘッセ行列をそれぞれ算出し、この複数のヘッセ行列に基づいて、それぞれ誤判別の原因となるレスポンスを抑制する補正を行って、各画素位置における線状構造らしさおよび面状構造らしさの評価値を算出することが好ましい。この場合には、各カーネルサイズに対応するサイズの構造物を好適に評価することができる。またこのために、複数の解像度の画像に対して一定の半径Ｒの第１の中実球モデル関数を用いて第２の実施形態の画像処理を行ってもよく、１種類の解像度の画像に対して、第１の半径Ｒを異ならせて、複数サイズの第１の半径Ｒの第１の中実球モデル関数を用いて第２の実施形態の画像処理を行ってもよい。
　上記各実施形態による画像処理方法は、以下のような方向判別のための画像処理装置にも適用可能である。例えば、方向判別のための画像処理装置は、第１の実施形態における分割部４０と評価部３４と、表示部５０を備えておらず、代わりに方向判別部を備えたものとできる。このような画像処理装置の一例として、画像取得部１０による画像取得処理と、検出領域設定部２０による検出領域設定処理と、フィルタリング部３２によるフィルタリング処理については、第１の実施形態と同じであり、各機能ブロックの機能も共通したものとして、方向判別部を以下のように構成できる。
　方向判別部は、フィルタリング部３２が算出した評価行列を固有値解析して固有値および固有ベクトルｅ１，ｅ２，ｅ３を算出し、算出された固有ベクトルの向きに基づいて線状構造物の主軸方向または面状構造の法線方向を判別する。詳細には、フィルタリング部３２が得たヘッセ行列の固有値λ１，λ２，λ３および固有ベクトルｅ１，ｅ２，ｅ３を取得し、式（６）に示す関係を満たしていれば、固有値λ１に対応する固有ベクトルｅ１の指す方向を線状構造の軸方向と判別し、式（７）に示す関係を満たしていれば、固有値λ３に対応する固有ベクトルｅ３の指す方向を面状構造の法線方向と判別する。なお、方向判別部は、ヘッセ行列の固有値と固有ベクトルに基づいて、式（６）に基づく線状構造の軸方向、または、式（７）に基づく面状構造の法線方向のいずれか一方のみを判別してもよい。
　上記のような方向を判別するための画像処理装置により得られた線状構造の軸方向や面状構造の法線方向は、例えば血管や大腸など管状構造のＣＰＲ(Curved Planer Reformation / Reconstruction)法によって３次元医用画像から再構成されたＣＰＲ画像を作成する際の軸方向として用いることができ、線状構造の軸方向や面状構造の法線方向を必要とする周知の種々の処理に好適に用いることができる。なお、本発明により得られた線状構造の軸方向や面状構造の法線方向は、例えば、「Armin Kanitsar, Dominik Fleischmann, Rainer Wegenkittl, Petr Felkel, Eduard Groller: CPR - Curved Planar Reformation. IEEE Visualization 2002.」など周知の種々のＣＰＲ画像生成手法に用いることができる。
　また、上記各実施形態においては、画像が医用画像であり、構造物が血管であるため、様々な直径の血管を含むために多数誤判別が生じていた医用画像の構造物判別の精度を顕著に向上させることができる。
　また、上記各実施形態においては、フィルタリング部は、中実球を表す関数の２階偏微分行列を用いたフィルタリング、中実球を表す関数の１階偏微分行列を用いたフィルタリングおよび中空球を表す関数の１階偏微分行列を用いたフィルタリングをフーリエ空間において実施したため、実空間によりフィルタリング処理を行う場合よりも処理コストおよび処理速度が低減でき好ましい。ただし、フィルタリング部は、中実球を表す関数の２階偏微分行列を用いたフィルタリング、中実球を表す関数の１階偏微分行列を用いたフィルタリングおよび中空球を表す関数の１階偏微分行列を用いたフィルタリングを、実空間で行ってもよい。
　また、上記各実施形態に示したように、評価部は、固有ベクトルおよび補正された固有値に基づいて、構造物の局所的な点状構造、線状構造、面状構造の少なくとも１つを好適に判別することができる。ただし、これらの固有ベクトルや固有値は、点状構造、線状構造、面状構造だけでなく、あらゆる形状の構造物の判別に用いてもよい。
　なお、上記各実施形態では線状構造物として血管の判別を例に挙げたが、気管支等の他の線状構造物の判別にも適用できる。また、また骨のみならず、皮膚、葉間胸膜等の面状構造物の判別にも利用することができる。
　なお、上記各実施形態においては、３次元画像Ｍ０に含まれる線状構造物および面状構造物を対象に各評価処理または方向判別処理を行っているが、２次元画像に含まれる線状構造物および面状構造物を対象に各評価処理または方向判別処理を行ってもよい。
　また、上記各実施形態においては、Graph Cut領域分割方法を用いて３次元画像Ｍ０に含まれる線状構造物と面状構造物とを分割しているが、Watershedアルゴリズム等の他の領域分割手法を用いてもよいことはもちろんである。ここで、Watershedアルゴリズムは、画像の画素値情報を高度とみなした地形において水を満たしていく際に、異なるくぼみに溜まる水の間で境界ができるように画像を分割する手法である。このため、３次元画像Ｍ０上の評価値Ｌ０，Ｐ０に対し、適当な平滑化を行った後でWatershedアルゴリズムを実行することによって、線状構造物および面状構造物を分割することが可能である。
　上記の各実施形態はあくまでも例示であり、上記のすべての説明が本発明の技術的範囲を限定的に解釈するために利用されるべきものではない。
　この他、上記の実施形態におけるシステム構成、ハードウェア構成、処理フロー、モジュール構成、ユーザインターフェースや具体的処理内容等に対して、本発明の趣旨から逸脱しない範囲で様々な改変を行ったものも、本発明の技術的範囲に含まれる。

Claims (12)

  1. 　画像中の各画素位置において、中実球を表す関数の２階偏微分を用いたフィルタリングを施してヘッセ行列を算出するフィルタリング部と、
    　前記算出されたヘッセ行列を固有値解析して得られた固有値および固有ベクトルを用いて、画像中に含まれる構造物を判別する評価部とを備え、
    　前記フィルタリング部が、前記画像中の前記各画素位置において、前記中実球の半径と同じ半径を有する中空球を表す関数の１階偏微分を用いたフィルタリングを施して１階偏微分ベクトルを取得し、前記取得した１階偏微分ベクトルを前記固有ベクトルの方向に射影した値を用いて、前記中実球を表す関数の各方向の２階偏微分の、前記中実球の中心から対称に離間した２つの位置に表れる応答波形のうち一方の応答波形を打ち消す補正を行う補正部を備えたものであることを特徴とする画像処理装置。
  2. 　前記評価行列の固有値をλ1、λ２、λ３とし、固有ベクトルｅ１＝（ｘ１,ｙ１、ｚ１），ｅ２＝（ｘ２,ｙ２、ｚ２），ｅ３＝（ｘ３,ｙ３、ｚ３）とし、前記中空球を表す関数の前記１階偏微分ベクトルを（ρ１,ρ２,ρ３）とすると、前記補正部は、下記式（１２）で算出されるρ１ ’,ρ２ ’,ρ３ ’および所定の係数αを用いて、下記式（１３）に示すように、前記固有値を補正することにより、前記一方の応答波形を打ち消す補正を行うものであることを特徴とする請求項１記載の画像処理装置。
    Figure JPOXMLDOC01-appb-M000017
    Figure JPOXMLDOC01-appb-M000018
  3. 　前記中空球を表す関数は、下記式（１０）により表されるものであることを特徴とする請求項１または２記載の画像処理装置。
    Figure JPOXMLDOC01-appb-M000019
     ここで、ｘ、ｙ、ｚは３次元空間の座標、ｒはその極座標表現であり、Ｒ４は前記中空球の半径とする。
  4. 　画像中の各画素位置において、中実球を表す関数の２階偏微分を用いたフィルタリングを施してヘッセ行列を算出するフィルタリング部と、
    　前記算出されたヘッセ行列を固有値解析して得られた固有値および固有ベクトルを用いて、画像中に含まれる構造物を判別する評価部とを備え、
    　前記フィルタリング部が、前記画像中の各画素位置において、前記第１の中実球の半径である第１の半径より大きい第２の半径を有する第２の中実球を表す関数の１階偏微分を用いたフィルタリングを施して１階偏微分ベクトルを算出し、かつ、前記画像中の各画素位置において、前記第１の半径より小さい第３の半径を有する第３の中実球を表す関数の１階偏微分を用いたフィルタリングを施して１階偏微分ベクトルを算出し、前記第２の中実球を表す関数の１階偏微分ベクトルを前記固有ベクトルの方向に射影した値および第３の中実球を表す関数の１階偏微分ベクトルを前記固有ベクトルの方向に射影した値を用いて、前記中実球を表す関数の各方向の２階偏微分の、前記中実球の中心から対称に離間した２つの位置に表れる応答波形のうち一方の位置の応答波形を打ち消す補正を行う補正部を備えたものであることを特徴とする画像処理装置。
  5. 　前記一方の位置の応答波形は、１つの正のピークと１つ負のピークが隣接する波形であり、
    　前記第２の半径は、前記中実球の中心から前記正のピークまでの長さと前記中実球の中心から前記負のピークまでの長さのうち、長い方の長さと一致するものであり、
    　前記第３の半径は、前記中実球の中心から前記正のピークまでの長さと前記中実球の中心から前記負のピークまでの長さのうち、短い方の長さと一致するものであることを特徴とする請求項４に記載の画像処理装置。
  6. 　前記フィルタリング部は、複数のサイズの前記中実球を表す関数について、該各中実球を表す関数の２階偏微分行列によるフィルタリングを施した評価行列をそれぞれ算出するものであることを特徴とする請求項１から３のいずれか１項記載の画像処理装置。
  7. 　前記フィルタリング部は、複数のサイズの前記第１の中実球を表す関数について、該各第１の中実球を表す関数の２階偏微分行列によるフィルタリングを施した評価行列をそれぞれ算出するものであることを特徴とする請求項４または５のいずれか１項記載の画像処理装置。
  8. 　前記フィルタリング部は、前記中実球を表す関数の２階偏微分行列を用いたフィルタリングを、フーリエ空間において実施するものであることを特徴とする請求項１から７のいずれか１項記載の画像処理装置。
  9. 　前記評価部は、前記構造物の局所的な点状構造、線状構造、面状構造の少なくとも１つを判別するものであることを特徴とする請求項１から８のいずれか1項記載の画像処理装置。
  10. 　前記画像が医用画像であり、前記構造物が血管であることを特徴とする請求項１から９のいずれか１項記載の画像処理装置。
  11. 　画像中の各画素位置において、中実球を表す関数の２階偏微分を用いたフィルタリングを施してヘッセ行列を算出するフィルタリング工程と、
    　前記算出されたヘッセ行列を固有値解析して得られた固有値および固有ベクトルを用いて、画像中に含まれる構造物を判別する評価工程とを備えた画像処理方法であって、
    　前記フィルタリング工程が、前記画像中の前記各画素位置において、前記中実球の半径と同じ半径を有する中空球を表す関数の１階偏微分を用いたフィルタリングを施して１階偏微分ベクトルを取得し、前記取得した１階偏微分ベクトルを前記固有ベクトルの方向に射影した値を用いて、前記中実球を表す関数の各方向の２階偏微分の、前記中実球の中心から対称に離間した２つの位置に表れる応答波形のうち一方の応答波形を打ち消す補正を行うものであることを特徴とする画像処理方法。
  12. 　コンピュータに、
    　画像中の各画素位置において、中実球を表す関数の２階偏微分を用いたフィルタリングを施してヘッセ行列を算出するフィルタリング工程と、
    　前記算出されたヘッセ行列を固有値解析して得られた固有値および固有ベクトルを用いて、画像中に含まれる構造物を判別する評価工程とを実行させる画像処理プログラムであって、
    　　前記フィルタリング工程が、前記画像中の前記各画素位置において、前記中実球の半径と同じ半径を有する中空球を表す関数の１階偏微分を用いたフィルタリングを施して１階偏微分ベクトルを取得し、前記取得した１階偏微分ベクトルを前記固有ベクトルの方向に射影した値を用いて、前記中実球を表す関数の各方向の２階偏微分の、前記中実球の中心から対称に離間した２つの位置に表れる応答波形のうち一方の応答波形を打ち消す補正を行うものであることを特徴とする画像処理プログラム。
PCT/JP2013/001636 2012-03-14 2013-03-13 画像処理装置および方法並びにプログラム WO2013136784A1 (ja)

Priority Applications (4)

Application Number Priority Date Filing Date Title
CN201380014282.4A CN104168820B (zh) 2012-03-14 2013-03-13 图像处理装置及方法
EP13760382.5A EP2826415B1 (en) 2012-03-14 2013-03-13 Image processing device, method, and program
US14/484,046 US9183637B2 (en) 2012-03-14 2014-09-11 Image processing apparatus, method, and program
IN8510DEN2014 IN2014DN08510A (ja) 2012-03-14 2014-10-11

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2012056855A JP5833958B2 (ja) 2012-03-14 2012-03-14 画像処理装置および方法並びにプログラム
JP2012-056855 2012-03-14

Related Child Applications (1)

Application Number Title Priority Date Filing Date
US14/484,046 Continuation US9183637B2 (en) 2012-03-14 2014-09-11 Image processing apparatus, method, and program

Publications (1)

Publication Number Publication Date
WO2013136784A1 true WO2013136784A1 (ja) 2013-09-19

Family

ID=49160712

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2013/001636 WO2013136784A1 (ja) 2012-03-14 2013-03-13 画像処理装置および方法並びにプログラム

Country Status (6)

Country Link
US (1) US9183637B2 (ja)
EP (1) EP2826415B1 (ja)
JP (1) JP5833958B2 (ja)
CN (1) CN104168820B (ja)
IN (1) IN2014DN08510A (ja)
WO (1) WO2013136784A1 (ja)

Families Citing this family (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2017019059A1 (en) * 2015-07-29 2017-02-02 Perkinelmer Health Sciences, Inc. Systems and methods for automated segmentation of individual skeletal bones in 3d anatomical images
US10977764B2 (en) 2015-12-29 2021-04-13 Dolby Laboratories Licensing Corporation Viewport independent image coding and rendering
KR102139801B1 (ko) * 2017-09-21 2020-07-30 (주)릴리커버 피부진단방법, 피부진단장치 및 기록매체
WO2019069867A1 (ja) 2017-10-03 2019-04-11 株式会社根本杏林堂 血管抽出装置および血管抽出方法
US10573060B1 (en) * 2018-06-14 2020-02-25 Kilburn Live, Llc Controller binding in virtual domes
US10740957B1 (en) * 2018-06-14 2020-08-11 Kilburn Live, Llc Dynamic split screen
CN109998681B (zh) * 2019-03-16 2022-02-01 哈尔滨理工大学 一种区分镜面反射区域和血管的内腔图像预处理方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20110064289A1 (en) * 2009-09-14 2011-03-17 Siemens Medical Solutions Usa, Inc. Systems and Methods for Multilevel Nodule Attachment Classification in 3D CT Lung Images
JP2011104206A (ja) * 2009-11-19 2011-06-02 Hiroshima City Univ 動脈瘤候補の検出支援装置および検出方法
JP2011206531A (ja) 2010-03-11 2011-10-20 Fujifilm Corp 画像処理装置および方法並びにプログラム

Family Cites Families (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6448968B1 (en) * 1999-01-29 2002-09-10 Mitsubishi Electric Research Laboratories, Inc. Method for rendering graphical objects represented as surface elements
US6342886B1 (en) * 1999-01-29 2002-01-29 Mitsubishi Electric Research Laboratories, Inc Method for interactively modeling graphical objects with linked and unlinked surface elements
US9451932B2 (en) * 2004-12-30 2016-09-27 Crystalview Medical Imaging Limited Clutter suppression in ultrasonic imaging systems
CN100357701C (zh) * 2005-07-12 2007-12-26 北京航空航天大学 一种栅格型标定点亚像素提取方法
CN100405004C (zh) * 2006-08-25 2008-07-23 北京航空航天大学 光条图像特征高精度快速提取装置及方法
JP5161845B2 (ja) * 2009-07-31 2013-03-13 富士フイルム株式会社 画像処理装置及び方法、データ処理装置及び方法、並びにプログラム
CN101877063A (zh) * 2009-11-25 2010-11-03 中国科学院自动化研究所 一种基于亚像素级特征点检测的图像匹配方法
US8724428B1 (en) * 2012-11-15 2014-05-13 Cggveritas Services Sa Process for separating data recorded during a continuous data acquisition seismic survey

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20110064289A1 (en) * 2009-09-14 2011-03-17 Siemens Medical Solutions Usa, Inc. Systems and Methods for Multilevel Nodule Attachment Classification in 3D CT Lung Images
JP2011104206A (ja) * 2009-11-19 2011-06-02 Hiroshima City Univ 動脈瘤候補の検出支援装置および検出方法
JP2011206531A (ja) 2010-03-11 2011-10-20 Fujifilm Corp 画像処理装置および方法並びにプログラム

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
A. F. FRANGI ET AL.: "Multiscale vessel enhancement filtering", PROCEEDINGS OF MICCAI, 1998, pages 130 - 137, XP002345875, DOI: doi:10.1007/BFb0056195
AKIHIRO ODA ET AL.: "Development of a lymph node detection method from 3D abdominal CT images using multi shape ellipsoidal structure detection filter", IEICE TECHNICAL REPORT MI2011-123, January 2012 (2012-01-01), pages 251 - 256, XP008174743 *
ARMIN KANITSAR; DOMINIK FLEISCHMANN; RAINER WEGENKITTL; PETR FELKEL; EDUARD GROLLER: "CPR - Curved Planar Reformation", IEEE VISUALIZATION, 2002
M. LAW ET AL.: "Three Dimensional Curvilinear Structure Detection Using Optimally Oriented Flux", PROCEEDINGS OF ECCV, 2008, pages 368 - 382, XP019109329
YURI Y. BOYKOV; MARIE-PIERRE JOLLY: "Interactive Graph Cuts for Optimal Boundary and Region Segmentation of Objects in N-D images", PROCEEDINGS OF ''INTERNATIONAL CONFERENCE ON COMPUTER VISION, vol. I, July 2001 (2001-07-01), pages 105 - 112
ZHAOYANG JIN ET AL.: "Enhancement of Venous Vasculature in Susceptibility Weighted Images of the Brain Using Multi-Scale Vessel Enhancement Filtering", 3RD INTERNATIONAL CONFERENCE ON BIOMEDICAL ENGINEERING AND INFORMATION (BMEI 2010), 2010, pages 226 - 230, XP031803974 *

Also Published As

Publication number Publication date
IN2014DN08510A (ja) 2015-05-15
EP2826415A4 (en) 2016-02-24
CN104168820A (zh) 2014-11-26
EP2826415A1 (en) 2015-01-21
JP5833958B2 (ja) 2015-12-16
US9183637B2 (en) 2015-11-10
US20150016686A1 (en) 2015-01-15
JP2013191007A (ja) 2013-09-26
CN104168820B (zh) 2016-05-11
EP2826415B1 (en) 2018-01-10

Similar Documents

Publication Publication Date Title
WO2013136784A1 (ja) 画像処理装置および方法並びにプログラム
JP6267710B2 (ja) 医用画像中の肺結節を自動検出するためのシステム及び方法
JP6570145B2 (ja) 画像を処理する方法、プログラム、代替的な投影を構築する方法および装置
Liu et al. Deep iterative reconstruction estimation (DIRE): approximate iterative reconstruction estimation for low dose CT imaging
JP5037705B2 (ja) 画像処理装置および方法並びにプログラム
Hogeweg et al. Suppression of translucent elongated structures: applications in chest radiography
CN101779222A (zh) 对高对比度对象进行的基于投影的去除
WO2017086433A1 (ja) 医用画像処理方法及び装置及びシステム及びプログラム
Oh et al. Learning bone suppression from dual energy chest X-rays using adversarial networks
Lee et al. No-reference perceptual CT image quality assessment based on a self-supervised learning framework
JP4709290B2 (ja) 画像処理装置および方法並びにプログラム
JP5748693B2 (ja) 画像処理装置および方法並びにプログラム
Klinder et al. Lobar fissure detection using line enhancing filters
Czajkowska et al. 4d segmentation of ewing’s sarcoma in MR images
Lecron et al. Points of interest detection in cervical spine radiographs by polygonal approximation
Gill et al. Segmentation of lungs with interstitial lung disease in CT Scans: A TV-L 1 based texture analysis approach
You et al. CT Super-resolution GAN Constrained by the Identical, Residual, and Cycle Learning Ensemble
Shi et al. A combinational filtering method for enhancing suspicious structures in chest X-rays
US20220335602A1 (en) Method and system for image normalisation
Al-Samaraie et al. Medical colored image enhancement using wavelet transform followed by image sharpening
Han et al. Reconstruction of Partially Broken Vascular Structures in X-ray Images via Vesselness-loss-based Multi-scale Generative Adversarial Networks
Schwier et al. An object-based image analysis approach to spine detection in CT images
CN114565711A (zh) 基于深度学习的心脏图像重建方法及系统
Bu et al. Method of CT pulmonary vessel image enhancement based on structure tensor
Ruppert et al. Symmetry-based Method for Anomalies Detection on Morphological Neuroimages

Legal Events

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

Ref document number: 201380014282.4

Country of ref document: CN

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

Ref document number: 13760382

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

WWE Wipo information: entry into national phase

Ref document number: 2013760382

Country of ref document: EP

REG Reference to national code

Ref country code: BR

Ref legal event code: B01A

Ref document number: 112014022824

Country of ref document: BR

ENP Entry into the national phase

Ref document number: 112014022824

Country of ref document: BR

Kind code of ref document: A2

Effective date: 20140915