WO2014156089A1 - 画像処理装置、画像処理プログラムおよび画像処理装置の作動方法 - Google Patents

画像処理装置、画像処理プログラムおよび画像処理装置の作動方法 Download PDF

Info

Publication number
WO2014156089A1
WO2014156089A1 PCT/JP2014/001618 JP2014001618W WO2014156089A1 WO 2014156089 A1 WO2014156089 A1 WO 2014156089A1 JP 2014001618 W JP2014001618 W JP 2014001618W WO 2014156089 A1 WO2014156089 A1 WO 2014156089A1
Authority
WO
WIPO (PCT)
Prior art keywords
contour
energy
image processing
object region
pixels
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Ceased
Application number
PCT/JP2014/001618
Other languages
English (en)
French (fr)
Inventor
嘉郎 北村
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Fujifilm Corp
Original Assignee
Fujifilm Corp
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Fujifilm Corp filed Critical Fujifilm Corp
Priority to DE112014001697.7T priority Critical patent/DE112014001697T5/de
Publication of WO2014156089A1 publication Critical patent/WO2014156089A1/ja
Priority to US14/865,616 priority patent/US9965698B2/en
Anticipated expiration legal-status Critical
Ceased legal-status Critical Current

Links

Images

Classifications

    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/02Arrangements for diagnosis sequentially in different planes; Stereoscopic radiation diagnosis
    • A61B6/03Computed tomography [CT]
    • A61B6/032Transmission computed tomography [CT]
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/50Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment specially adapted for specific body parts; specially adapted for specific clinical applications
    • A61B6/503Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment specially adapted for specific body parts; specially adapted for specific clinical applications for diagnosis of the heart
    • GPHYSICS
    • G06COMPUTING OR CALCULATING; COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00Pattern recognition
    • G06F18/20Analysing
    • G06F18/24Classification techniques
    • GPHYSICS
    • G06COMPUTING OR CALCULATING; COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/10Segmentation; Edge detection
    • G06T7/11Region-based segmentation
    • GPHYSICS
    • G06COMPUTING OR CALCULATING; COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/10Segmentation; Edge detection
    • G06T7/12Edge-based segmentation
    • GPHYSICS
    • G06COMPUTING OR CALCULATING; COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/10Segmentation; Edge detection
    • G06T7/149Segmentation; Edge detection involving deformable models, e.g. active contour models
    • GPHYSICS
    • G06COMPUTING OR CALCULATING; COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/10Segmentation; Edge detection
    • G06T7/162Segmentation; Edge detection involving graph-based methods
    • GPHYSICS
    • G06COMPUTING OR CALCULATING; COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/60Analysis of geometric attributes
    • G06T7/68Analysis of geometric attributes of symmetry
    • GPHYSICS
    • G06COMPUTING OR CALCULATING; COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/70Determining position or orientation of objects or cameras
    • G06T7/73Determining position or orientation of objects or cameras using feature-based methods
    • GPHYSICS
    • G06COMPUTING OR CALCULATING; COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2200/00Indexing scheme for image data processing or generation, in general
    • G06T2200/04Indexing scheme for image data processing or generation, in general involving 3D image data
    • GPHYSICS
    • G06COMPUTING OR CALCULATING; COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10072Tomographic images
    • G06T2207/10081Computed x-ray tomography [CT]
    • GPHYSICS
    • G06COMPUTING OR CALCULATING; COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/20Special algorithmic details
    • G06T2207/20021Dividing image into blocks, subimages or windows
    • GPHYSICS
    • G06COMPUTING OR CALCULATING; COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/20Special algorithmic details
    • G06T2207/20072Graph-based image processing
    • GPHYSICS
    • G06COMPUTING OR CALCULATING; COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/20Special algorithmic details
    • G06T2207/20112Image segmentation details
    • GPHYSICS
    • G06COMPUTING OR CALCULATING; COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/20Special algorithmic details
    • G06T2207/20112Image segmentation details
    • G06T2207/20116Active contour; Active surface; Snakes
    • GPHYSICS
    • G06COMPUTING OR CALCULATING; COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/30Subject of image; Context of image processing
    • G06T2207/30004Biomedical image processing
    • G06T2207/30101Blood vessel; Artery; Vein; Vascular
    • GPHYSICS
    • G06COMPUTING OR CALCULATING; COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/30Subject of image; Context of image processing
    • G06T2207/30172Centreline of tubular or elongated structure

Definitions

  • the present invention relates to an image processing apparatus, an image processing program, and an operation method of the image processing apparatus that divide an image into a plurality of regions using a discrete optimization method such as graph cut or QPBO (Quadratic Pseudo-Boolean Optimization).
  • a discrete optimization method such as graph cut or QPBO (Quadratic Pseudo-Boolean Optimization).
  • Non-Patent Document 1 and Non-Patent Document 2 show a method for obtaining a solution by converting high-order energy into secondary energy.
  • Patent Document 1 specifies a specific energy by setting energy based on the contour. A method of segmenting an object region having an outline and other regions is shown.
  • Non-Patent Document 2 as one application example in which high-order energy is used for segmentation, high-order energy is set for pixels in a rectangular window. This method has a problem that it is difficult to sufficiently sparse high-order energy and the amount of calculation is large.
  • an object of the present invention is to provide an image processing apparatus, an image processing program, and an operation method of the image processing apparatus that use high-order energy effective for segmentation.
  • the image processing apparatus of the present invention uses an energy minimization method to set each pixel constituting image data as a first variable representing a target region or a second variable representing a non-target region excluding the target region.
  • An image processing apparatus for labeling A first contour existing on the object region side having a shape similar to the contour of the object region, and a non-object region side existing at a position facing the first contour across the contour of the object region Contour specifying means for specifying the existing second contour, i ( ⁇ 1) pixels forming the entire first contour, and (Ni) (N ⁇ 4, forming the entire second contour) (N ⁇ i ⁇ 1) pixel selection means for selecting pixels, the i pixels of the first contour all belong to the object area, and all (Ni) pixels of the second contour are all Energy setting means for setting the Nth-order energy when satisfying the condition of belonging to the non-target region to be smaller than the Nth-order energy when not satisfying the condition, and labeling by minimizing the Nth-order energy And labeling means for attaching.
  • An image processing program uses a computer to store a first variable representing a target area for each pixel constituting image data or a second non-target area excluding the target area using an energy minimization method.
  • a program for functioning as an image processing device for labeling a variable of A first contour existing on the object region side having a shape similar to the contour of the object region, and a non-object region side existing at a position facing the first contour across the contour of the object region Contour specifying means for specifying the existing second contour, i ( ⁇ 1) pixels forming the entire first contour, and (Ni) (N ⁇ 4, forming the entire second contour) (Ni ⁇ 1) pixel selection means for selecting pixels, i pixels of the first contour all belong to the object area, and (Ni) pixels of the second contour are all non-pixels
  • Energy setting means for setting the Nth order energy when satisfying the condition of belonging to the target region to be smaller than the Nth order energy when not satisfying the condition, and labeling by minimizing the Nth order energy It is also intended to function as
  • the operation method of the present invention includes an outline specifying unit, a pixel selecting unit, an energy setting unit, and a labeling unit, and represents each object constituting the image data by using an energy minimization method.
  • An operation method of an image processing apparatus for labeling a first variable or a second variable representing a non-target region excluding the target region The contour specifying means has a shape similar to the contour of the object region, the first contour existing on the object region side, and the non-existing position existing at a position facing the first contour across the contour of the object region
  • a contour specifying step for specifying the second contour existing on the object region side, i ( ⁇ 1) pixels forming the entire first contour by the pixel selection means, and the entire second contour are formed ( Ni) (N ⁇ 4, Ni ⁇ 1) pixel selection step, and the energy setting means all the i pixels of the first contour belong to the object region, and Energy that is set so that the Nth order energy when the condition that (Ni) pixels of two contours all belong to the non-target region is smaller than
  • the “shape similar to the contour of the object region” means a shape having a shape similar to the contour shape of the object region, and the first contour is separated from the contour of the object region by a substantially constant distance. What is located inside the region is preferable, and the second contour is preferably separated from the contour of the target region by a substantially constant distance and positioned outside the target region (non-target region side).
  • the contour specifying means specifies the first contour and the second contour so that the first contour and the second contour are separated by a predetermined distance.
  • the image processing apparatus further includes contour estimation means for estimating the contour of the object region using image data, and the contour specifying means specifies the first contour and the second contour based on the estimated contour. May be.
  • the contour estimation means may estimate the contour of the object region from the image data using differential filtering capable of detecting an edge.
  • the contour estimation means may estimate the contour of the target area from the image data using a shape model expressed by a predetermined parameter.
  • the energy setting means sets all the Nth-order energies when the conditions are not satisfied to an equal value.
  • the labeling means may use a QPBO algorithm for energy minimization.
  • Block diagram showing the configuration of the first image processing apparatus Flow chart showing the flow of image processing Diagram for explaining a method to generate a three-dimensional CPR image from a contrast CT image of the heart
  • Example of graph when classifying into 2 classes An example of a graph representing a quadratic expression in which higher-order energy is converted by a pseudo-Boolean function Diagram for explaining the cross-sectional shape of a blood vessel Example of shape where blood vessel region appears on image
  • the figure for demonstrating the pattern in which the blood vessel region is separated from the background region Diagram for explaining the setting of secondary energy The figure which shows the relationship between the outline of a target object area
  • the figure for demonstrating the method of selecting the pixel on a 1st outline, and the pixel on a 2nd outline A diagram for explaining how to set a circle and Nth order energy that are likely to
  • the image processing apparatus of the present invention is realized by loading an image processing program into a computer and executing it.
  • the image processing program is stored and distributed in a storage medium such as a CD-ROM, and is installed in the computer from the storage medium such as a CD-ROM.
  • the program is distributed via a network such as the Internet and installed in the computer.
  • FIG. 1 is a block diagram showing a configuration of a first image processing apparatus 1 according to the present invention
  • FIG. 2 is a flowchart showing a processing flow of the image processing apparatus 1 according to the present invention.
  • the image processing apparatus 1 includes an image data input receiving unit 11, a preprocessing unit 12, a contour specifying unit 13, a pixel selecting unit 14, an energy setting unit 15, a labeling unit 16, Display means 17.
  • the image data input receiving means 11 receives the input of the image data P to be processed and stores it in the storage device (S1).
  • the image data P is a simple X-ray image, a tomographic image captured by a CT apparatus or an MRI apparatus, and the like. In addition, a plurality of organs, tissues, or lesions are photographed in the image data P.
  • the processing for detecting the contour of the coronary artery from the contrast CT image of the heart where the image data P is the contrast CT image (voxel data) of the heart, will be described.
  • the preprocessing means 12 extracts a coronary artery centerline l (broken line) from the contrast CT image by a known coronary artery extraction process (S2).
  • a coronary artery extraction process S2
  • Various methods have been proposed for coronary artery extraction processing, including the methods disclosed in Japanese Patent No. 4709290 and Japanese Patent No. 4717935.
  • Cross-sectional images S 1 , S 2 , S 3 perpendicular to the blood vessel center line are generated from the contrast CT image (left figure in FIG. 3), and the cross-section is such that the central coordinates of the blood vessel in this cross-sectional image are located on the straight line L
  • the images S 1 , S 2 , S 3 are arranged (the right diagram in FIG. 3), and a three-dimensional CPR image is generated (S3).
  • a contrast CT image of the heart (A in the figure, processing is performed on a white blood vessel) to a three-dimensional CPR image (B and C in the figure, where C is a cross-sectional view of the blood vessel)
  • a in the figure processing is performed on a white blood vessel
  • B and C in the figure, where C is a cross-sectional view of the blood vessel An example of generating is shown.
  • the three-dimensional CPR image generated by the preprocessing means 12 is segmented into two areas, a background area and a blood vessel area.
  • Image segmentation is represented by the following expression (1) when defined as variables ⁇ x 1 , x 2 ,..., X n ⁇ and x ⁇ ⁇ 1, 0 ⁇ corresponding to each pixel of the image. It can be solved by minimizing the energy function consisting of the next energy E i and the second energy E ij .
  • V is a set of pixels constituting the image
  • Ni is a set of pixels adjacent to the pixel i.
  • the primary energy E i (xi) in the equation (1) is a value that depends only on the label assigned to each pixel, and has a direct influence that depends on which label is assigned to each pixel.
  • the sum of the second items defines energy so as to reflect prior knowledge of how the labels given to adjacent pixels should be related.
  • a directed graph representing Expression (1) is set.
  • a primary energy E i (x i ) is defined for the edge from the vertex s to the vertex x i and the edge from the vertex x i to the vertex t, and the edge connecting the vertices corresponding to the adjacent pixels is defined.
  • a secondary energy E ij (x i , x j ) is defined.
  • the energy (weight) of each edge is determined so that the label that gives the minimum energy corresponds to the minimum cut. Good.
  • the third-order or higher energy function is converted into a second-order energy function, and the converted second-order energy function is optimized by a conventionally known QPBO algorithm.
  • the higher order energy function of the original third order or higher is optimized by a conventionally known QPBO algorithm.
  • Rother et al. Disclose a method of converting a graph by adding auxiliary variables z 0 and z 1 .
  • auxiliary variables z 0 and z 1 For example, assume that a cubic energy function is given by the following equation.
  • the object is a blood vessel and the non-object is a background.
  • the cross-sectional shape of the blood vessel is substantially circular.
  • a circle (broken line) is applied to the blood vessel contour on the image, the inside of the circle contour is the blood vessel region, and the outside is the background region. is there.
  • a graph cut is performed by setting a small energy only for a combination of two pixels in which the pixel located in the background region takes 0 and the pixel located in the blood vessel region takes 1 Can be separated into an object area and a background area.
  • FIG. 9A when the left or right circle is separated from the background (black part) using the inside (white part) of the left or right circle as a blood vessel, and a part of the circle is missing as shown in FIG. 9B.
  • FIG. 9C There is a case where either one of the shapes (white portion) is separated as a blood vessel, or a shape shown in FIG. 9C.
  • the secondary energy is set to a small energy when two pixels straddle the outline of a circle (see FIG.
  • both patterns A and B in FIG. 9 the sum of energy decreases in proportion to the length of the contour of the circle, so that the separation of both patterns A and B in FIG. 9 having the same contour length is the same. May happen with probability. Further, even in the pattern as shown in FIG. 9C, the sum of the energy is relatively small, so it is one of the solution candidates.
  • the blood vessels have a shape close to a circle, it is not preferable that the blood vessels are separated in a pattern such as B and C in FIG.
  • N-order energy using the pixel values of N ( ⁇ 4) pixels as variables is used.
  • N the number of pixels located in the background around a circular outline having a size close to that of a blood vessel take 0 ( ⁇ in FIG. 7), and pixels located in a blood vessel around the outline become 1 (see FIG. 7).
  • a small energy is set only for the combination that takes 7). Since the circles overlap each other, the N-order energy of the two circles does not decrease when one decreases. Since one circle is selected exclusively in this way, the image is segmented so that the blood vessel region has a shape close to a circle.
  • An energy function E obtained by adding higher-order energy of 4th order or higher is defined as the following equation.
  • c is a set of pixels constituting higher-order energy
  • Xc is a vector composed of binary variables (0 or 1) of related pixels.
  • C represents a set of c. Note that the order of each c may be different, for example, the number of pixels associated with a large circle is large and small for a small circle.
  • the contour specifying means 13 specifies the inner contour and the outer contour of the blood vessel (object) region. Specifically, as shown in FIG. 11, the first contour (broken line) existing on the blood vessel (object) region side across the blood vessel contour (solid line) and the first contour across the blood vessel contour The second contour (one-dot chain line) existing on the background (non-object) region side existing at the specified position is specified.
  • the first contour and the second contour are shaped similar to the blood vessel contour. Since the contour of the blood vessel is substantially circular, a circle close to the size of the blood vessel is set, a circle with the first contour is provided inside, and a circle with the second contour is provided outside.
  • a gap is provided at a predetermined distance between the first contour on the blood vessel region side (inner circle) and the second contour on the background region side (outer circle) (S4).
  • a plurality of arcs may be set by dividing the circle of the first contour and the second contour. Dividing in this way has the effect of reducing energy even if only some of the patterns of the circles match, so a more flexible segmentation result can be obtained.
  • the pixel selection means 14 searches for a place where a blood vessel-like shape appears from the three-dimensional CPR image, and as shown in FIG. 13, i pixels on the first contour and N ⁇ on the second contour. Select i pixels.
  • the energy setting means 15 sets high-order energy together with primary energy and secondary energy in the graph.
  • the secondary energy is set at the edge between the vertices x i and x j corresponding to the adjacent pixels, and the smaller energy is set as the difference between the adjacent pixel values is larger (S6).
  • the energy setting means 15 specifies a position where a circular shape appears in the three-dimensional CPR image and sets the higher order energy in order to efficiently set the higher order energy. If all the i pixels in the first contour have pixel values that are likely to be blood vessels, and all the (Ni) pixels in the second contour have pixel values that seem to be the background, the condition as the contour is satisfied. If the condition is satisfied, the high-order energy is set to a small energy (its weight is a negative coefficient and the absolute value is large). If this condition is not satisfied, the Nth-order energy satisfies this condition Set to be larger than the energy of.
  • the three-dimensional OOF filter has a spherical filter shape with a radius R (circular in a two-dimensional cross section), and a response appears where a circle with a radius R (dashed line) appears as shown in FIG. Therefore, a three-dimensional OOF filter is applied to each position of the three-dimensional CPR image P to identify a position where a circular shape appears from the three-dimensional CPR image. Further, the three-dimensional OOF filter can calculate the second-order partial differential on the surface of the sphere, and can calculate the likelihood of a cylindrical structure from three eigenvalues obtained by eigenvalue decomposition of the second-order partial differential matrix.
  • FIGS. 15A and 15B show a state in which a circle is detected at a plurality of positions by applying a three-dimensional OOF filter to each position of the three-dimensional CPR image.
  • the pixel selecting means 14 selects N pixels on the plane Q, and the energy setting means 15 sets the Nth order energy (S8).
  • the higher order is set so that the value is relatively small. Set energy.
  • the N pixels selected on the plane Q at the position where the response of the three-dimensional OOF filter is low satisfy the contour condition
  • the N pixels selected on the plane Q at the position where the response is high meet the contour condition.
  • the Nth-order energy is set so as to be larger than the case where it is satisfied, and energy smaller than the higher-order energy when the condition of the circle is not satisfied is set. For example, at a position where the response of the 3D OOF filter is high, a large negative energy (coefficient is negative and the absolute value is large) is set as the Nth order energy when the contour condition is satisfied, and the response of the 3D OOF filter Is set to a small negative energy (coefficient is negative and absolute value is small) as the Nth order energy when the contour condition is satisfied.
  • the reason why the energy magnitude is adaptively changed as described above is because it is known in advance that the higher the response of the filter, the higher the probability that a circle will exist.
  • the energy weight may be set according to the probability. This probability can be calculated as a ratio at which the combination of labels in which the Nth-order energy that is set by the above-described processing decreases matches the label of correct data.
  • the energy setting means 15 may set all the Nth-order energies when the contour condition is not satisfied to an equal value. In this way, if the energy when the contour condition is not satisfied is set to an equal value, the energy becomes sparse, so that the higher order can be converted to the lower order with a few additional variables, and the solution can be obtained relatively quickly.
  • the labeling means 16 converts the set higher-order (third-order or higher) energy into a second-order energy function, and then optimizes the second-order energy function after conversion with a conventionally known QPBO algorithm. Then, segmentation into two regions, a blood vessel region and a background region (S9).
  • the blood vessel region is segmented so that the extracted circles (A in FIG. 15) are three-dimensionally continuous as shown in FIGS. FIG. 15B). Further, segmentation with a shape deviating from a circle can be prevented.
  • the display means 17 displays an image of only the blood vessel region on the display device according to the result of the labeling means 16 (S10).
  • higher-order energy is set in a circular shape such as a blood vessel
  • shape of the contour set as higher-order energy may be any shape. Pixels on the first contour on the object region side and the second contour on the non-object region side are selected across various contours, and all the pixels on the first contour belong to the object region, and the second contour
  • the higher-order energy when all the pixels belong to the non-target region is set to be smaller than the higher-order energy when this condition is not satisfied.
  • the first contour is always away from the contour of the target region by a substantially constant distance inside (target region side) so as to be similar to the contour of the target region
  • the second contour is It is preferable that the distance from the contour of the object region is always a constant distance outside (non-target region side) so as to be similar to the contour, but the distance from the contour of the target region is within a predetermined error range. May be specified to be different.
  • FIGS. 17A and 17B four or more pixels are used to form the contour, and if there is at least one pixel in at least the first contour or the second contour, the shape of the contour is specified. can do.
  • FIG. 17A shows an example in which the object area and the non-object area are separated along a straight line
  • FIG. 17B shows an example in which the object area and the non-object area are separated in a corner shape.
  • various shapes such as a triangle and a rectangle can be specified.
  • the shape of the surface or curved surface may be specified, and the surface or curved surface in the three-dimensional image may be detected and separated.
  • the contour estimation means 18 can use differential filtering capable of detecting edges as a method for estimating the contour.
  • the response of the differential filter changes from positive to negative in the vicinity of the contour. Therefore, the position having the positive response of the differential filter is the first contour (object area) and the position having the negative response is the second.
  • the pixel is selected as a contour (non-object region), the pixel selecting unit 14 selects pixels on the first contour and the second contour, and the energy setting unit 15 has a first contour at a position having a positive response and a negative response.
  • the energy is set so as to be minimized with respect to the combination in the second contour at the position having.
  • an appropriate image is generated such as an image of only the target region and the image is displayed on the display unit 17. Display on the device.
  • the contour estimation means 18 can also estimate the contour of the target region from the image data using a shape model expressed by a predetermined parameter.
  • a shape model expressed by a predetermined parameter.
  • ASM Active Shape Model
  • the shape of the object area is subjected to principal component analysis, independent component analysis, etc., and the average shape and a vector for transforming the average shape are obtained, and the parameters of each vector are determined.
  • high-order energy may be set at the corresponding contour position while changing various parameters. At this time, if the contour of the first contour on the non-target region (background) and the second contour on the target region are selected at a certain interval so that the contour may change flexibly.
  • the shape model expressed by a predetermined parameter may be divided into a plurality of line segments, and the energy may be set for each line segment as if the circle was divided into a plurality of arcs.
  • Non-Patent Document 2 it is possible to improve the accuracy of area division of a structure having a specific shape by using high-order energy.
  • the present invention is characterized in that pixels are flexibly selected without being fixed like a rectangular window.
  • the higher order energy to be set can be made sparse, and a solution can be obtained at high speed with a small amount of calculation.
  • the method for selecting pixels is preferably performed along the contour shape to be segmented. For example, since the blood vessel contour in the three-dimensional image has a cylindrical shape, a combination of pixels that are circular is selected.
  • the target contour shape is learned as a statistical model (for example, ASM), pixels along the contour of the model are selected.
  • ASM statistical model
  • energy minimization has been described using a graph cut technique, but any method that minimizes a discrete function may be used.
  • Specific methods for energy minimization include graph cut, QPBO, Tree Reweighted Message Passing (for details, see V. Kolmogorov, Convergent Tree-reweighted Message Passing for Energy Minimization, IEEE PAMI, 28 (10), October 2006 See).
  • Also, for how to label each object area and non-object area using graph cut see 'Y. Boykov and V. Kolmogorov, An Experimental Comparison of Min-Cut / Max-Flow Algorithms for Energy Minimization in Vision, IEEE IEEE PAMI, 26 (9), pp. 1124-1137, 2004.
  • a method of labeling each of the object area and the non-object area using QPBO is described in detail in Non-Patent Document 2.

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Health & Medical Sciences (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Medical Informatics (AREA)
  • General Health & Medical Sciences (AREA)
  • Biophysics (AREA)
  • Optics & Photonics (AREA)
  • Pathology (AREA)
  • Radiology & Medical Imaging (AREA)
  • Biomedical Technology (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Molecular Biology (AREA)
  • Surgery (AREA)
  • Animal Behavior & Ethology (AREA)
  • High Energy & Nuclear Physics (AREA)
  • Public Health (AREA)
  • Veterinary Medicine (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Geometry (AREA)
  • Cardiology (AREA)
  • Oral & Maxillofacial Surgery (AREA)
  • Dentistry (AREA)
  • Software Systems (AREA)
  • Pulmonology (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Evolutionary Computation (AREA)
  • Evolutionary Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Artificial Intelligence (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Apparatus For Radiation Diagnosis (AREA)
  • Image Analysis (AREA)

Abstract

【解決手段】2値のラベル付けをする際に、対象物領域の輪郭に相似した形状を有する、輪郭特定手段(13)で対象物領域側に存在する第1輪郭、および、非対象物領域側に存在する第2輪郭を特定し、画素選択手段(14)で第1輪郭と第2輪郭の全体を形成する画素をN個選択して、エネルギー設定手段(15)で第1輪郭の画素が全て前記対象物領域に属し、かつ、第2輪郭の画素が全て非対象領域に属するという条件を満たす場合のN次のエネルギーに小さいエネルギーを設定した上で、エネルギーを最小化してラベル付けをする。

Description

画像処理装置、画像処理プログラムおよび画像処理装置の作動方法
 本発明は、グラフカットやQPBO(Quadratic Pseudo-Boolean Optimization)などの離散最適化手法を用いて複数の領域に画像を分割する画像処理装置、画像処理プログラムおよび画像処理装置の作動方法に関する。
 近年、グラフの最小切断(グラフカット)アルゴリズムを使ったエネルギー最小化が、画像処理に盛んに応用されるようになってきた。特に画像のセグメンテーション(領域分割)の問題をエネルギー最小化問題として効率的に解く方法が提案されている。
 グラフカットは、2変数に依存する2次のエネルギーを基本として発展してきたが、Kolmogorovらは3変数に依存する劣モジュラな3次のエネルギーを最小化する方法を提案している(例えば、非特許文献1)。
 さらに、4次以上の高次のエネルギー関数を2次のエネルギー関数に変換して解く方法が提案されている。変換後のエネルギー関数の2次項の係数が全て負の場合をサブモジュラと呼び、Max-Flow/Min-Cutアルゴリズムにより高速に最小解を得ることができる。一方、変換後のエネルギー関数の2次項の係数に正の値が含まれる場合(非サブモジュラ)は、QPBOアルゴリズムを用いて最小化または近似解を得ることで、元の高次エネルギー関数の最適化が可能である。このエネルギー関数の変換は複数知られており、例えばRotherらによりこの変換が提案されている(例えば、非特許文献2)。
 また、グラフカットの手法を用いて、抽出対象となる対象物の領域を抽出するために、対象物の輪郭を表す閉曲線との距離に基づいてエネルギーを設定し、オブジェクト領域とそれ以外の領域とにセグメンテーションする手法を提案したものがある(例えば、特許文献1)。
V. Kolmogorov and R. Zabih, What energy functions can be minimized via Graph Cuts?, IEEE PAMI, 26(2), 2004 C. Rother, P. Kohli, W. Feng, J. Jia, Minimizing Sparse Higher Order Energy Functions of Discrete Variables, CVPR, 2009
特開2012-027713号公報
 非特許文献1および非特許文献2には、高次エネルギーを2次のエネルギーに変換して解を得る方法が示され、特許文献1には、輪郭に基づいてエネルギーを設定することで特定の輪郭を持つ対象物領域とそれ以外の領域とにセグメンテーションする手法が示されている。非特許文献2には、高次のエネルギーを利用してセグメンテーションする一つの応用例として、矩形のウィンドウ内の画素に対して高次のエネルギーを設定している。この方法では高次エネルギーを十分に疎にすることが難しく、計算量が大きいという問題があった。
 そこで、本発明では、セグメンテーションに有効な高次のエネルギーを利用した画像処理装置、画像処理プログラムおよび画像処理装置の作動方法を提供することを目的とする。
 本発明の画像処理装置は、エネルギー最小化手法を用いて、画像データを構成する各画素を対象物領域を表す第1の変数または前記対象物領域を除く非対象領域を表す第2の変数にラベル付けをする画像処理装置であって、
 対象物領域の輪郭に相似した形状を有する、対象物領域側に存在する第1輪郭、および、前記対象物領域の輪郭を挟んで第1輪郭と対向した位置に存在する非対象物領域側に存在する第2輪郭、を特定する輪郭特定手段と、第1輪郭の全体を形成するi(≧1)個の画素と、第2輪郭の全体を形成する(N-i)(N≧4、N-i≧1)個の画素を選択する画素選択手段と、第1輪郭の前記i個の画素が全て対象物領域に属し、かつ、第2輪郭の(N-i)個の画素が全て非対象領域に属するという条件を満たす場合のN次のエネルギーを、前記条件を満たさない場合のN次のエネルギーよりも小さくなるように設定するエネルギー設定手段と、N次のエネルギーを最小化してラベル付けをするラベリング手段とを備えることを特徴とするものである。
 本発明の画像処理プログラムは、コンピュータを、エネルギー最小化手法を用いて、画像データを構成する各画素を対象物領域を表す第1の変数または前記対象物領域を除く非対象領域を表す第2の変数にラベル付けをする画像処理装置として機能させるためのプログラムであって、
 対象物領域の輪郭に相似した形状を有する、対象物領域側に存在する第1輪郭、および、前記対象物領域の輪郭を挟んで第1輪郭と対向した位置に存在する非対象物領域側に存在する第2輪郭、を特定する輪郭特定手段と、第1輪郭の全体を形成するi(≧1)個の画素と、第2輪郭の全体を形成する(N-i)(N≧4、N-i≧1)個の画素を選択する画素選択手段と、第1輪郭のi個の画素が全て対象物領域に属し、かつ、第2輪郭の(N-i)個の画素が全て非対象領域に属するという条件を満たす場合のN次のエネルギーを、条件を満たさない場合のN次のエネルギーよりも小さくなるように設定するエネルギー設定手段と、N次のエネルギーを最小化してラベル付けをするラベリング手段として機能させるものであることを特徴とするものである。
 本発明の作動方法は、輪郭特定手段と、画素選択手段と、エネルギー設定手段と、ラベリング手段とを有し、エネルギー最小化手法を用いて、画像データを構成する各画素を対象物領域を表す第1の変数または前記対象物領域を除く非対象領域を表す第2の変数にラベル付けをする画像処理装置の作動方法であって、
 輪郭特定手段により、対象物領域の輪郭に相似した形状を有する、対象物領域側に存在する第1輪郭、および、対象物領域の輪郭を挟んで第1輪郭と対向した位置に存在する前記非対象物領域側に存在する第2輪郭、を特定する輪郭特定ステップと、画素選択手段により第1輪郭の全体を形成するi(≧1)個の画素と、第2輪郭の全体を形成する(N-i)(N≧4、N-i≧1)個の画素を選択する画素選択ステップと、エネルギー設定手段により第1輪郭の前記i個の画素が全て対象物領域に属し、かつ、第2輪郭の(N-i)個の画素が全て非対象領域に属するという条件を満たす場合のN次のエネルギーを、前記条件を満たさない場合のN次のエネルギーよりも小さくなるように設定するエネルギー設定ステップと、ラベリング手段により、N次のエネルギーを最小化してラベル付けをするラベリングステップを備えたことを特徴とするものである。
 「対象物領域の輪郭に相似した形状」とは、対象物領域の輪郭の形状に似通った形状を有するものをいい、第1輪郭は対象物領域の輪郭から略一定の距離分離れて対象物領域内側に位置するものが好ましく、第2輪郭は対象物領域の輪郭から略一定の距離分離れて対象物領域外側(非対象領域側)に位置するものが好ましい。
 また、輪郭特定手段が、第1輪郭と第2輪郭とが所定の距離離間するように第1輪郭および第2輪郭を特定するものが好ましい。
 また、画像データを用いて前記対象物領域の輪郭を推定する輪郭推定手段をさらに備え、輪郭特定手段が、推定した輪郭に基づいて、前記第1輪郭および前記第2輪郭を特定するものであってもよい。
 また、輪郭推定手段が、エッジを検出可能な微分フィルタリングを用いて画像データから対象物領域の輪郭を推定するものであってもよい。
 また、輪郭推定手段が、所定のパラメータで表現された形状モデルを用いて前記画像データから前記対象領域の輪郭を推定するものであってもよい。
 さらに、エネルギー設定手段が、条件を満たさない場合のN次のエネルギーを全て等しい値に設定するものが望ましい。
 さらにまた、ラベリング手段が、エネルギーの最小化のためにQPBOアルゴリズムを用いるものであってもよい。
第1の画像処理装置の構成を示したブロック図 画像処理の流れを示すフローチャート 心臓の造影CT画像から三次元のCPR画像を生成する手法を説明するための図 心臓の造影CT画像から三次元のCPR画像を生成した例 2クラスに分類する場合のグラフの一例 高次のエネルギーを擬ブール関数で変換した2次式を表すグラフの一例 血管の断面形状を説明するための図 血管領域が画像上にあらわれる形状の例 血管領域が背景領域から分離されるパターンを説明するための図 2次のエネルギーの設定を説明するための図 対象物領域の輪郭と第1輪郭と第2輪郭の関係を示す図 第1輪郭と第2輪郭を離間させた一例 第1輪郭上の画素と第2輪郭上の画素を選択する手法を説明するための図 三次元OOFフィルタを用いて検出した血管らしい円とN次のエネルギーの設定方法を説明するための図 三次元CPR画像中から検出された複数の円から血管領域をセグメンテーションする手法を説明するための図 従来の手法による血管領域の抽出と本発明の手法による血管領域の抽出の違いを説明するための図 N次のエネルギーを設定する直線の形状とその形状を構成する画素のとり方を説明するための図 N次のエネルギーを設定する角の形状とその形状を構成する画素のとり方を説明するための図 第2の画像処理装置の構成を示したブロック図 対象領域の輪郭の円を複数の円弧で定義する手法を説明するための図
 以下、図面を参照して本発明の実施形態を詳細に説明する。本発明の実施形態では、グラフカットの手法を用いて画像中の各画素を対象物領域を表す第1の変数または対象物領域を除く非対象領域を表す第2の変数に領域分割する手法について説明する。本発明の画像処理装置は、画像処理プログラムがコンピュータにロードされて実行されることにより実現する。また、画像処理プログラムはCD-ROM等の記憶媒体に記憶されて配布され、CD-ROM等の記憶媒体からコンピュータにインストールされる。あるいは、インターネット等のネットワークを介してプログラムが配布されて、コンピュータにインストールされる。
 図1に、本発明に係る第1の画像処理装置1の構成を示したブロック図を示し、図2に、本発明に係る画像処理装置1の処理の流れを示したフローチャートを示す。
 図1に示すように、画像処理装置1は、画像データ入力受付手段11と、前処理手段12と、輪郭特定手段13と、画素選択手段14と、エネルギー設定手段15と、ラベリング手段16と、表示手段17とを備えている。
 画像データ入力受付手段11は、画像処理対象の画像データPの入力を受け付け、記憶装置に記憶する(S1)。画像データPは、単純X線撮影画像や、CT装置やMRI装置で撮影された断層画像等である。また、画像データPには、複数の臓器や組織あるいは病変部が撮影されている。
 以下、画像データPが心臓の造影CT画像(ボクセルデータ)であり、この心臓の造影CT画像から冠動脈の輪郭を検出する処理について説明する。
 前処理手段12は、既知の冠動脈抽出処理により、図3に示すように、造影CT画像から冠動脈の中心線l(破線)を抽出する(S2)。冠動脈抽出処理は、特許第4709290号や特許第4717935号に示される方法をはじめ、種々の方法が提案されている。造影CT画像から血管の中心線に直行する断面画像S,S,Sを生成し(図3の左図)、この断面画像の血管の中心座標が直線L上に位置するように断面画像S,S,Sを配置して(図3の右図)、三次元のCPR画像を生成する(S3)。図4のA~Cに心臓の造影CT画像(同図のA、白い血管に対して処理を行う)から三次元のCPR画像(同図のB、C、ただしCは血管の断面の図)を生成した例を示す。
 前処理手段12で生成した三次元のCPR画像を背景領域、血管領域の2つに領域にセグメンテーションする。
 まず、画像を2つに領域に分割するために用いるグラフカットの手法について説明する。画像のセグメンテーションは、画像の各画素に対応する変数{x,x,・・・,x}、x∈{1,0}と定義したとき、次式(1)で表される1次のエネルギーEと2次のエネルギーEijからなるエネルギー関数を最小化することで解くことができる。
Figure JPOXMLDOC01-appb-M000001
 Vは画像を構成する画素の集合であり、Niは画素iと隣接関係がある画素の集合である。
 式(1)の1次のエネルギーEi(xi)は、各画素に割り当てられるラベルのみに依存する値であり、各画素にどのラベルを割り当てたかによって決まる直接的な影響が現れる。第2項目の和は、隣接する画素に与えられるラベルがどのような関係にあるべきかという事前知識を反映するようにエネルギーを定義したものである。
 図5に示すように、式(1)を表した有向グラフを設定する。グラフは各画素に対応する頂点xi(または、頂点x)と、さらに各ラベルに対応する2つの頂点s(=0)と頂点t(=1)を持ち、頂点xi、頂点xそれぞれについて、頂点sから頂点xiへ向かうエッジと頂点xiから頂点tへ向かうエッジと、隣接する画素に対応する頂点x、x間のエッジで構成される。さらに、頂点sから頂点xiへ向かうエッジと頂点xiから頂点tへ向かうエッジに対して、1次のエネルギーE(x)を定義し、隣接する画素に対応する頂点を結ぶエッジに対しては、2次のエネルギーEij(x,x)を定義する。式(1)であらわされるエネルギーE(x)の最小化により各画素に対するラベルを得るためには、最小エネルギーを与えるラベルと最小切断が対応するように各エッジのエネルギー(重み)を定めればよい。
 例えば、この2次のエネルギーEij(x,x)が下式で与えられるとする。
Figure JPOXMLDOC01-appb-M000002
 このとき2次のエネルギーEijを擬ブール式で書くと、下のような2次式となる。
Figure JPOXMLDOC01-appb-M000003
 擬ブール式において、2次の項の係数が全て負の場合をサブモジュラと呼び、Max-Flow/Min-Cutアルゴリズムにより高速に最小解を得ることができる。一方、2次の項の係数に正の値が含まれる場合(つまり、非サブモジュラ)には、QPBOアルゴリズムを用いて最小化または準最適解を得ることができる(詳細は、C. Rother, V. Kolmogorov, V. Lempitsky, M. Szummer, Optimizing Binary MRFs via Extended Roof Duality, CVPR, 2007を参照)。
 さらに、3次以上のエネルギー関数については、3次以上のエネルギー関数を2次のエネルギー関数に変換した上で、変換後の2次のエネルギー関数を従来から知られているQPBOアルゴリズムで最適化することにより、元の3次以上の高次エネルギー関数の最適化が可能である。
 例えば、Rotherらは補助変数z,zを追加してグラフを変換する手法を開示している。例えば、3次のエネルギー関数が、下式で与えられているとする。
Figure JPOXMLDOC01-appb-M000004
 この3次のエネルギーEijkを擬ブール式で書くと、下のような2次式となる。
Figure JPOXMLDOC01-appb-M000005
 この変換をグラフで表すと、図6のようになる。このように変換した後、補助変数を含めて最小化することで解を得る。その結果得られる解は、もとの高次関数の解と等しい(詳細は、C. Rother, P. Kohli, W. Feng, J. Jia, Minimizing Sparse Higher Order Energy Functions of Discrete Variables, CVPR, 2009を参照)。
 4次以上のエネルギーについても同様の変換を行うことで解を得ることが可能である。この変換は、高次エネルギーの組合せ表のうち値が異なる要素ごとに、変換を行う必要がある。複雑な場合には追加する補助変数の数が多くなり最適化が難しくなるが、エネルギーが疎の場合、例えば、式(4)のように特定の値のみが小さい場合は、相対的にエネルギーが小さい要素だけを変換すればよい。従って変換後の問題が比較的単純で、実用的な計算時間で解が得られる。
 以下、4次以上のエネルギーを用いて、形状が既知である対象物が存在する対象物領域とそれ以外の領域(非対象物領域)にセグメンテーションする手法について説明する。以下、対象物を血管とし、非対象物を背景とする。
 図7に示すように、血管の断面形状は概ね円形であり、画像上の血管の輪郭に円(破線)をあてはめたとき、その円の輪郭の内側が血管領域であり、外側が背景領域である。
 2次のエネルギーで表現する場合は、2つの画素のうち背景領域に位置する画素が0をとり、血管領域に位置する画素が1を取る組合せのみ、小さなエネルギーを設定してグラフカットを行うことで対象物領域と背景領域に分離することができる。
 図8に示すように、エネルギーを2つの円の輪郭が互いに重なった状態で設定する場合を考える。このとき、円の輪郭に2次のエネルギーを設定して血管領域と背景領域を分離すると、図9のA~Cに示す6通りのパターンに分離される可能性がある。図9のAのように、左または右の円のどちらかの内側(白い部分)を血管として背景(黒い部分)から分離される場合と、図9のBのように円が一部欠けた形状のどちらかの内側(白い部分)を血管として分離される場合と、図9のCに示す形状で分離される場合がある。2次のエネルギーは、2つの画素が円の輪郭を跨ぐ(図10参照)場合に小さいエネルギーを設定したとする。図9のA、Bのいずれのパターンでも、エネルギーの和は円の輪郭の長さに比例して小さくなるので、輪郭の長さが等しい図9のA、Bのいずれのパターンの分離も同じ確率で起こる可能性がある。さらに、図9のCのようなパターンでもエネルギーの和は相対的に小さいので、解の候補の一つである。しかし、血管は円に近い形状をもつことを考慮すると、図9のB、Cのようなパターンで分離されるのは好ましくない。
 そこで、円に近い形状で血管領域をセグメンテーションするために、1次のエネルギーと2次のエネルギーに加えて、N(≧4)個の画素の画素値を変数とするN次のエネルギーを用いる。N次のエネルギーには、血管に近い大きさを持った円形の輪郭の周囲の背景に位置する画素が0(図7の●)をとり、輪郭の周囲の血管に位置する画素が1(図7の○)をとる組合せのみ小さいエネルギーを設定する。円は互いに重なりあう部分があるので、2つの円のN次エネルギーは一方が小さくなるともう一方が小さくなることはない。このように排他的に一方の円が選ばれるように作用するので、血管領域が円に近い形状になるように画像がセグメンテーションされる。
 4次以上の高次エネルギーを加えたエネルギー関数Eを次式のように定義する。
Figure JPOXMLDOC01-appb-M000006
 ここでcは高次エネルギーを構成する画素の集合であり、Xcは関係する画素のバイナリ変数(0または1)からなるベクトルである。Cはcの集合を表す。なお、各cの次数はそれぞれ異なってもよく、例えば大きな円に関連する画素の数は多く、小さな円に対しては少ない。
 まず、輪郭特定手段13で、血管(対象物)領域の内側の輪郭と外側の輪郭を特定する。具体的には、図11に示すように、血管の輪郭(実線)を挟んだ血管(対象物)領域側に存在する第1輪郭(破線)と、血管の輪郭を挟んで第1輪郭と対向した位置に存在する背景(非対象物)領域側に存在する第2輪郭(一点鎖線)を特定する。第1輪郭と第2輪郭は、血管の輪郭に相似した形状にする。血管の輪郭は略円形であるので、血管の大きさに近い円を設定し、さらに、その内側に第1輪郭の円を設け、外側に第2輪郭の円を設ける。
 しかし、完全に血管の大きさをもった円形に固定してセグメンテーションが行われると、柔軟性がなく都合が悪い。実際の画像にあらわれる血管は、完全な円形ではなく少し歪んだ円になることもあり、実際の画像に即した輪郭線が得られるようにしたい。そこで、図12に示すように、血管領域側の第1輪郭(内側の円)と背景領域側の第2輪郭(外側の円)との間を所定の距離離して隙間を設ける(S4)。こうすることで、その間であれば自由な輪郭(破線)をとっても設定した高次エネルギーの条件は成立し、エネルギーが小さくなる。すなわち高次が成立する範囲で、2次のエネルギーが最小となるように柔軟に輪郭線が選択される。また、完全な円を設定するのではなく、図19に示すように、第1輪郭と第2輪郭の円を分割して複数の円弧を設定するようにしてもよい。このように分割すると、円の一部のパターンが一致するだけでもエネルギーを小さくする効果があるため、より柔軟なセグメンテーション結果を得ることができる。
 次に、画素選択手段14は、三次元CPR画像から血管らしい形状が現れている場所を探して、図13に示すように、第1輪郭上のi個の画素と第2輪郭上のN-i個の画素を選択する。
 エネルギー設定手段15は、グラフに1次のエネルギーと2次のエネルギーとともに高次のエネルギーを設定する。1次のエネルギーは頂点s(=0:背景)から頂点xiへ向かうエッジと頂点xiから頂点t(=1:血管)へ向かうエッジとに設定され、各画素の値が血管らしい値を持っている場合は、頂点xiから頂点tへ向かうエッジより頂点sから頂点xiへ向かエッジに小さいエネルギーを設定し、各画素の値が背景らしい値を持っている場合は、頂点sから頂点xiへ向かうエッジより頂点xiから頂点tへ向かうエッジに小さいエネルギーを設定する(S5)。2次のエネルギーは隣接する画素に対応する頂点x、x間のエッジに設定され、隣接画素値の差分が大きいほど小さいエネルギーを設定する(S6)。
 さらに、エネルギー設定手段15は、高次のエネルギーを効率的に設定するために、三次元CPR画像中の円らしい形状があらわれている位置を特定して高次のエネルギーを設定する。第1輪郭のi個の画素が全て血管らしい画素値を持ち、かつ、第2輪郭の(N-i)個の画素が全て背景らしい画素値を持つ場合は輪郭としての条件を満たすので、この条件を満たす場合は、高次のエネルギーは小さいエネルギー(その重みは係数が負で絶対値が大きい)を設定し、この条件を満たさない場合には、N次のエネルギーはこの条件を満たした場合のエネルギーより大きくなるように設定する。
 円形のN次のエネルギーを設定する位置を特定するために、三次元OOFフィルタを用いる。三次元OOFフィルタは半径Rの球状(2次元断面では円形)のフィルタ形状を持ち、図14のAに示すように、半径Rの円(破線)の形状が現れているところにレスポンスがあらわれる。そこで、三次元CPR画像Pの各位置に三次元OOFフィルタを施して、三次元CPR画像から円形らしい形状があらわれている位置を特定する。また、三次元OOFフィルタは球の表面における2階偏微分を算出し、その2階偏微分行列を固有値分解して得た3つの固有値から、円柱構造物らしさを算出することができる。三次元OOFフィルタを施した位置に円柱構造物が存在すれば、3つの固有値のうち2つが大きく、1つがゼロに近い値をとる。同時に得られる3つ固有ベクトルのうち、最もゼロに近い固有値を持つ固有ベクトルが円柱構造物の主軸t(走行方向)に相当する。この主軸と直交する平面Q上に、所定の大きさを持つ血管の輪郭らしい画像が存在している(三次元OOFフィルタの詳細は、M. Law et. al., “Three Dimensional Curvilinear Structure Detection Using Optimally Oriented Flux”, Proceedings of ECCV, pp. 368-382, 2008を参照)。フィルタリング処理によって、このような円が現れている平面Qを検出する(S7)。また血管の正確な大きさは事前にわからないので、様々な大きさの円に対応したOOFフィルタで検出を行う。図15のA、Bに三次元CPR画像の各位置に三次元OOFフィルタを施して複数の位置で円が検出された様子を示す。
 図14のBに示すように、画素選択手段14はこの平面Q上でN個の画素を選択して、エネルギー設定手段15でN次のエネルギーを設定する(S8)。三次元OOFフィルタで三次元CPR画像Pをスキャンし、そのレスポンスが高い位置の平面Qにおいて、選択したN個の画素が輪郭の条件を満たしたとき相対的に小さい値になるように高次のエネルギーを設定する。一方、三次元OOFフィルタのレスポンスが低い位置の平面Qで選択したN個の画素が輪郭の条件を満たす場合には、レスポンスが高い位置の平面Qで選択したN個の画素が輪郭の条件を満たす場合より大きくなるようにN次のエネルギーを設定し、かつ、円の条件を満たしていないときの高次のエネルギーよりは小さいエネルギーを設定する。例えば、三次元OOFフィルタのレスポンスが高い位置では、輪郭の条件を満足したときのN次のエネルギーに大きい負のエネルギー(係数が負で絶対値が大きい)を設定し、三次元OOFフィルタのレスポンスが低い位置では、輪郭の条件を満足したときのN次のエネルギーに小さい負のエネルギー(係数が負で絶対値が小さい)を設定する。
 以上のように適応的にエネルギーの大きさを変える理由は、フィルタのレスポンスが高いほど円が存在する可能性が高いと事前に知っているからである。一方で、データから統計的に算出することも可能である。例えば、血管がどのようにセグメンテーションされるべきかを示す正解のデータがある場合には、三次元OOFフィルタのレスポンスに対してN次のエネルギーが小さくなるケースの確率を統計的に算出(学習)しておき、その確率に応じてエネルギーの重みを設定するようにしてもよい。この確率は、上述の処理によって設定するN次のエネルギーが小さくなるラベルの組み合わせが、正解のデータのラベルと一致する割合として算出できる。
 また、エネルギー設定手段15が、輪郭の条件を満たさない場合のN次のエネルギーを全て等しい値に設定するようにしてもよい。このように輪郭の条件を満たさない場合のエネルギーを等しい値にすれば、エネルギーが疎になるので少ない付加変数で高次を低次に変換することができ、解が比較的高速に得られる。
 ラベリング手段16は、設定された高次(3次以上)のエネルギーを2次のエネルギー関数に変換した上で、変換後の2次のエネルギー関数を従来から知られているQPBOアルゴリズムで最適化して、血管領域と背景領域の2つの領域にセグメンテーションする(S9)。
 ここで、円形の高次エネルギーを設定した時の効果について説明する。図16のAの2次元の原画像に対して、従来の2次元のエネルギーを設定してグラフカット処理を行った場合には、図16のBに示すように、かなり歪んだ円形で抽出される。一方、図16のCに示すような円形の高次エネルギーを設定した場合には、図16のDに示すような円形に近い形状にセグメンテーションが行われる。
 このように円形の高次エネルギーを設定することにより、図15のA、Bに示すように抽出された円形(図15のA)が、3次元的に連続するよう血管領域がセグメンテーションされる(図15のB)。また円形から逸脱した形状でセグメンテーションされることを防ぐことができる。
 表示手段17は、ラベリング手段16の結果に応じて、血管領域のみの画像を表示装置上に表示する(S10)。
 上述では、血管のように輪郭の形状が円形であるものに高次のエネルギーを設定する場合について説明したが、高次エネルギーとして設定する輪郭の形はどんな形でも良い。様々な輪郭を挟んで対象物領域側の第1輪郭と非対象物領域側の第2輪郭上の画素を選択して、第1輪郭の画素が全て対象物領域に属し、かつ、第2輪郭の画素が全て非対象領域に属するという場合の高次のエネルギーを、この条件を満たさない場合の高次のエネルギーよりも小さくなるように設定する。
 第1輪郭は、対象物領域の輪郭と相似になるように対象物領域の輪郭から常に略一定の距離内側(対象部領域側)に離れたものが好ましく、第2輪郭は、対象物領域の輪郭と相似になるように対象物領域の輪郭から常に略一定の距離外側(非対象部領域側)に離れたものが好ましいが、予め決められた誤差の範囲内で対象物領域の輪郭から距離が異なるように特定されても良い。
  また、図17A及びBに示すように、輪郭を構成するためには4個以上の画素を用い、少なくとも第1輪郭または第2輪郭に少なくとも1つ以上の画素があれば、輪郭の形状を特定することができる。図17Aに、直線に沿って対象物領域と非対象物領域に分離される例を示し、図17Bに、角のある形状で対象物領域と非対象物領域に分離される例を示している。他にも三角形や矩形など様々な形状を特定することができる。あるいは、面や曲面の形状を特定して、3次元画像中の面や曲面を検出し、分離するようにしてもよい。
 また、上述では、血管のように予め形状が分かっている場合について説明したが、図18に示す第2の画像処理装置1aのように、画像データPを用いて対象物領域の輪郭を推定する輪郭推定手段18を設ける構成にしても良い。
 輪郭推定手段18は、輪郭を推定する手法として、エッジが検出可能な微分フィルタリングを利用することができる。輪郭特定手段13では、輪郭の近傍で微分フィルタのレスポンスが正から負に変わるので、微分フィルタの正のレスポンスを持つ位置を第1輪郭(対象物領域)とし負のレスポンスを持つ位置を第2輪郭(非対象物領域)として特定し、画素選択手段14では第1輪郭と第2輪郭上の画素を選択して、エネルギー設定手段15で正のレスポンスを持つ位置の第1輪郭と負のレスポンスを持つ位置の第2輪郭に組み合わせに対して最小になるようにエネルギーを設定する。以下、前述の通りラベリング手段16を用いてエネルギー最小化により対象領域と非対象領域にセグメンテーションを行なった後、対象領域のみの画像にするなど適切な画像を生成して表示手段17で画像を表示装置上に表示する。
 あるいは、ある程度形状が分かっている場合には、輪郭推定手段18は、所定のパラメータで表現された形状モデルを用いて画像データから対象領域の輪郭を推定することもできる。例えば、ASM(Active Shape Model)のように、対象物領域の形状を主成分分析や独立成分分析等を行って平均形状と平均形状を変形するためのベクトルを求めておき、各ベクトルのパラメータを変えることで、平均形状から様々な形状に変形する手法が知られている。この形状モデルにおいて、各種パラメータを変化させながら、対応する輪郭位置に高次エネルギーを設定するようにしてもよい。このとき、輪郭が柔軟に変化してもよいように、非対象領域(背景)上にある第1輪郭と対象領域上にある第2輪郭の画素を一定の間隔をあけて選択するようにするとよりよい。
 なお、上述の実施例において円を複数の円弧に分割したように、所定のパラメータで表現された形状モデルを複数の線分に分割し、線分ごとにエネルギーを設定するようにしてもよい。
 非特許文献2のRotherらの手法では高次のエネルギーを用いることにより、特定の形状をもつ構造物の領域分割の精度を向上させることを可能にしている。一方、本発明は、上記で説明したように、矩形のウィンドウのように固定せず、画素を柔軟に選択する点に特徴がある。これにより、設定する高次エネルギーを疎にすることができ、少ない計算量で高速に解を得ることが可能である。画素を選択する方法は、セグメンテーション対象の輪郭形状に沿って行うことが好ましい。例えば3次元画像中の血管輪郭は円筒形状をしているから、円形となる画素の組み合わせを選択する。また対象の輪郭形状が統計的なモデル(例えばASM)として学習されている場合は、そのモデルの輪郭に沿った画素を選択する。このように対象形状の輪郭に沿った画素を選択することが肝要であり、対象構造はなんでもよい。
 本実施の形態では、グラフカットの手法をもちいて、エネルギー最小化について説明したが、離散関数を最小化するものであればよい。また、エネルギー最小化の具体的な手法にはグラフカットやQPBO、Tree Reweighted Message Passing(詳細は、V. Kolmogorov, Convergent Tree-reweighted Message Passing for Energy Minimization, IEEE PAMI, 28(10), October 2006を参照)などがある。また、グラフカットを用いて対象物領域と非対象物領域のそれぞれについてラベル付けする方法については、‘Y. Boykov and V. Kolmogorov, An Experimental Comparison of Min-Cut/Max-Flow Algorithms for Energy Minimization in Vision, IEEE PAMI, 26(9), pp. 1124-1137, 2004’に詳細に説明されている。QPBOを用いて対象物領域と非対象物領域のそれぞれについてラベル付けする方法は、非特許文献2に詳細に説明されている。

Claims (9)

  1.  エネルギー最小化手法を用いて、画像データを構成する各画素を対象物領域を表す第1の変数または前記対象物領域を除く非対象領域を表す第2の変数にラベル付けをする画像処理装置であって、
     前記対象物領域の輪郭に相似した形状を有する、前記対象物領域側に存在する第1輪郭、および、前記対象物領域の輪郭を挟んで前記第1輪郭と対向した位置に存在する前記非対象物領域側に存在する第2輪郭、を特定する輪郭特定手段と、
     前記第1輪郭の全体を形成するi(≧1)個の画素と、前記第2輪郭の全体を形成する(N-i)(N≧4、N-i≧1)個の画素を選択する画素選択手段と、
     前記第1輪郭の前記i個の画素が全て前記対象物領域に属し、かつ、前記第2輪郭の(N-i)個の画素が全て前記非対象領域に属するという条件を満たす場合のN次のエネルギーを、前記条件を満たさない場合のN次のエネルギーよりも小さくなるように設定するエネルギー設定手段と、
     前記N次のエネルギーを最小化してラベル付けをするラベリング手段とを備える画像処理装置。
  2.  前記輪郭特定手段が、前記第1輪郭と前記第2輪郭とが所定の距離離間するように第1輪郭および第2輪郭を特定するものであることを特徴とする請求項1記載の画像処理装置。
  3.  前記画像データを用いて前記対象物領域の輪郭を推定する輪郭推定手段をさらに備え、
     前記輪郭特定手段が、前記推定した輪郭に基づいて、前記第1輪郭および前記第2輪郭を特定することを特徴とする請求項1記載の画像処理装置。
  4.  前記輪郭推定手段が、エッジを検出可能な微分フィルタリングを用いて前記画像データから対象物領域の輪郭を推定することを特徴とする請求項3記載の画像処理装置。
  5.  前記輪郭推定手段が、所定のパラメータで表現された形状モデルを用いて前記画像データから前記対象領域の輪郭を推定することを特徴とする請求項3記載の画像処理装置。
  6.  前記エネルギー設定手段が、前記条件を満たさない場合のN次のエネルギーを全て等しい値に設定することを特徴とする請求項1記載の画像処理装置。
  7.  前記ラベリング手段が、エネルギーの最小化のためにQPBOアルゴリズムを用いることを特徴とする請求項1記載の画像処理装置。
  8.  コンピュータを、エネルギー最小化手法を用いて、画像データを構成する各画素を対象物領域を表す第1の変数または前記対象物領域を除く非対象領域を表す第2の変数にラベル付けをする画像処理装置として機能させるための画像処理プログラムであって、
     前記対象物領域の輪郭に相似した形状を有する、前記対象物領域側に存在する第1輪郭、および、前記対象物領域の輪郭を挟んで前記第1輪郭と対向した位置に存在する前記非対象物領域側に存在する第2輪郭、を特定する輪郭特定手段と、
     前記第1輪郭の全体を形成するi(≧1)個の画素と、前記第2輪郭の全体を形成する(N-i)(N≧4、N-i≧1個の画素を選択する画素選択手段と、
     前記第1輪郭の前記i個の画素が全て前記対象物領域に属し、かつ、前記第2輪郭の(N-i)個の画素が全て前記非対象領域に属するという条件を満たす場合のN次のエネルギーを、前記条件を満たさない場合のN次のエネルギーよりも小さくなるように設定するエネルギー設定手段と、
     前記N次のエネルギーを最小化してラベル付けをするラベリング手段として機能させるものであることを特徴とする画像処理プログラム。
  9.  輪郭特定手段と、画素選択手段と、エネルギー設定手段と、ラベリング手段とを有し、エネルギー最小化手法を用いて、画像データを構成する各画素を対象物領域を表す第1の変数または前記対象物領域を除く非対象領域を表す第2の変数にラベル付けをする画像処理装置の作動方法であって、
     前記輪郭特定手段により、前記対象物領域の輪郭に相似した形状を有する、前記対象物領域側に存在する第1輪郭、および、前記対象物領域の輪郭を挟んで前記第1輪郭と対向した位置に存在する前記非対象物領域側に存在する第2輪郭、を特定する輪郭特定ステップと、
     前記画素選択手段により、前記第1輪郭の全体を形成するi(≧1)個の画素と、前記第2輪郭の全体を形成する(N-i)(N≧4、N-i≧1)個の画素を選択する画素選択ステップと、
     前記エネルギー設定手段により、前記第1輪郭の前記i個の画素が全て前記対象物領域に属し、かつ、前記第2輪郭の(N-i)個の画素が全て前記非対象領域に属するという条件を満たす場合のN次のエネルギーを、前記条件を満たさない場合のN次のエネルギーよりも小さくなるように設定するエネルギー設定ステップと、
     前記ラベリング手段により、前記N次のエネルギーを最小化してラベル付けをするラベリングステップを備えた画像処理装置の作動方法。
PCT/JP2014/001618 2013-03-27 2014-03-20 画像処理装置、画像処理プログラムおよび画像処理装置の作動方法 Ceased WO2014156089A1 (ja)

Priority Applications (2)

Application Number Priority Date Filing Date Title
DE112014001697.7T DE112014001697T5 (de) 2013-03-27 2014-03-20 Bildverarbeitungsvorrichtung, Bildverarbeitungsprogramm und Betriebsverfahren für Bildverarbeitungsvorrichtung
US14/865,616 US9965698B2 (en) 2013-03-27 2015-09-25 Image processing apparatus, non-transitory computer-readable recording medium having stored therein image processing program, and operation method of image processing apparatus

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2013-066092 2013-03-27
JP2013066092A JP5879291B2 (ja) 2013-03-27 2013-03-27 画像処理装置、画像処理プログラムおよび画像処理装置の作動方法

Related Child Applications (1)

Application Number Title Priority Date Filing Date
US14/865,616 Continuation US9965698B2 (en) 2013-03-27 2015-09-25 Image processing apparatus, non-transitory computer-readable recording medium having stored therein image processing program, and operation method of image processing apparatus

Publications (1)

Publication Number Publication Date
WO2014156089A1 true WO2014156089A1 (ja) 2014-10-02

Family

ID=51623103

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2014/001618 Ceased WO2014156089A1 (ja) 2013-03-27 2014-03-20 画像処理装置、画像処理プログラムおよび画像処理装置の作動方法

Country Status (4)

Country Link
US (1) US9965698B2 (ja)
JP (1) JP5879291B2 (ja)
DE (1) DE112014001697T5 (ja)
WO (1) WO2014156089A1 (ja)

Families Citing this family (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP6626344B2 (ja) * 2015-09-29 2019-12-25 キヤノン株式会社 画像処理装置、画像処理装置の制御方法およびプログラム
US10002428B2 (en) * 2015-12-01 2018-06-19 Ottawa Hospital Research Institute Method and system for identifying bleeding
JP7401739B2 (ja) * 2019-10-04 2023-12-20 富士通株式会社 次数変換装置、次数変換方法、および次数変換プログラム
CN113450368B (zh) * 2020-03-27 2025-03-11 阿里巴巴集团控股有限公司 图像处理方法及装置
JP7509585B2 (ja) 2020-06-16 2024-07-02 株式会社日立製作所 制約用イジングモデル生成システム、組合せ最適化計算システム及び制約用イジングモデル計算方法
CN112766272B (zh) * 2021-01-15 2024-10-15 北京迈格威科技有限公司 目标检测方法、装置和电子系统

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2006053919A (ja) * 2004-08-06 2006-02-23 Microsoft Corp 画像データ分離システム及びその方法
JP2010121945A (ja) * 2008-11-17 2010-06-03 Nippon Hoso Kyokai <Nhk> 3次元形状生成装置
JP2013020600A (ja) * 2011-06-17 2013-01-31 Denso Corp 画像処理装置

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
GB0305315D0 (en) * 2003-03-07 2003-04-09 Weber Martin Image processing system
JP4184842B2 (ja) * 2003-03-19 2008-11-19 富士フイルム株式会社 画像判別装置、方法およびプログラム
JP4717935B2 (ja) 2009-03-23 2011-07-06 富士フイルム株式会社 画像処理装置および方法並びにプログラム
JP4709290B2 (ja) 2009-03-03 2011-06-22 富士フイルム株式会社 画像処理装置および方法並びにプログラム
JP5505164B2 (ja) 2010-07-23 2014-05-28 ソニー株式会社 画像処理装置および方法、並びにプログラム
JP5716170B2 (ja) * 2010-07-26 2015-05-13 石川 博 情報処理方法および情報処理装置

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2006053919A (ja) * 2004-08-06 2006-02-23 Microsoft Corp 画像データ分離システム及びその方法
JP2010121945A (ja) * 2008-11-17 2010-06-03 Nippon Hoso Kyokai <Nhk> 3次元形状生成装置
JP2013020600A (ja) * 2011-06-17 2013-01-31 Denso Corp 画像処理装置

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
HARUYUKI IWAMA ET AL.: "Foreground, Shadow and Background Segmentation Based on Homography- Correspondence Pair", THE TRANSACTIONS OF THE INSTITUTE OF ELECTRONICS, INFORMATION AND COMMUNICATION ENGINEERS, vol. J94-D, no. 8, 1 August 2011 (2011-08-01), pages 1300 - 1313 *
PHAM VIET-QUOC ET AL.: "Image Segmentation from Bounding Shape without Using Appearance Model", THE TRANSACTIONS OF THE INSTITUTE OF ELECTRONICS, INFORMATION AND COMMUNICATION ENGINEERS, vol. J94-D1, no. 8, 1 August 2011 (2011-08-01), pages 1183 - 1193 *

Also Published As

Publication number Publication date
JP2014191564A (ja) 2014-10-06
JP5879291B2 (ja) 2016-03-08
US9965698B2 (en) 2018-05-08
DE112014001697T5 (de) 2015-12-17
US20160019435A1 (en) 2016-01-21

Similar Documents

Publication Publication Date Title
US8131076B2 (en) Editing of pre-segmented images using seeds derived from contours
US7400767B2 (en) System and method for graph cuts image segmentation using a shape prior
JP5879291B2 (ja) 画像処理装置、画像処理プログラムおよび画像処理装置の作動方法
US8787642B2 (en) Method, device and computer-readable recording medium containing program for extracting object region of interest
JP6539303B2 (ja) 3d医用画像中の対象物を分割するための3d対象物の変換
US20100266175A1 (en) Image and data segmentation
CN112927353A (zh) 基于二维目标检测和模型对齐的三维场景重建方法、存储介质及终端
Shapovalov et al. Spatial inference machines
JP2020109614A (ja) 画像処理装置、画像処理システム、画像処理方法、プログラム
CN107545579B (zh) 一种心脏分割方法、设备和存储介质
US8290266B2 (en) Model-based method and system for image segmentation and modelling
CN102663762A (zh) 医学图像中对称器官的分割方法
Hosseini-Asl et al. Lung segmentation based on nonnegative matrix factorization
EP2902970B1 (en) Image processing device, image processing program, and operating method for image processing device
JP5897445B2 (ja) 分類装置、分類プログラムおよび分類装置の動作方法
Choong et al. Multistage image clustering and segmentation with normalised cuts
EP3871148B1 (en) Method and device for training a neural network to specify landmarks on 2d and 3d images
CN116630212B (zh) 一种基于条件gan网络的自适应特征融合的数据合成方法
CN110751658A (zh) 一种基于互信息和点扩散函数的抠图方法
Benzian et al. 3D Mesh Segmentation by Region growing based on discrete curvature
Michel et al. Arm-nms: Shape based non-maximum suppression for instance segmentation in large scale imagery
Chen et al. Pulmonary nodule segmentation in computed tomography with an encoder-decoder architecture
Denk et al. Feature line detection of noisy triangulated CSGbased objects using deep learning
CN121305097B (zh) 一种适用于桥梁点云的构件自动化分割识别方法
JP2006277713A (ja) 3次元メッシュモデルの特徴稜線抽出装置、プログラム及び方法

Legal Events

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

Ref document number: 14774489

Country of ref document: EP

Kind code of ref document: A1

WWE Wipo information: entry into national phase

Ref document number: 1120140016977

Country of ref document: DE

Ref document number: 112014001697

Country of ref document: DE

122 Ep: pct application non-entry in european phase

Ref document number: 14774489

Country of ref document: EP

Kind code of ref document: A1