WO2019151302A1 - 画像処理方法、プログラムおよび記録媒体 - Google Patents

画像処理方法、プログラムおよび記録媒体 Download PDF

Info

Publication number
WO2019151302A1
WO2019151302A1 PCT/JP2019/003106 JP2019003106W WO2019151302A1 WO 2019151302 A1 WO2019151302 A1 WO 2019151302A1 JP 2019003106 W JP2019003106 W JP 2019003106W WO 2019151302 A1 WO2019151302 A1 WO 2019151302A1
Authority
WO
WIPO (PCT)
Prior art keywords
image
image processing
embryo
dimensional image
processing method
Prior art date
Application number
PCT/JP2019/003106
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 US16/966,488 priority Critical patent/US11710237B2/en
Priority to EP19747001.6A priority patent/EP3748580B1/en
Priority to CN201980009540.7A priority patent/CN111630564B/zh
Publication of WO2019151302A1 publication Critical patent/WO2019151302A1/ja

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/0002—Inspection of images, e.g. flaw detection
    • G06T7/0012—Biomedical image inspection
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N15/00—Investigating characteristics of particles; Investigating permeability, pore-volume or surface-area of porous materials
    • G01N15/10—Investigating individual particles
    • G01N15/14—Optical investigation techniques, e.g. flow cytometry
    • G01N15/1429—Signal processing
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
    • G01N21/17—Systems in which incident light is modified in accordance with the properties of the material investigated
    • 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/12—Edge-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/155—Segmentation; Edge detection involving morphological operators
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/60—Analysis of geometric attributes
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/0059—Measuring for diagnostic purposes; Identification of persons using light, e.g. diagnosis by transillumination, diascopy, fluorescence
    • A61B5/0062—Arrangements for scanning
    • A61B5/0066—Optical coherence imaging
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2200/00—Indexing scheme for image data processing or generation, in general
    • G06T2200/04—Indexing scheme for image data processing or generation, in general involving 3D image data
    • 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/10101—Optical tomography; Optical coherence tomography [OCT]
    • 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/30024—Cell structures in vitro; Tissue sections in vitro
    • 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/30044—Fetus; Embryo

Definitions

  • the present invention relates to an image processing technique for providing useful information to a user who evaluates an embryo from data corresponding to a three-dimensional image of the embryo obtained by tomographic imaging of the embryo.
  • Patent Document 1 proposes a technique for evaluating embryo quality by applying a non-invasive tomographic imaging technique such as optical coherence tomography (Optical Coherence Tomography; OCT).
  • OCT optical Coherence Tomography
  • a three-dimensional image of an embryo is reconstructed from a tomographic image obtained by OCT imaging. Based on the three-dimensional image, the quality of the embryo is evaluated using criteria such as the number of blastomeres, regularity, and fragmentation rate as the morphological criteria of the embryo.
  • the condition of a good quality embryo is that the size and shape of the cells contained in the embryo are well aligned.
  • a quantitative index for individually evaluating the size and shape of each individual cell constituting the embryo is obtained from the imaging data.
  • the above prior arts are exemplified by a display mode in which a transparent body on the surface of the embryo is erased to expose internal cells, and a display mode in which it is rotated so that it can be observed from various directions (for example, FIG. 5A to FIG. 5I) do not contribute to the above-described quantification.
  • the tomographic image obtained by OCT imaging and the three-dimensional image reconstructed therefrom have the potential to provide a lot of information useful for the evaluation of the embryo, such as the internal structure of the embryo. doing.
  • a specific data processing method for that purpose has not been established so far.
  • the present invention has been made in view of the above problems, and provides a data processing method suitable for obtaining quantitative information from data obtained by OCT imaging, thereby effectively supporting an embryo evaluation work by a user.
  • the purpose is to do.
  • Image processing for obtaining an index value for indexing the size, dividing the three-dimensional image into an area where the index value is larger than a predetermined first threshold and a smaller area, and dividing each of the divided areas by a watershed algorithm Is the method.
  • Another aspect of the present invention is a program for causing a computer to execute the above steps.
  • Another aspect of the present invention is a computer-readable recording medium storing the above program.
  • embryos particularly mammalian fertilized eggs, have a generally spherical outer shape, and the inside has a structure in which a plurality of cells that are generally spherical or ellipsoidal are gathered.
  • a spatial region in other words, a boundary between cells in a three-dimensional space.
  • Such demarcation can be performed by region segmentation processing for a three-dimensional image such as an OCT image.
  • a watershed algorithm that is a known image processing algorithm can be used.
  • the result of area division by the Watershed algorithm often deviates from the cell boundary determined by an expert from a three-dimensional image.
  • the boundary of the cells gradually become less clear, and the minute structures (such as polar bodies and incomplete cell division) ( For example, fragmentation) and image noise are considered to cause misdivision.
  • a method of image processing called LocalnessThickness calculation is used (particularly, a three-dimensional image).
  • 3D Local Thickness calculation the object is expressed in gradation based on the radius of the sphere inscribed in the object in the three-dimensional image. For example, an object having a relatively large structure is represented by a high gradation. On the other hand, relatively small objects are represented with lower gradations.
  • a value that depends on the size of an object in this case, a gradation value
  • a relatively large object corresponding to a cell or the like and a minute object corresponding to a fragment or the like can be distinguished. It becomes.
  • the boundary between cells becomes clear. For this reason, if the region division by the Watershed algorithm is further performed on the region thus divided, a division result that matches the boundary of an actual cell or the like can be obtained.
  • index values that quantitatively indicate the shape, size, etc. are obtained by individually evaluating relatively large ones of the areas divided by the above image processing method as individual cells. Is possible. For example, a region corresponding to a cell can be distinguished from other regions depending on whether the index value obtained by the Local Thickness calculation is larger or smaller than a predetermined first threshold value.
  • a three-dimensional image of an embryo is divided into a plurality of regions by the above-described image processing method, and the index value of the plurality of regions is larger than the first threshold value.
  • This is an image processing method for determining an area as an aggregate of blastomeres. This makes it possible to distinguish cells that are relatively large structures in the embryo from relatively small structures such as polar bodies and fragments. In addition, it is possible to quantitatively determine the characteristics of each three-dimensional shape, etc., and use it for embryo evaluation.
  • This embodiment embodies an image processing method according to the present invention.
  • an embryo cultured in a medium (culture medium) carried on a culture vessel is tomographically imaged by optical coherence tomography (OCT) technology, and based on the data obtained thereby, Quantitative information on the cells that make up the embryo is output.
  • OCT optical coherence tomography
  • the purpose is to assist the user (specifically, a doctor or embryo cultivator) with the embryo evaluation work.
  • the image processing method of this embodiment can be applied for the purpose of obtaining knowledge for determining whether or not the culture is progressing well.
  • FIG. 1 is a principle diagram showing a configuration example of an image processing apparatus for performing tomographic imaging.
  • the image processing apparatus 1 performs tomographic imaging using an embryo cultured in a medium as an imaging object, and performs image processing on the obtained tomographic image to create a three-dimensional (3D) image of the imaging object.
  • 3D three-dimensional
  • xyz orthogonal coordinate axes are set as shown in FIG.
  • the xy plane represents a horizontal plane.
  • the z axis represents the vertical axis, and more specifically, the ( ⁇ z) direction represents the vertical downward direction.
  • the image processing apparatus 1 includes a holding unit 10.
  • the holding unit 10 holds the culture vessel 11 in a substantially horizontal posture with its opening surface facing upward.
  • the culture container (hereinafter simply referred to as “container”) 11 is a dish or plate having a transparent bottom made of glass or resin.
  • a predetermined amount of an appropriate medium M is injected into the container 11 in advance, and a sample Sp (embryo in this case) is cultured on the bottom 111 of the container 11 in the medium.
  • a sample Sp is illustrated in FIG. 1, a plurality of samples Sp may be cultured in one container 11.
  • the imaging unit 20 is disposed below the container 11 held by the holding unit 10.
  • the imaging unit 20 is an optical coherence tomography (OCT) apparatus that can capture a tomographic image of an object to be imaged in a non-contact and non-destructive (non-invasive) manner.
  • OCT optical coherence tomography
  • an imaging unit 20 that is an OCT apparatus includes a light source 21 that generates illumination light for an object to be imaged, a beam splitter 22, an object optical system 23, a reference mirror 24, a spectroscope 25, And a photodetector 26.
  • the image processing apparatus 1 further includes a control unit 30 that controls the operation of the apparatus and a drive control unit 40 that controls the movable mechanism of the imaging unit 20.
  • the control unit 30 includes a central processing unit (CPU) 31, an A / D converter 32, a signal processing unit 33, a 3D restoration unit 34, an interface (IF) unit 35, an image memory 36, and a memory 37.
  • the CPU 31 controls the operation of the entire apparatus by executing a predetermined control program.
  • the control program executed by the CPU 31 and data generated during processing are stored in the memory 37.
  • the A / D converter 32 converts a signal output from the photodetector 26 of the imaging unit 20 according to the amount of received light into digital data.
  • the signal processing unit 33 performs signal processing to be described later based on the digital data output from the A / D converter 32 to create a tomographic image of the object to be imaged.
  • the 3D reconstruction unit 34 has a function of creating a captured three-dimensional image (3D image) of the embryo based on the image data of the plurality of captured tomographic images.
  • the image data of the tomographic image created by the signal processing unit 33 and the image data of the three-dimensional image created by the 3D restoration unit 34 are stored and saved as appropriate by the image memory 36.
  • the interface unit 35 is responsible for communication between the image processing apparatus 1 and the outside. Specifically, the interface unit 35 has a communication function for communicating with an external device, and a user interface function for receiving an operation input from the user and notifying the user of various information.
  • an input device 351 and a display unit 352 are connected to the interface unit 35.
  • the input device 351 includes, for example, a keyboard, a mouse, a touch panel, and the like that can accept operation inputs related to device function selection, operation condition setting, and the like.
  • the display unit 352 includes a display device that displays various processing results such as a tomographic image created by the signal processing unit 33 and a stereoscopic image created by the 3D restoration unit 34, such as a liquid crystal display.
  • the CPU 31 gives a control command to the drive control unit 40.
  • the drive control unit 40 causes the movable mechanism of the imaging unit 20 to perform a predetermined operation.
  • the sample Sp (specifically, the embryo) is captured by a combination of the scanning movement of the imaging unit 20 executed by the drive control unit 40 and the detection of the amount of received light by the photodetector 26. ) Is acquired.
  • FIG. 2A and 2B are diagrams for explaining the imaging principle of the image processing apparatus. More specifically, FIG. 2A is a diagram showing an optical path in the imaging unit 20, and FIG. 2B is a diagram schematically showing a state of tomographic imaging of the sample Sp.
  • the imaging unit 20 functions as an optical coherence tomography (OCT) apparatus.
  • OCT optical coherence tomography
  • the sample Sp is displayed as a substantially spherical shape.
  • the shape of the sample Sp is not particularly limited in actual imaging.
  • a low coherence light beam L1 including a broadband wavelength component is emitted from a light source 21 having a light emitting element such as a light emitting diode or a super luminescent diode (SLD).
  • the light beam L1 enters the beam splitter 22 and branches. A part of the light L2 is directed to the container 11 as indicated by a broken line arrow, and a part of the light L3 is directed to the reference mirror 24 as indicated by a one-dot chain line arrow.
  • the light L2 directed toward the container 11 enters the container 11 through the object optical system 23. More specifically, the light L ⁇ b> 2 emitted from the beam splitter 22 enters the container bottom 111 through the object optical system 23.
  • the object optical system 23 has a function of converging light L ⁇ b> 2 from the beam splitter 22 toward the container 11 onto the sample Sp in the container 11, and a function of condensing the reflected light emitted from the sample Sp toward the beam splitter 22.
  • the object optical system 23 is typically represented by a single objective lens, but may be a combination of a plurality of optical elements.
  • the object optical system 23 is movable in the z direction by a focus adjustment mechanism 41 provided in the drive control unit 40. Thereby, the focal position of the object optical system 23 with respect to the object to be imaged can be changed in the z direction.
  • the optical axis of the object optical system 23 is parallel to the vertical direction, and is therefore perpendicular to the planar container bottom 111.
  • the incident light is incident on the object optical system 23 in parallel with the optical axis, and the arrangement of the object optical system 23 is determined so that the optical center of the illumination light coincides with the optical axis.
  • the sample Sp is not transparent to the light L2, the light L2 incident through the container bottom 111 is reflected by the surface of the sample Sp.
  • the sample Sp has a certain degree of transparency with respect to the light L2, the light L2 enters the sample Sp and is reflected by the internal structure. If, for example, near infrared rays are used as the light L2, incident light can reach the inside of the sample Sp. The reflected light from the sample Sp is emitted in various directions as scattered light. Of these, the light L4 emitted within the condensing range of the object optical system 23 is condensed by the object optical system 23 and sent to the beam splitter 22.
  • the reference mirror 24 is supported on its reflecting surface in a posture perpendicular to the incident direction of the light L3. Further, the mirror drive mechanism 42 provided in the drive control unit 40 can move in the direction along the incident direction of the light L3 (the y direction in the figure). The light L3 incident on the reference mirror 24 is reflected by the reflecting surface and travels toward the beam splitter 22 as light L5 that travels in the reverse direction along the incident optical path. This light L5 becomes the reference light. When the position of the reference mirror 24 is changed by the mirror driving mechanism 42, the optical path length of the reference light changes.
  • interference due to the phase difference occurs between the reflected light L4 and the reference light L5, but the spectral spectrum of the interference light varies depending on the depth of the reflecting surface. That is, the spectral spectrum of the interference light has information in the depth direction of the object to be imaged. Therefore, the reflected light intensity distribution in the depth direction of the object to be imaged can be obtained by spectrally detecting the interference light for each wavelength to detect the amount of light and performing Fourier transform on the detected interference signal.
  • the OCT imaging technique based on such a principle is called Fourier domain OCT (FD-OCT).
  • a spectroscope 25 is provided on the optical path of the interference light from the beam splitter 22 to the photodetector 26.
  • the spectroscope 25 for example, one using a prism, one using a diffraction grating, or the like can be used.
  • the interference light is spectrally separated for each wavelength component by the spectroscope 25 and received by the photodetector 26.
  • the interference signal output from the light detector 26 according to the interference light detected by the light detector 26 is Fourier-transformed to reflect the sample Sp in the depth direction at the incident position of the light beam L2, that is, the z-direction reflection.
  • a light intensity distribution is required.
  • a reflected light intensity distribution in a plane parallel to the xz plane is obtained. From the result, a tomographic image of the sample Sp having a cross section of the plane can be created.
  • a series of operations for acquiring one tomographic image It in a cross section parallel to the xz plane by beam scanning in the x direction is referred to as one imaging.
  • a tomographic image is taken each time the beam incident position in the y direction is changed in multiple stages. By doing so, as shown in FIG. 2B, it is possible to obtain a large number of tomographic images It obtained by tomographically imaging the sample Sp in a cross section parallel to the xz plane. If the scanning pitch in the y direction is reduced, image data with sufficient resolution to grasp the three-dimensional structure of the sample Sp can be obtained.
  • the beam scanning in the x direction and the y direction can be realized by various methods. For example, a method of changing the beam incident position in the xy direction using an optical component that changes an optical path such as a galvanometer mirror (not shown) can be applied. Further, it is also possible to apply a method in which one of the container 11 carrying the sample Sp and the imaging unit 20 is moved in the xy direction to change the relative positions thereof.
  • the imaging unit 20 has a demultiplexing function for branching light from the light source 21 into illumination light and reference light, and a function for combining signal light and reference light to generate interference light. This is realized by the beam splitter 22.
  • optical fiber couplers as exemplified below may be used as those having such a demultiplexing / multiplexing function in the OCT apparatus.
  • 3A and 3B are diagrams showing a specific configuration example of the OCT apparatus.
  • the same reference numerals are given to configurations that are the same as or correspond to the configurations in the principle diagram described above. Its structure and function are basically the same as those in the above principle diagram unless otherwise specified, and detailed description thereof will be omitted. Further, since the OCT imaging principle for detecting interference light by the optical fiber coupler is basically the same as described above, detailed description is omitted.
  • the imaging unit 20a includes an optical fiber coupler 220 as a demultiplexer / multiplexer in place of the beam splitter 22.
  • One optical fiber 221 constituting the optical fiber coupler 220 is connected to the light source 21, and the low coherence light emitted from the light source 21 is branched by the optical fiber coupler 220 into light to the two optical fibers 222 and 223.
  • the optical fiber 222 constitutes an object system optical path. More specifically, the light emitted from the end portion of the optical fiber 222 enters the object optical system 23 via the collimator lens 223. Reflected light (signal light) from the imaging object enters the optical fiber 222 via the object optical system 23 and the collimator lens 223.
  • the other optical fiber 224 constitutes a reference system optical path. More specifically, light emitted from the end of the optical fiber 224 enters the reference mirror 24 via the collimator lens 225. The reflected light (reference light) from the reference mirror 24 enters the optical fiber 224 via the collimator lens 225. The signal light propagating through the optical fiber 222 interferes with the reference light propagating through the optical fiber 224 in the optical fiber coupler 220. The interference light enters the photodetector 26 via the optical fiber 226 and the spectroscope 25. As described above, the intensity distribution of the reflected light on the object to be imaged is obtained from the interference light received by the photodetector 26.
  • the optical fiber coupler 220 is provided in the imaging unit 20b.
  • the optical fiber 224 is not used, and a collimator lens 223 and a beam splitter 227 are provided for the optical path of light emitted from the optical fiber 222.
  • the object optical system 23 and the reference mirror 24 are arrange
  • the focus adjustment mechanism 41 can adjust the depth of focus of the object optical system 23 relative to the object to be imaged by moving the object optical system 23 in the approaching / separating direction with respect to the container 11. .
  • the optical path length of the reference light can be changed by the mirror driving mechanism 42 moving the reference mirror 24 along the light incident direction.
  • the configuration of the image processing apparatus 1 can be applied to either the above-described beam splitter or the one using an optical fiber coupler.
  • the imaging apparatus for capturing a tomographic image is not limited to the above-described FD-OCT imaging apparatus.
  • a device based on another imaging principle such as a time domain OCT (TD-OCT) imaging device is applicable.
  • TD-OCT time domain OCT
  • reflected light of illumination light incident on each position in the three-dimensional space including the sample Sp and its periphery is used as signal light.
  • strength of the signal corresponding to the light quantity of the signal light from the said position are used for a process. Therefore, any imaging method that can acquire such data may be used.
  • FIGS. 4A to 4D are schematic views illustrating typical structures of embryos to be processed.
  • the image processing apparatus 1 When the image processing apparatus 1 is used, for example, for the purpose of infertility treatment, an in vitro fertilized human fertilized egg at the initial stage of culture becomes the imaging object of the image processing apparatus 1.
  • cleavage is started when an egg is fertilized, and a blastocyst is formed through a state called a morula.
  • the image processing method of this embodiment provides a user with information suitable for observing an embryo in the blastocyst stage immediately after fertilization, for example.
  • FIG. 4A schematically shows the structure of an embryo at an early stage (for example, from the 4-cell stage to the morula stage).
  • the outer shape of the embryo E1 is generally spherical. Its surface is covered with a jelly-like glycoprotein layer Z called zona pellucida.
  • the inside contains a plurality of cells C generated by cell division of a fertilized egg.
  • the inside of the zona pellucida Z is occupied by a plurality of cells C of equal size, and the number increases as the cleavage proceeds. According to the Veeck classification, it is desirable that the cells are aligned.
  • FIG. 4B shows a state in which the variation in the size of the cell C and the fragment F are generated, but only one of them may appear.
  • the state of the embryo is inferior to that shown in FIG. 4A. Therefore, in early embryos up to morula, whether the size of each cell C formed by cleavage is uniform, whether the surface of each cell C is smooth, whether fragment F is generated, etc. It can be an index for evaluating the state.
  • FIG. 4C and 4D schematically show the internal structure of the embryo at the blastocyst stage.
  • FIG. 4C shows embryo E3 in a good state
  • FIG. 4D shows embryo E4 in a poorer state.
  • the cells C having advanced cleavage are arranged as a thin layer on the surface of the embryo to form a trophectoderm T.
  • a cavity called blastocoel B is generated in the internal space surrounded by the trophectoderm T.
  • an inner cell mass I in which a large number of cells C are densely formed is formed in a part of the internal space.
  • the trophectoderm T is densely formed by many cells C, and the inner cell mass I is also densely constituted by relatively many cells C.
  • the vegetative ectoderm T may be formed sparsely with fewer cells C, the inner cell mass I may be composed of small cells and become small. Therefore, the culture state of the embryo can be indicated by the thickness and density of the vegetative ectoderm T and the size and density of the inner cell mass I.
  • FIG. 5 is a flowchart showing the image processing method of this embodiment.
  • step S101 is performed in an appropriate incubator prepared separately from the image processing apparatus 1.
  • the processing after step S102 is realized by the CPU 31 provided in the image processing apparatus 1 executing a control program recorded in advance in the memory 37, but is not limited thereto.
  • step S103 and the subsequent steps in FIG. 5 can be executed by a computer device having a general calculation function and an image output function.
  • These processes may be configured such that a computer apparatus different from the image processing apparatus 1 receives and processes OCT imaging data from the image processing apparatus 1.
  • the CPU 31 of the image processing apparatus 1 only needs to execute processing specialized for imaging. As a result, the processing load on the CPU 31 is reduced.
  • an embryo to be evaluated is cultured in an appropriate culture environment (step S101).
  • the cultured embryo is subjected to OCT imaging by the image processing apparatus 1 at an appropriate timing (step S102).
  • Data obtained by imaging specifically, data representing reflected light intensity from each position in the three-dimensional space including the embryo and its surroundings is stored in the memory 37 or the image memory 36 as original signal data.
  • Using this original signal data by creating pixel data (voxel data) for each coordinate position in the three-dimensional space, a three-dimensional image of the embryo as shown in Patent Document 1 is constructed and output as an image Is possible.
  • the created voxel data is stored in the image memory 36 as original image data.
  • the three-dimensional image is binarized for later processing (step S103). For example, among the pixels represented by the voxel data, a pixel whose pixel value (luminance value) exceeds a predetermined threshold is represented by “1”, and a pixel whose pixel value is smaller than the threshold is represented by “0”. Thereby, a three-dimensional image is binarized according to the brightness
  • a region having a high reflected light intensity that is, a high density of a substance that reflects light, is expressed with high luminance in the three-dimensional space included in the imaging range.
  • a region where the density of the substance is low and the reflected light intensity is low is expressed with low luminance.
  • the three-dimensional image of the embryo can include structures such as zona pellucida, trophectoderm, and fragments in addition to the blastomeres formed by cleavage. These need to be identified in order to evaluate the embryo in detail.
  • boundaries between structures included in the three-dimensional image are defined.
  • an image processing algorithm that executes region division in a three-dimensional space, for example, a known Watershed algorithm.
  • the boundary between the convex part and the concave part of the binarized three-dimensional image is specified as the boundary of the region.
  • the image of the embryo E2 includes a plurality of cells C formed by cleavage, a fragment F distributed around the cell C, and a zona pellucida Z surrounding them.
  • structure internal space the space area occupied by the structure in the binarized three-dimensional image
  • structure surface The surface of the structure in the three-dimensional image, that is, the curved surface that is the boundary between the structure internal space and the surrounding space.
  • processing will be described using a two-dimensional image obtained by cutting a three-dimensional image along one cross section. However, actual image processing is performed on a three-dimensional image represented by voxel data.
  • FIG. 6A to 6D are diagrams showing the principle of processing for separating a large structure and a small structure.
  • FIG. 6A is a two-dimensional image corresponding to the binarized image of embryo E2 shown in FIG. 4B.
  • the binarized image is modeled as a collection of inscribed spheres (inscribed circle in the two-dimensional image).
  • An area represented by an inscribed sphere having a radius larger than an appropriate threshold (the “first threshold” of the present invention) hereinafter referred to as “first area” and an inscribed sphere having a radius smaller than the threshold.
  • first threshold an appropriate threshold
  • second regions hereinafter referred to as “second regions”.
  • the radius of the inscribed sphere is equal to the threshold is included in the first region or the second region.
  • relatively large structures such as cells C and their aggregates are represented as a collection of relatively large radius inscribed spheres.
  • small structures such as zona pellucida Z and fragment F are represented as a collection of inscribed spheres with smaller radii. Therefore, the object is achieved by dividing the structure according to the radius of the inscribed sphere. For example, if a structure represented by an inscribed sphere having a radius equal to or smaller than the threshold value Rth shown in FIG. 6B is deleted, only a region corresponding to an aggregate of cells having a large structure is extracted as shown in FIG. 6C. Can do.
  • a structure can be modeled by the following operations. Among the points in the structure internal space, the point where the radius of the inscribed sphere inscribed in the structure surface around the point becomes the maximum, in other words, the point having the maximum distance to the structure surface is found.
  • the inscribed sphere at this time can be specified by the position of the center and the size of the radius. Such a point is not necessarily one.
  • a point where the radius of the inscribed sphere inscribed in the structure surface around the point is maximized is found.
  • Local Thickness algorithm has been proposed as a computation method that can obtain the same result in a shorter time (for example, “Computing Local Thickness of 3D Structures with ImageJ”, Robert P. Dougherty and Karl-Heinz Kunzelmann, Microscopy & Microanalysis 2007 Meeting, August 5-9, Broward County, Convention Center, Ft. Lauderdale, Florida), and this embodiment can also be applied.
  • this calculation for each point of the object in the image, a sphere including the point and inscribed in the surface of the structure is searched for having the maximum radius. The higher the inscribed sphere radius, the higher the score is given to the point.
  • the calculation result is output in such a display manner that the higher the score, the higher the gradation value (luminance value).
  • an object in a three-dimensional image is represented by a gradation value corresponding to the radius of a sphere inscribed in the object. Therefore, gradation conversion is performed on a relatively large object to a relatively high gradation value (luminance value), and a relatively small object is converted to a relatively low gradation value (luminance value). If this gradation value (luminance value) is used as an index value for distinguishing objects, it is possible to distinguish relatively large objects such as cells and minute objects such as fragments. In addition to this, the boundary between cells becomes clear, and appropriate division between regions becomes possible.
  • the maximum gradation value Gmax is assigned to the point where the radius of the inscribed sphere is the maximum value Rmax, and the gradation value becomes smaller as the radius of the inscribed sphere becomes smaller.
  • a gradation value is assigned to each point in the body internal space.
  • the cells C formed by cleavage and the aggregates thereof, or the blastomeres and the aggregates thereof have a relatively large structure. That is, it can be determined that the region divided into the first region by the above processing is a region occupied by cells, blastomeres, or aggregates thereof.
  • the region divided into the second regions can be determined as a region occupied by a structure different from the cell aggregate, that is, the zona pellucida Z or the fragment F.
  • the zona pellucida Z is a continuous thin layer covering the surface of the embryo, while the fragment F is a small mass released from the blastomere and its aggregates. From this, among the regions divided into the second regions, for example, a region having a sphericity greater than an appropriate threshold (the “second threshold” in the present invention) can be determined as the fragment F. Further, it is possible to determine that a region having a certain volume or more in a region having a small sphericity is a transparent zone Z, and the other is a pseudo structure due to image noise. In this way, by separating a minute structure from a blastomer and its aggregate, quantitative information such as shape and size can be obtained from the three-dimensional image.
  • the number and size of fragments are information for effectively supporting embryo classification work based on Veeck classification.
  • FIG. 7A and FIG. 7B are diagrams showing examples of application to blastocysts.
  • FIG. 4C in the blastocyst E3, an inner cell mass I which is a relatively large structure is present inside the trophectoderm T which is a thin cell layer.
  • the trophectoderm T is represented as a collection of inscribed spheres with a small radius
  • the inner cell mass I is represented as a collection of larger inscribed spheres. Therefore, as shown by a solid line in FIG. 7B, if a region represented by an inscribed sphere having a large radius is cut out from the binarized image thus modeled, a region corresponding to the inner cell mass I can be extracted. .
  • a region represented by an inscribed sphere having a small radius is cut out, a region corresponding to the trophectoderm T can be extracted.
  • step S104 Arithmetic processing based on the above principle, for example, Local Thickness calculation is performed on the binarized image (step S104).
  • the object in the binarized image is divided into a region (first region) corresponding to the aggregate of blastomeres and another region (second region).
  • first region is the processing target.
  • step S105 by removing the second area where the radius of the inscribed sphere is smaller than the threshold value Rth from the binarized image (step S105), only the first area corresponding to the aggregate of blastomeres remains.
  • Region division processing by the Watershed algorithm is executed on the first region (step S106). Thereby, the individual cells constituting the aggregate of blastomeres are separated. Since the fragment F causing the error is removed, the boundary between individual cells can be accurately defined by the Watershed algorithm. Note that both the LocalnessThickness calculation algorithm and the Watershed algorithm have already been put into practical use as general-purpose image processing software libraries.
  • the size, for example, volume and surface area of each cell can be easily obtained (step S107).
  • the variation in the volume of each cell constituting the aggregate of blastomeres can be used as an index for determining the quality of the embryo based on the Veeck classification.
  • each cell in order to evaluate each cell in more detail, it is possible to introduce an index value corresponding to the shape of the cell.
  • Individual cells are considered to be generally spherical when in an independent state.
  • deformation occurs particularly in a portion in contact with other cells.
  • each region specified as one cell is approximated to an ellipsoid (step S108).
  • the degree of divergence between the surface of each cell and the surface of the approximate ellipsoid is quantified and expressed so as to meet the above requirement. Specifically, it is as follows. Here, consider two cells in contact with each other. However, it is possible to similarly consider a blastomere assembly in which more cells are in contact with each other.
  • FIGS. 8B and 8C are diagrams schematically showing separation of individual cells by area division.
  • 9A to 9D are diagrams schematically showing ellipsoid approximation of cells.
  • FIG. 8A two cells C1 and C2 are an integrated object in the binarized image in a state before region division.
  • the boundary between the two cells C1 and C2 is defined as shown in FIGS. 8B and 8C.
  • the solid line indicates the outline of the surface of the cell exposed without contacting other cells (hereinafter referred to as “exposed surface”)
  • the dotted line indicates the boundary between the cells defined by the region division process.
  • the cells thus divided are considered to be spherical or elliptical as a whole.
  • ellipsoid approximation is performed using only the information of the exposed surface that is not in contact with other cells in the surface of the divided region. By doing so, errors due to cell deformation can be suppressed.
  • the exposed surface is extracted. That is, the binarized image of the divided cell region is expanded by one pixel, and corresponds to the logical product of each pixel with the surrounding region of the aggregate of the blastomeres in the binarized image before the division. Create an image. In this way, an image of the exposed surface having a thickness of one pixel is obtained.
  • the boundary surface created by the division process is included in the original aggregate of blastomeres and therefore does not appear in the result of the logical product.
  • the approximate ellipsoid Ea1 corresponding to the cell C1 approximates the shape of the exposed surface of the cell C1 indicated by the solid line, and is a boundary created by region division indicated by the dotted line. Does not necessarily match. That is, the approximate ellipsoid Ea1 is not affected by the deformation of the cell C1 due to contact with other cells. The same applies to the approximate ellipsoid Ea2 corresponding to the cell C2.
  • the shape of the cell can be quantitatively represented by the parameter representing the shape.
  • the parameters for example, the length of the major axis and the minor axis, the volume, the flatness, the direction of the major axis, and the like can be used. This makes it possible to easily and quantitatively evaluate the shape of individual cells, compare shapes among a plurality of cells, and the like.
  • a numerical value indicating the degree of deviation between the exposed surface of the cell and the surface of the approximate ellipsoid can be a value that quantitatively indicates the state of the cell.
  • each position in the three-dimensional space is represented by the xyz coordinate system.
  • the approximate shape of the cell can be represented by an ellipsoid, it is preferable to represent each position in spherical coordinates with the center of the cell as the origin in order to more express the shape of the cell. Therefore, as shown in FIG. 9B, coordinate conversion from the xyz orthogonal coordinate system to the r ⁇ spherical coordinate system using the radius r and the two declinations ⁇ and ⁇ as coordinate variables is performed.
  • the center of gravity of the approximate ellipsoid obtained previously is set as the origin O of the spherical coordinate system.
  • This origin O does not have to coincide with the origin of the xyz orthogonal coordinate system in the original image data.
  • coordinate conversion from orthogonal coordinates to spherical coordinates is performed by an appropriate conversion process (step S109).
  • each position in the three-dimensional space specified by the xyz orthogonal coordinate system is represented by the r ⁇ spherical coordinate system.
  • the center of gravity of the ellipsoid can be, for example, the midpoint of a line segment representing the major axis or minor axis of the ellipsoid.
  • the coordinates of each point on the exposed surface of the cell C1 (C2) and the surface of the approximate ellipsoid Ea1 (Ea2) are represented by polar coordinates.
  • this distance ⁇ r is This value represents the degree of divergence between the cell and the ellipsoid. That is, if the distance ⁇ r is small, the divergence between the surface of the cell and the surface of the ellipsoid is small, and it can be said that the divergence is remarkable as the distance ⁇ r is large.
  • the distance ⁇ r is obtained in the same manner as described above for the entire exposed surface in the radial direction, and the maximum value, the difference between the maximum value and the minimum value, At least one of the absolute integral values and the like can be used as an index value (step S110).
  • index values represent the degree of cell distortion. The smaller the index value, the smoother the surface of the cell and the less the distortion, and it can be estimated that the cell is a good cell.
  • the evaluation on the shape surface of each structure can be performed based on a binarized image obtained by binarizing a three-dimensional image obtained by OCT imaging.
  • information relating to the density of the three-dimensional object due to the difference in density of the structure is required.
  • the binarized image obtained by dividing individual cells and other structures into individual regions as described above is used as a mask, and this is applied to a three-dimensional image including light and shade represented by the original image data. . Then, the three-dimensional object included in the original image can be cut out individually for each region. That is, by creating an image corresponding to the logical product of each pixel of the original image and the binarized image divided into regions (step S111), an image obtained by dividing the original image into regions can be obtained. A similar result can be obtained by dividing a three-dimensional object in the original image at the boundary defined by the region division described above.
  • the shading is represented by the pixel value (luminance value) of each pixel in the area.
  • the internal structure can be individually evaluated.
  • the average concentration in the divided area is an index value representing the average density of the substance in the area, and is information useful for identifying, for example, a cell containing a nucleus and a fragment containing no nucleus.
  • a density histogram is created based on the internal luminance value (step S112). This makes it possible to obtain information useful for estimating the internal structure. If a density histogram is created for an entire region, the properties of the entire region can be expressed. In addition, by creating a density histogram for a part of one area, the structure inside the area can be shown in more detail.
  • FIG. 10A and 10B are diagrams showing examples of density histograms.
  • FIG. 10A shows an example of a density histogram for two regions R1 and R2 inside the cell C1 defined by region division.
  • a narrow peak appears in a relatively low density region. It can be assumed that this region is a low material density, for example, a cavity.
  • a broad peak appears in a relatively high concentration region, indicating that the substance is contained at a high concentration and in a localized manner. For example, such a tendency is seen in a region where intracellular organelles are concentrated or a region where movement for cell division is active.
  • the density histogram there is a big difference in the density histogram between the nucleus area in the cell and other areas.
  • the density histogram for example, it becomes possible to provide the user with information that is a clue regarding the internal structure of the cell, such as the peak position, height, and width.
  • Such a technique for obtaining a concentration histogram is also useful for evaluating blastocysts. That is, as shown in FIG. 10B, an inner cell mass I is formed inside a thin trophectoderm T in a blastocyst. In an embryo in a good state, the inner cell mass I is densely composed of a large number of cells, whereas in an inner cell mass of an embryo that is not so, there are cases where the number of cells is low or the density is low. Therefore, when the density histogram is obtained for the inner cell mass I specified by the region division, in the embryo in a good state, the peak spread is small as shown by the solid line in the lower graph of FIG. 10B.
  • the peak spread becomes large as shown by the broken line.
  • information indicating the peak position, height, width, and the like can be used for embryo evaluation.
  • index values for example, by applying these index values to the Gardner classification, it is possible to quantitatively perform classification that conventionally relies on subjective judgment.
  • the structures included in the three-dimensional image may include various types of structures such as blastomeres, cells, and fragments.
  • an appropriate threshold value is set in advance for the index value obtained in this way, it can be implemented as an image determination method for automatically determining what kind of structure the type is.
  • the present invention is not limited to the above-described embodiment, and various modifications other than those described above can be made without departing from the spirit of the present invention.
  • the image processing apparatus 1 of the above embodiment has a function of performing OCT imaging of the sample Sp and a function of generating and outputting an output image from the imaging data.
  • the image processing method of the present invention itself does not have an imaging function, and can also be executed by a computer apparatus that has acquired imaging data obtained by imaging with another apparatus having an imaging function.
  • the present invention may be implemented as a software program for causing a computer apparatus to execute steps S102 to S112 among the processing steps of FIG.
  • Distribution of such a program can be performed in the form of downloading via a telecommunication line such as the Internet, for example. It is also possible to distribute a computer-readable recording medium in which the program is recorded. In addition, by causing an existing OCT imaging apparatus to read this program via an interface, the present invention can be implemented by the imaging apparatus.
  • the process of erasing the second region corresponding to the fragment or the like among the regions divided by the Local Thickness calculation Is running is also possible to erase the first area leaving the second area, for example for the purpose of observing the fragment in detail.
  • an image that individually represents the first area and the second area may be created.
  • processing for dividing a three-dimensional image of an embryo steps S104 to S106), size calculation of the divided region (step S107), and individual cells (Step S108), processing to express the surface of the cell with polar coordinates (step S109), processing to determine the degree of divergence between the cell and the approximate ellipsoid (step S110), processing to determine the concentration distribution in the region (Steps S111 to S112) are executed in this order.
  • processing for dividing a three-dimensional image of an embryo steps S104 to S106
  • size calculation of the divided region step S107
  • individual cells steps S108
  • processing to express the surface of the cell with polar coordinates step S109
  • processing to determine the degree of divergence between the cell and the approximate ellipsoid step S110
  • processing to determine the concentration distribution in the region Steps S111 to S112
  • the execution order is not limited to the above.
  • various quantitative information useful for evaluation work by the user is presented from a three-dimensional image obtained by OCT imaging.
  • the usage form of the information by the user is not limited to those exemplified above, but is arbitrary.
  • the gradation value given to each point in the three-dimensional image by the Local ⁇ Thickness calculation is used as the index value. Also good.
  • the size of an object is expressed by a gradation value.
  • Such a gradation value can be suitably applied as an index value in the above-described processing.
  • the region dividing step may be configured to demarcate a boundary between a plurality of regions and divide the three-dimensional image represented by the original image data into a plurality of demarcated boundaries. By dividing the object in the three-dimensional image at these boundaries, individual cells and other structures of the blastomer can be separated from each other and individually evaluated.
  • the divided area may be configured such that an exposed surface that is not in contact with another area among the surfaces of the area is approximated with an ellipsoid.
  • a cell is indefinite but generally has an ellipsoidal shape. Therefore, by approximating a cell to an ellipsoid, the shape can be expressed as a specific numerical value. In this case, the surface of the cell is deformed in a portion in contact with another cell. Therefore, by approximating such an exposed surface with little deformation to an ellipsoid, errors due to surface deformation can be reduced.
  • an index value indicating the degree of divergence between the surface of the approximate ellipsoid obtained by ellipsoid approximation and the exposed surface may be further calculated. Since it can be said that the approximate ellipsoid shows an ideal shape of the cell, the magnitude of the deviation from the approximate ellipsoid shows the magnitude of the distortion of the cell and the size of the surface irregularities. If there is an index value that quantitatively indicates such a shape-related feature, it will be useful for the user's evaluation work.
  • the coordinates of each position of the exposed surface and the surface of the approximate ellipsoid are represented by polar coordinates with the center of gravity of the approximate ellipsoid as the origin, and the point on the exposed surface on the same radius and the surface of the approximate ellipsoid
  • the index value may be obtained based on the distance to the point.
  • the distance between the cell and the approximate ellipsoid viewed from the center of gravity of the ellipsoid can be used as information that directly represents the degree of divergence between them.
  • a histogram of pixel values of each pixel in the divided area may be calculated. According to such a configuration, it is possible to create information in which the internal structure of the cell is quantified.
  • the image processing method divides a three-dimensional embryo image into a plurality of regions by the above method, and determines a first region having an index value greater than a first threshold as an aggregate of blastomeres.
  • a first threshold as an aggregate of blastomeres.
  • Normal blastomeres in embryos are considered to have a larger structure than other structures. Therefore, it is extremely reasonable to determine the first region estimated to be a relatively large structure as an aggregate of blastomeres.
  • one region further divided from the first region may be determined as one cell.
  • the first region is regarded as a blastomer consisting of a plurality of cells, it is appropriate to consider that each region obtained by further dividing the first region corresponds to one cell.
  • a region having a sphericity greater than the second threshold among regions divided into the second regions whose index values are smaller than the first threshold may be determined as fragments. Fragments are smaller in structure than cells but have a certain sphericity. Accordingly, a region having a relatively high sphericity among the regions divided into the second regions can be regarded as a fragment. In this way, it is possible to distinguish fragments from other structures.
  • the present invention can be applied for the purpose of supporting the work of evaluating the state of the cultured embryo, and can be used, for example, to increase the success probability of in vitro fertilization and artificial insemination in infertility treatment.

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • General Health & Medical Sciences (AREA)
  • Health & Medical Sciences (AREA)
  • Chemical & Material Sciences (AREA)
  • Radiology & Medical Imaging (AREA)
  • Immunology (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Pathology (AREA)
  • Quality & Reliability (AREA)
  • Medical Informatics (AREA)
  • Biochemistry (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Dispersion Chemistry (AREA)
  • Signal Processing (AREA)
  • Geometry (AREA)
  • Investigating Or Analysing Materials By Optical Means (AREA)
  • Image Analysis (AREA)

Abstract

本発明は、ＯＣＴ撮像により得られるデータから定量的な情報を求めるのに適したデータ処理方法を提供する。本発明に係る画像処理方法は、培養された胚を光干渉断層撮像して得られる、胚の三次元像に対応する原画像データを取得する工程（Ｓ１０２）と、原画像データに基づき、三次元像を複数の領域に領域分割する工程（Ｓ１０３～Ｓ１０５）とを備え、領域分割する工程では、三次元像に対しLocal Thickness演算を実行して三次元像に含まれるオブジェクトの大きさを指標する指標値を求め、該指標値が所定の第１閾値より大きい領域と小さい領域とに三次元像を分割し（Ｓ１０４）、分割された領域の各々を、Watershedアルゴリズムにより分割する（Ｓ１０５）。

Description

画像処理方法、プログラムおよび記録媒体
　この発明は、胚を断層撮像することで得られる胚の三次元像に対応するデータから、胚を評価するユーザに有用な情報を提供するための画像処理技術に関するものである。
　例えば不妊治療を目的とした生殖補助医療においては、体外で受精させ一定期間培養した胚（受精卵）を胎内に戻すことが行われる。しかしながら、その成功率は必ずしも高くなく、患者の精神的および経済的な負担も大きい。この問題を解決するために、培養される胚の状態を的確に判断する方法が模索されている。
　従来、胚培養が良好に進行しているか否かの評価については、例えば顕微鏡観察により医師や胚培養士が目視で行うことが一般的である。その判断指標として、例えばVeeck分類やGardner分類などが広く用いられている。ただし、これらは胚の形態学的特徴に対するおおよその判断基準を示したものにすぎない。最終的な評価は、評価者の主観的判断に依存しているのが現状である。このため、客観的かつ定量的な評価を可能とするための技術が求められる。
　例えば特許文献１では、光干渉断層撮像（光コヒーレンストモグラフィ、Optical Coherence Tomography；ＯＣＴ）のような非侵襲の断層撮像技術を応用して胚の品質を評価する技術が提案されている。この技術では、ＯＣＴ撮像により得られる断層画像から胚の三次元像が再構成される。そして、その三次元像に基づき、胚の形態学的基準としての割球の数、その規則性、フラグメンテーション率等の基準を用いて胚の品質が評価される。
特開２０１７－５２１０６７号公報
　上記文献には、ＯＣＴ撮像で得られる三次元像に基づき胚の評価が可能であることが記載されている。また、胚を種々の方向から見た三次元像の例が開示されている。しかしながら、具体的な像の加工や定量的な指標の算出については述べられていない。すなわち、この従来技術は、これまでの目視観察に供されていた顕微鏡像のような二次元像を、単に三次元像に置き換えたものにすぎない。このため、従来の二次元像を用いるケースに比べれば、胚の評価作業をより効果的に支援することができるが、ＯＣＴ撮像結果から得られる情報が十分有効に活用されているとは言えない。
　例えば、上記したVeeck分類やGardner分類などでは、胚に含まれる細胞の大きさや形状がよく揃っていることが良質の胚の条件とされている。このためには、胚を構成する個々の細胞のそれぞれについて、その大きさや形状を個別に評価するための定量的な指標が撮像データから求められれば、ユーザにとって便宜である。しかしながら、上記従来技術は、例えば胚表面の透明体を消去して内部の細胞を露出させた表示態様や、それを回転させ種々の方向から観察できるようにした表示態様が例示されるものの（例えば図５Ａ～図５Ｉ）、上記した評価のための定量化に資するものではない。
　このように、ＯＣＴ撮像による得られる断層画像およびそれから再構成される三次元像は、例えば胚の内部構造等、胚の評価のために有用な多くの情報を提供することができる可能性を有している。しかしながら、そのための具体的なデータ処理方法についてはこれまで確立されるに至っていない。
　この発明は上記課題に鑑みなされたものであり、ＯＣＴ撮像により得られるデータから定量的な情報を求めるのに適したデータ処理方法を提供することで、ユーザによる胚の評価作業を効果的に支援することを目的とする。
　この発明の一の態様は、上記目的を達成するため、培養された胚を光干渉断層撮像して得られる、前記胚の三次元像に対応する原画像データを取得する工程と、前記原画像データに基づき、前記三次元像を複数の領域に領域分割する工程とを備え、前記領域分割する工程では、前記三次元像に対しLocal Thickness演算を実行して前記三次元像に含まれるオブジェクトの大きさを指標する指標値を求め、該指標値が所定の第１閾値より大きい領域と小さい領域とに前記三次元像を分割し、分割された領域の各々を、Watershedアルゴリズムにより分割する画像処理方法である。
　また、この発明の他の一の態様は、上記各工程をコンピュータに実行させるためのプログラムである。また、この発明の他の一の態様は、上記プログラムを記憶したコンピュータ読み取り可能な記録媒体である。
　詳しくは後述するが、胚、特に哺乳類の受精卵は、その外形が概ね球形であり、またその内部も、概ね球形ないし楕円体である細胞が複数集まった構造となっている。このような胚を構成する個々の細胞を個別に評価するためには、各細胞がどのように空間領域を占めているか、言い換えれば三次元空間における細胞間の境界を画定する必要がある。
　このような境界の画定は、ＯＣＴ画像のような三次元像に対する領域分割処理によって実施することが可能である。例えば、公知の画像処理アルゴリズムであるWatershed（分水嶺）アルゴリズムを利用することができる。しかしながら、本願発明者の知見によれば、Watershedアルゴリズムによる領域分割の結果が、三次元像から熟練者が判定する細胞の境界と乖離したものとなるケースが多い。胚の三次元像においては、卵割により形成される細胞の数が多くなるにしたがって細胞の境界が次第に明瞭でなくなってくること、極体や不完全な細胞分裂によって生じた微小な構造体（例えばフラグメント）や画像ノイズなどが含まれること等が、誤分割の原因になっているものと考えられる。
　そこで、本発明では、細胞同士の境界や、細胞と細胞以外の比較的小さな構造体との境界を明確にするために、Local Thickness演算と呼ばれる画像処理の一手法が用いられる（特に三次元像に適用されるとき、3D Local Thickness演算ともいう）。この3D Local Thickness演算では、三次元像中のオブジェクトに内接する球の半径の大きさに基づいて、オブジェクトが階調表現される。例えば、比較的大きな構造を有するオブジェクトは高い階調で表される。一方、比較的小さなオブジェクトはより低い階調で表される。このようにオブジェクトの大きさに依存する値（この場合は階調値）を指標値として利用すれば、細胞などに対応する比較的大きなオブジェクトとフラグメントなどに対応する微小なオブジェクトとの区別が可能となる。また細胞同士の境界が明瞭となる。このため、こうして分割される領域に対しさらにWatershedアルゴリズムによる領域分割を実行すれば、実際の細胞等の境界により合致した分割結果を得ることができる。
　このため、上記の画像処理方法による分割された領域のうち、比較的大きなものを個々の細胞として個別に評価することで、その形状や大きさ等を定量的に示す種々の指標値を求めることが可能になる。例えば、Local Thickness演算により得られる指標値が所定の第１閾値より大きいか小さいかで、細胞に対応する領域とそれ以外の領域とを区別することが可能である。
　また、この発明の他の一の態様は、上記した画像処理方法により、胚の三次元像を複数の領域に分割し、前記複数の領域のうち前記指標値が前記第１閾値より大きい第１領域を割球の集合体と判定する画像処理方法である。これにより、胚において比較的大きな構造体である細胞と、極体やフラグメント等の比較的小さな構造体とを区別することが可能となる。また、それぞれの三次元形状等の特徴を定量的に求めて、胚の評価に供することが可能となる。
　上記のように、本発明によれば、培養された胚をＯＣＴ撮像することにより得られるデータから、胚を構成する個々の細胞の評価に必要な定量的な情報を求めることが可能になる。
　この発明の前記ならびにその他の目的と新規な特徴は、添付図面を参照しながら次の詳細な説明を読めば、より完全に明らかとなるであろう。ただし、図面は専ら解説のためのものであって、この発明の範囲を限定するものではない。
断層撮像を行うための画像処理装置の構成例を示す原理図である。 この画像処理装置における撮像原理を説明する第１の図である。 この画像処理装置における撮像原理を説明する第２の図である。 ＯＣＴ装置の具体的構成例を示す第１の図である。 ＯＣＴ装置の具体的構成例を示す第２の図である。 処理対象となる胚の代表的な構造を例示する第１の図である。 処理対象となる胚の代表的な構造を例示する第２の図である。 処理対象となる胚の代表的な構造を例示する第３の図である。 処理対象となる胚の代表的な構造を例示する第４の図である。 この実施形態の画像処理方法を示すフローチャートである。 大小の構造体を分離する処理の原理を示す第１の図である。 大小の構造体を分離する処理の原理を示す第２の図である。 大小の構造体を分離する処理の原理を示す第３の図である。 大小の構造体を分離する処理の原理を示す第４の図である。 胚盤胞への適用例を示す第１の図である。 胚盤胞への適用例を示す第２の図である。 領域分割による個々の細胞の分離を模式的に示す第１の図である。 領域分割による個々の細胞の分離を模式的に示す第２の図である。 領域分割による個々の細胞の分離を模式的に示す第３の図である。 細胞の楕円体近似を模式的に示す第１の図である。 細胞の楕円体近似を模式的に示す第２の図である。 細胞の楕円体近似を模式的に示す第３の図である。 細胞の楕円体近似を模式的に示す第４の図である。 濃度ヒストグラムの例を示す第１の図である。 濃度ヒストグラムの例を示す第２の図である。
　以下、本発明に係る画像処理方法の具体的な実施形態について説明する。本実施形態は、本発明に係る画像処理方法を具現化したものである。この実施形態では、培養容器に担持された培地（培養液）中で培養された胚が、光干渉断層撮像（Optical Coherence Tomography；ＯＣＴ）技術によって断層撮像され、これにより得られたデータに基づき、胚を構成する細胞の定量的情報が出力される。その目的は、ユーザ（具体的には医師または胚培養士）による胚の評価作業を支援することである。例えば不妊治療を目的とする受精卵の培養において、培養が良好に進行しているか否かを判断するための知見を得ることを目的として、本実施形態の画像処理方法を適用することができる。
　最初に、断層画像を撮像するための画像処理装置の構成と、当該画像処理装置による撮像原理とについて説明する。その後に、撮像により取得されたデータを用いた画像処理について説明する。
　図１は断層撮像を行うための画像処理装置の構成例を示す原理図である。この画像処理装置１は、培地中で培養された胚を被撮像物として断層撮像し、得られた断層画像を画像処理して、被撮像物の三次元（３Ｄ）像を作成する。以下の各図における方向を統一的に示すために、図１に示すようにｘｙｚ直交座標軸を設定する。ここでｘｙ平面が水平面を表す。また、ｚ軸が鉛直軸を表し、より詳しくは（－ｚ）方向が鉛直下向き方向を表している。
　画像処理装置１は保持部１０を備えている。保持部１０は、培養容器１１をその開口面を上向きにして略水平姿勢に保持する。培養容器（以下、単に「容器」という）１１は、ガラス製または樹脂製の透明な底部を持つ、ディッシュやプレートと呼ばれるものである。容器１１には予め適宜の培地Ｍが所定量注入されており、培地中では容器１１の底部１１１に試料Ｓｐ（ここでは胚）が培養されている。図１では１つの試料Ｓｐのみが記載されているが、１つの容器１１内で複数の試料Ｓｐが培養されていてもよい。
　保持部１０により保持された容器１１の下方に、撮像ユニット２０が配置される。撮像ユニット２０には、被撮像物の断層画像を非接触、非破壊（非侵襲）で撮像することが可能な、光干渉断層撮像（Optical Coherence Tomography；ＯＣＴ）装置が用いられる。詳しくは後述するが、ＯＣＴ装置である撮像ユニット２０は、被撮像物への照明光を発生する光源２１と、ビームスプリッタ２２と、物体光学系２３と、参照ミラー２４と、分光器２５と、光検出器２６とを備えている。
　また、画像処理装置１はさらに、装置の動作を制御する制御ユニット３０と、撮像ユニット２０の可動機構を制御する駆動制御部４０とを備えている。制御ユニット３０は、ＣＰＵ（Central Processing Unit）３１、Ａ／Ｄコンバータ３２、信号処理部３３、３Ｄ復元部３４、インターフェース（ＩＦ）部３５、画像メモリ３６およびメモリ３７を備えている。
　ＣＰＵ３１は、所定の制御プログラムを実行することで装置全体の動作を司る。ＣＰＵ３１が実行する制御プログラムや処理中に生成したデータは、メモリ３７に保存される。Ａ／Ｄコンバータ３２は、撮像ユニット２０の光検出器２６から受光光量に応じて出力される信号をデジタルデータに変換する。信号処理部３３は、Ａ／Ｄコンバータ３２から出力されるデジタルデータに基づき、後述する信号処理を行って、被撮像物の断層画像を作成する。３Ｄ復元部３４は、撮像された複数の断層画像の画像データに基づいて、撮像された胚の三次元像（３Ｄ像）を作成する機能を有する。信号処理部３３により作成された断層画像の画像データおよび３Ｄ復元部３４により作成された三次元像の画像データは、画像メモリ３６により適宜記憶保存される。
　インターフェース部３５は、画像処理装置１と外部との通信を担う。具体的には、インターフェース部３５は、外部機器と通信を行うための通信機能と、ユーザからの操作入力を受け付け、また各種の情報をユーザに報知するためのユーザインターフェース機能とを有する。この目的のために、インターフェース部３５には、入力デバイス３５１と表示部３５２とが接続されている。入力デバイス３５１は、装置の機能選択や動作条件設定などに関する操作入力を受け付け可能な例えばキーボード、マウス、タッチパネルなどで構成される。また、表示部３５２は、信号処理部３３により作成された断層画像や３Ｄ復元部３４により作成された立体像など各種の処理結果を表示する表示デバイス、例えば液晶ディスプレイを有する。
　また、ＣＰＵ３１は駆動制御部４０に制御指令を与える。これに応じて駆動制御部４０は、撮像ユニット２０の可動機構に所定の動作を行わせる。次に説明するように、駆動制御部４０により実行される撮像ユニット２０の走査移動と、光検出器２６による受光光量の検出との組み合わせにより、被撮像物である試料Ｓｐ（具体的には胚）の断層画像が取得される。
　図２Ａおよび図２Ｂはこの画像処理装置における撮像原理を説明する図である。より具体的には、図２Ａは撮像ユニット２０における光路を示す図であり、図２Ｂは試料Ｓｐの断層撮像の様子を模式的に示す図である。前記したように、撮像ユニット２０は光干渉断層撮像（ＯＣＴ）装置として機能するものである。なお、ここでは説明のために試料Ｓｐを略球形として表示している。しかしながら、実際の撮像において試料Ｓｐの形状は特に限定されない。
　撮像ユニット２０では、例えば発光ダイオードまたはスーパールミネッセントダイオード（ＳＬＤ）などの発光素子を有する光源２１から、広帯域の波長成分を含む低コヒーレンス光ビームＬ１が出射される。光ビームＬ１はビームスプリッタ２２に入射して分岐する。破線矢印で示すように一部の光Ｌ２が容器１１に向かい、一点鎖線矢印で示すように一部の光Ｌ３が参照ミラー２４に向かう。
　容器１１に向かった光Ｌ２は、物体光学系２３を経て容器１１に入射する。より具体的には、ビームスプリッタ２２から出射される光Ｌ２は、物体光学系２３を介して容器底部１１１に入射する。物体光学系２３は、ビームスプリッタ２２から容器１１に向かう光Ｌ２を容器１１内の試料Ｓｐに収束させる機能と、試料Ｓｐから出射される反射光を集光してビームスプリッタ２２に向かわせる機能とを有する。図では物体光学系２３は単一の対物レンズにより代表的に表されているが、複数の光学素子が組み合わされたものであってもよい。
　物体光学系２３は、駆動制御部４０に設けられた焦点調整機構４１により、ｚ方向に移動可能となっている。これにより、被撮像物に対する物体光学系２３の焦点位置がｚ方向に変更可能である。物体光学系２３の光軸は鉛直方向と平行であり、したがって平面状の容器底部１１１に垂直である。また、物体光学系２３への照明光の入射方向は光軸と平行であり、照明光の光中心が光軸と一致するように、物体光学系２３の配置が定められている。
　試料Ｓｐが光Ｌ２に対する透過性を有するものでなければ、容器底部１１１を介して入射した光Ｌ２は試料Ｓｐの表面で反射される。一方、試料Ｓｐが光Ｌ２に対してある程度の透過性を有するものである場合、光Ｌ２は試料Ｓｐ内まで進入してその内部の構造物により反射される。光Ｌ２として例えば近赤外線を用いれば、入射光を試料Ｓｐ内部まで到達させることが可能である。試料Ｓｐからの反射光は、散乱光として種々の方向に放射される。そのうち物体光学系２３の集光範囲内に放射された光Ｌ４が、物体光学系２３で集光されてビームスプリッタ２２へ送られる。
　参照ミラー２４は、その反射面を光Ｌ３の入射方向に対し垂直姿勢に支持されている。また、駆動制御部４０に設けられたミラー駆動機構４２により、光Ｌ３の入射方向に沿った方向（図ではｙ方向）に移動可能となっている。参照ミラー２４に入射した光Ｌ３は反射面で反射されて、入射光路を逆向きに辿るように進む光Ｌ５としてビームスプリッタ２２に向かう。この光Ｌ５が参照光となる。ミラー駆動機構４２により参照ミラー２４の位置が変更されることにより、参照光の光路長が変化する。
　試料Ｓｐの表面もしくは内部の反射面で反射された反射光Ｌ４と、参照ミラー２４で反射された参照光Ｌ５とは、ビームスプリッタ２２を介して光検出器２６に入射する。このとき、反射光Ｌ４と参照光Ｌ５との間で位相差に起因する干渉が生じるが、干渉光の分光スペクトルは反射面の深さにより異なる。つまり、干渉光の分光スペクトルは、被撮像物の深さ方向の情報を有している。したがって、干渉光を波長ごとに分光して光量を検出し、検出された干渉信号をフーリエ変換することにより、被撮像物の深さ方向における反射光強度分布を求めることができる。このような原理に基づくＯＣＴ撮像技術は、フーリエドメイン（Fourier Domain）ＯＣＴ（ＦＤ－ＯＣＴ）と称される。
　この実施形態の撮像ユニット２０では、ビームスプリッタ２２から光検出器２６に至る干渉光の光路上に、分光器２５が設けられている。分光器２５としては、例えばプリズムを利用したもの、回折格子を利用したもの等を用いることができる。干渉光は、分光器２５により波長成分ごとに分光されて光検出器２６に受光される。
　光検出器２６が検出した干渉光に応じて光検出器２６から出力される干渉信号をフーリエ変換することで、試料Ｓｐのうち、光ビームＬ２の入射位置における深さ方向、つまりｚ方向の反射光強度分布が求められる。容器１１に入射する光ビームＬ２をｘ方向に走査することで、ｘｚ平面と平行な平面における反射光強度分布が求められる。その結果から、当該平面を断面とする試料Ｓｐの断層画像を作成することができる。以下、本明細書では、ｘ方向へのビーム走査によってｘｚ平面と平行な断面における１つの断層画像Ｉｔを取得する一連の動作を、１回の撮像と称することとする。
　また、ｙ方向におけるビーム入射位置を多段階に変更しながら、その都度断層画像の撮像が行われる。こうすることで、図２Ｂに示すように、試料Ｓｐをｘｚ平面と平行な断面で断層撮像した多数の断層画像Ｉｔを得ることができる。ｙ方向の走査ピッチを小さくすれば、試料Ｓｐの立体構造を把握するのに十分な分解能の画像データを得ることができる。ｘ方向およびｙ方向へのビーム走査は、種々の方法により実現可能である。例えば、図示しないガルバノミラー等の光路を変化させる光学部品を用いてビーム入射位置をｘｙ方向に変化させる方法を適用可能である。また、試料Ｓｐを担持する容器１１と撮像ユニット２０とのいずれかをｘｙ方向に移動させてこれらの相対位置を変化させる方法も適用可能である。
　なお、上記の原理説明では、撮像ユニット２０において光源２１からの光を照明光と参照光とに分岐させる分波機能、および信号光と参照光とを合成して干渉光を生じさせる機能が、ビームスプリッタ２２により実現されている。一方、近年では、ＯＣＴ装置においてこのような分波・合波機能を担うものとして、以下に例示するような光ファイバカプラが用いられる場合がある。
　図３Ａおよび図３ＢはＯＣＴ装置の具体的構成例を示す図である。なお、理解を容易にするために、以下の説明では、上記した原理図の構成と同一のまたは相当する構成に同一符号を付すものとする。その構造および機能は、特に説明のない限り上記原理図のものと基本的に同じであり、詳しい説明は省略する。また、光ファイバカプラによる干渉光を検出するＯＣＴ撮像原理も基本的に上記と同じであるので、詳しい説明を省略する。
　図３Ａに示す構成例では、撮像ユニット２０ａは、ビームスプリッタ２２に代わる分波・合波器として光ファイバカプラ２２０を備えている。光ファイバカプラ２２０を構成する１つの光ファイバ２２１は光源２１に接続されており、光源２１から出射される低コヒーレンス光は、光ファイバカプラ２２０により２つの光ファイバ２２２，２２３への光に分岐される。光ファイバ２２２は物体系光路を構成する。より具体的には、光ファイバ２２２の端部から出射される光は、コリメータレンズ２２３を介して物体光学系２３に入射する。被撮像物からの反射光（信号光）は、物体光学系２３、コリメータレンズ２２３を介して光ファイバ２２２に入射する。
　他の光ファイバ２２４は参照系光路を構成する。より具体的には、光ファイバ２２４の端部から出射される光は、コリメータレンズ２２５を介して参照ミラー２４に入射する。参照ミラー２４からの反射光（参照光）は、コリメータレンズ２２５を介して光ファイバ２２４に入射する。光ファイバ２２２を伝搬する信号光と光ファイバ２２４を伝搬する参照光とが、光ファイバカプラ２２０において干渉する。干渉光は光ファイバ２２６および分光器２５を介して光検出器２６に入射する。光検出器２６により受光された干渉光から被撮像物における反射光の強度分布が求められることは上記原理通りである。
　図３Ｂに示す例でも、撮像ユニット２０ｂに光ファイバカプラ２２０が設けられる。ただし光ファイバ２２４は使用されず、光ファイバ２２２から出射される光の光路に対してコリメータレンズ２２３およびビームスプリッタ２２７が設けられる。そして、前述の原理通り、ビームスプリッタ２２７により分岐される２つの光路に、それぞれ物体光学系２３、参照ミラー２４が配置される。このような構成では、ビームスプリッタ２２７により信号光と参照光とが合成される。それにより生じた干渉光が、光ファイバ２２２，２２６を通って光検出器２６へ導かれる。
　これらの例では、図２Ａの原理図では空間中を進行する各光の光路の一部が光ファイバに置き換えられているが、動作原理は同じである。これらの例においても、焦点調整機構４１が物体光学系２３を容器１１に対し接近・離間方向に移動させることにより、被撮像物に対する物体光学系２３の焦点深さを調整することが可能である。また、ミラー駆動機構４２が参照ミラー２４を光の入射方向に沿って移動させることにより、参照光の光路長を変更可能である。
　以下、この画像処理装置１を用いた画像処理方法について詳説する。画像処理装置１の構成は、上記したビームスプリッタを用いるもの、光ファイバカプラを用いるもののいずれであっても適用可能である。また、断層画像を撮像するための撮像装置としては上記したＦＤ－ＯＣＴ撮像装置に限定されない。例えばタイムドメインＯＣＴ（ＴＤ－ＯＣＴ）撮像装置など他の撮像原理に基づくものも適用可能である。以下の画像処理方法では、試料Ｓｐおよびその周囲を含む三次元空間内の各位置に入射した照明光の反射光が信号光とされる。そして、三次元空間内の位置と、当該位置からの信号光の光量に対応する信号の強度とを対応付けたデータが処理に使用される。したがって、このようなデータを取得することのできる撮像方法であればよい。
　図４Ａないし図４Ｄは処理対象となる胚の代表的な構造を例示する模式図である。画像処理装置１が、例えば不妊治療を目的として利用されるとき、体外受精されたヒト受精卵の、培養の初期段階のものが画像処理装置１の撮像対象物となる。既に知られているように、卵子が受精すると卵割が開始され、桑実胚と呼ばれる状態を経て胚盤胞が形成される。本実施形態の画像処理方法は、例えば受精直後から胚盤胞期における胚を観察するのに好適な情報を、ユーザに提供するものである。
　図４Ａは初期段階（例えば４細胞期から桑実胚期）の胚の構造を模式的に示すものである。胚Ｅ１の外形は概ね球形である。その表面は透明帯と呼ばれるゼリー状の糖タンパク質の層Ｚで覆われている。その内部には、受精卵が細胞分裂することで生じる複数の細胞Ｃが含まれる。培養が良好に進行する状態では、図４Ａに示すように、透明帯Ｚの内部は大きさの均等な複数の細胞Ｃにより占められており、卵割が進むにつれてその数は増えてゆく。Veeck分類によると、それらの細胞は揃っていることが望ましいとされている。一方、培養の状態が良好でない場合、図４Ｂに示す胚Ｅ２のように、細胞Ｃの大きさが不揃いであったり、卵割により生成される細胞Ｃ以外にフラグメントと呼ばれる微小な断片が生じたりすることがある。なお、図４Ｂでは細胞Ｃの大きさのばらつきとフラグメントＦとが共に発生した状態を示しているが、これらの一方のみが出現することもある。
　いずれにしても、胚の状態としては図４Ａに示すものよりも劣っていると言える。したがって、桑実胚までの初期胚においては、卵割により形成される各細胞Ｃの大きさが揃っているか、各細胞Ｃの表面が滑らかであるか、フラグメントＦが生じているか等が、培養状態を評価するための指標となり得る。
　図４Ｃおよび図４Ｄは胚盤胞期における胚の内部構造を模式的に示すものである。このうち図４Ｃは状態が良好な胚Ｅ３を示し、図４Ｄはこれより状態の劣る胚Ｅ４を示している。胚盤胞期においては、図４Ｃに示すように、卵割の進んだ細胞Ｃが胚の表面に薄い層として並び栄養外胚葉Ｔを形成する。栄養外胚葉Ｔで囲まれる内部空間は、胞胚腔Ｂとよばれる空洞が生じる。また、内部空間の一部に、多数の細胞Ｃが密集した内細胞塊Ｉが形成される。
　状態のよい胚Ｅ３では、栄養外胚葉Ｔが多数の細胞Ｃにより密に形成され、また内細胞塊Ｉも比較的多くの細胞Ｃにより密に構成される。一方、これより状態の劣る胚Ｅ４では、栄養外胚葉Ｔがより少ない細胞Ｃで疎に形成されることや、内細胞塊Ｉが少ない細胞で構成され小さいものとなること等が生じ得る。したがって、栄養外胚葉Ｔの厚さや密度、内細胞塊Ｉの大きさや密度により、胚の培養状態を指標することができる。
　図５はこの実施形態の画像処理方法を示すフローチャートである。なお、以下に示す一連の処理のうち、ステップＳ１０１は、画像処理装置１とは別に準備された適宜のインキュベータ内で実施される。また、ステップＳ１０２以降の処理については、画像処理装置１に設けられたＣＰＵ３１が、メモリ３７に予め記録された制御プログラムを実行することにより実現されるが、これに限定されない。
　例えば、図５のステップＳ１０３以降については、一般的な演算機能および画像出力機能を有するコンピュータ装置によって実行することが可能である。これらの処理については、画像処理装置１とは別のコンピュータ装置が画像処理装置１からＯＣＴ撮像データを受け取って処理する構成であってもよい。このようにすれば、画像処理装置１のＣＰＵ３１は撮像に特化された処理のみを実行すればよいこととなる。これにより、ＣＰＵ３１の処理負荷が軽減される。
　以下、具体的な処理の内容を説明する。最初に、評価の対象となる胚が適宜の培養環境下で培養される（ステップＳ１０１）。培養された胚は、適宜のタイミングで、画像処理装置１によりＯＣＴ撮像される（ステップＳ１０２）。撮像により得られたデータ、具体的には、胚およびその周囲を含む三次元空間の各位置からの反射光強度を表すデータが、原信号データとしてメモリ３７または画像メモリ３６に保存される。この原信号データを用い、三次元空間における座標位置ごとの画素データ（ボクセルデータ）を作成することにより、特許文献１に示されたような胚の三次元像を構成し、画像として出力することが可能である。作成されたボクセルデータは、原画像データとして画像メモリ３６に保存される。
　後の処理のために三次元像が二値化される（ステップＳ１０３）。例えば、ボクセルデータにより表される各画素のうち、画素値（輝度値）が予め定められた閾値を超える画素が「１」、画素値が閾値より小さい画素が「０」で表される。これにより、三次元像は各位置の輝度に応じて二値化される。原画像データを可視化したＯＣＴ画像では、撮像範囲に含まれる三次元空間のうち反射光強度の高い、つまり光を反射する物質の密度が高い領域が高輝度で表される。一方、物質の密度が低く反射光強度の低い領域は低輝度で表される。三次元像を二値化することで、物質の密度が一定の値以上である領域とそれ以外の領域とを区分することができる。
　前記したように、胚の三次元像には、卵割により形成された割球の他に、透明帯、栄養外胚葉、フラグメント等の構造体が含まれ得る。胚を詳細に評価する上で、これらを識別する必要がある。この目的のために、まず三次元像に含まれる構造体間の境界を画定させる。このような境界の画定には、三次元空間における領域分割を実行する画像処理アルゴリズム、例えば公知のWatershed（分水嶺）アルゴリズムを利用することが可能である。Watershedアルゴリズムでは、二値化された三次元像の凸部と凹部との境界が領域の境界として特定される。
　しかしながら、ＯＣＴ撮像により得られた胚の三次元像においては、Watershedアルゴリズムによる分割結果が実際の胚構造と合致しないケースが多い。本願発明者の知見によれば、これは比較的大きな構造を有する割球の細胞間での境界の不明瞭さや、割球の周囲に分布する、透明帯、栄養外胚葉、フラグメント等のより小さな構造を有する構造体や、画像ノイズの影響と考えられる。そこで、Watershedアルゴリズムによる領域分割に先立って、三次元像における微小な構造体をより大きな構造体から分離するための処理を実行する。具体的には、以下の原理に基づく処理が実行される。
　ここでは、図４Ｂに示した胚Ｅ２を対象とする場合を例として説明する。図４Ｂに示すように、胚Ｅ２の像は、卵割により形成される複数の細胞Ｃ、その周囲に分布するフラグメントＦおよびこれらの取り囲む透明帯Ｚを含む。
　以下、二値化された三次元像において、構造体によって占められた空間領域を「構造体内部空間」と称する。また、三次元像における構造体の表面、つまり上記した構造体内部空間とその周囲の空間との境界である曲面を、「構造体表面」と称する。また、以下では、理解を容易にするために、三次元像を一の断面で切断した二次元画像を用いて処理の説明を行う。ただし実際の画像処理は、ボクセルデータにより表される三次元画像に対して行われる。
　図６Ａないし図６Ｄは大きな構造体と小さな構造体とを分離する処理の原理を示す図である。図６Ａは、図４Ｂに示す胚Ｅ２の二値化画像に相当する二次元画像である。この処理の基本原理では、図６Ｂに示すように、二値化画像を内接球（二次元画像では内接円）の集合体としてモデル化する。そして、適宜の閾値（本発明の「第１閾値」）より大きい半径の内接球で表される領域（以下、「第１領域」という）と、閾値より小さい半径の内接球で表される領域（以下、「第２領域」という）とに区分する。なお、内接球の半径が閾値と等しいケースを第１領域、第２領域のいずれに含めるかは任意である。このようなモデルでは、細胞Ｃおよびその凝集塊のように比較的大きな構造体は、比較的大きな半径の内接球の集合体として表される。一方、透明帯ＺやフラグメントＦ等の小さな構造体は、より小さな半径の内接球の集合体として表される。したがって、内接球の半径の大きさによって構造体を区分することで上記目的が達成される。例えば、図６Ｂに示す閾値Ｒｔｈ以下の半径の内接球で表される構造体を消去すれば、図６Ｃに示すように、大きな構造を有する細胞の集合体に対応する領域のみを抽出することができる。
　原理的には、次のような演算により構造体のモデル化が可能である。構造体内部空間の各点のうち、当該点を中心として構造体表面に内接する内接球の半径が最大となる、言い換えれば構造体表面までの距離が最大である点を見出す。中心の位置および半径の大きさにより、このときの内接球を特定することができる。このような点は必ずしも１つとは限らない。次に、既に特定された内接球に含まれない構造体内部空間の各点のうち、当該点を中心として構造体表面に内接する内接球の半径が最大となる点を見出す。このような処理を構造体内部空間の全ての点が少なくとも１つの内接球に含まれるまで繰り返すことで、構造体が内接球の集合体として表されることになる。
　このような演算を収束させるには長い時間が必要である。より短時間で同様の結果を得られる演算方法として、例えばLocal Thickness演算アルゴリズムが提唱されており（例えば、“Computing Local Thickness of 3D Structures with ImageJ”, Robert P. Dougherty and Karl-Heinz Kunzelmann, Microscopy & Microanalysis 2007 Meeting, August 5-9, Broward County Convention Center, Ft. Lauderdale, Florida）、本実施形態にもこれを適用することが可能である。この演算では、画像中のオブジェクトの各点について、当該点を含み、かつ構造体表面に内接する球のうち、半径が最大であるものが探索される。そして、内接球の半径が大きいほど当該点に対して高いスコアが与えられる。演算結果は、スコアの高い点ほど高い階調値（輝度値）となるような表示態様で出力される。
　このように、Local Thickness演算では、三次元像中のオブジェクトが、当該オブジェクトに内接する球の半径の大きさに応じた階調値で表現される。そのため、比較的大きなオブジェクトは比較的高い階調値（輝度値）に、比較的小さなオブジェクトは比較的低い階調値（輝度値）に階調変換される。この階調値（輝度値）を、オブジェクトを区別するための指標値として利用すれば、細胞などの比較的大きなオブジェクトとフラグメントなどの微小なオブジェクトとの区別が可能になる。これに加えて、細胞同士の境界が明瞭となり、領域間の適切な分割が可能となる。
　例えば図６Ｄに示すように、内接球の半径が最大値Ｒｍａｘである点に対して最大階調値Ｇｍａｘが割り当てられ、内接球の半径が小さいほど階調値も小さくなるように、構造体内部空間の各点に階調値が割り当てられる。Local Thickness演算によって求められた各点の階調値に対し、内接球半径が閾値Ｒｔｈであるときの階調値に対応する閾値Ｇｔｈを設定した場合を考える。そうすると、閾値Ｇｔｈより大きい階調値が割り当てられた点は、比較的大きな構造に含まれる第１領域の内部に存在するとみなすことができる。一方、閾値Ｇｔｈより小さい階調値が割り当てられた点は、比較的小さな構造に含まれる第２領域の内部に存在するとみなすことができる。
　前記したように、卵割により形成される細胞Ｃおよびそれらの凝集体、または、割球およびその集合体は、比較的大きな構造を有する。すなわち、上記処理で第１領域に区分された領域は、細胞または割球やその集合体が占める領域であると判定することができる。一方、第２領域に区分された領域は、細胞凝集塊とは異なる構造体、すなわち透明帯ＺあるいはフラグメントＦが占める領域であると判定することができる。
　透明帯Ｚは胚の表面を覆う連続した薄い層である一方、フラグメントＦは割球やその集合体から遊離した小さな塊である。このことから、第２領域に区分される領域のうち、例えば球形度が適宜の閾値（本発明の「第２閾値」）より大きい領域をフラグメントＦと判定することができる。また、球形度が小さい領域のうち一定以上の体積を有するものが透明帯Ｚ，それ以外を画像ノイズによる擬似的な構造体であると判定することができる。このように、微小な構造体を割球やその集合体から分離することで、その三次元像から形状やサイズ等の定量的情報を求めることが可能となる。
　これにより、それらの構造体についての評価に有用な情報をユーザに提供することが可能となる。例えば、フラグメントの数、大きさ等は、Veeck分類に基づく胚の分類作業を効果的に支援するための情報となる。
　図７Ａおよび図７Ｂは胚盤胞への適用例を示す図である。図４Ｃにも示したように、胚盤胞Ｅ３では、薄い細胞の層である栄養外胚葉Ｔの内部に比較的大きな構造体である内細胞塊Ｉが存在している。このため、図７Ａに示すように、栄養外胚葉Ｔは半径の小さな内接球の集合体として表される一方、内細胞塊Ｉはより大きな内接球の集合体として表される。したがって、図７Ｂに実線で示すように、こうしてモデル化された二値化画像から半径の大きな内接球で表される領域を切り出せば、内細胞塊Ｉに対応する領域を抽出することができる。一方、半径の小さな内接球で表される領域を切り出せば、栄養外胚葉Ｔに対応する領域を抽出することができる。
　図５に戻って、本実施形態の画像処理方法の説明を続ける。二値化画像に対し上記原理に基づく演算処理、例えばLocal Thickness演算を行う（ステップＳ１０４）。こうすることで、二値化画像中のオブジェクトは、割球の集合体に対応する領域（第１領域）とそれ以外の領域（第２領域）とに区分される。これにより、上記したようにそれぞれの領域を個別に定量的情報で表現することが可能となる。細胞の集合体に含まれる割球を個々に評価する目的においては、これらのうち第１領域が処理対象となる。
　すなわち、内接球の半径が閾値Ｒｔｈより小さい第２領域を二値化画像から除去することで（ステップＳ１０５）、割球の集合体に対応する第１領域のみが残る。この第１領域に対しWatershedアルゴリズムによる領域分割処理が実行される（ステップＳ１０６）。これにより、割球の集合体を構成する個々の細胞が分離される。誤差の原因となるフラグメントＦが除去されているため、Watershedアルゴリズムによって個々の細胞間の境界を精度よく画定することが可能である。なお、Local Thickness演算アルゴリズムおよびWatershedアルゴリズムは、いずれも汎用の画像処理ソフトウェアのライブラリとして既に実用化されたものである。
　こうして三次元空間において各細胞が占める領域が特定されれば、個々の細胞のサイズ、例えば体積や表面積を容易に求めることができる（ステップＳ１０７）。例えば、割球の集合体を構成する各細胞の体積のばらつきは、Veeck分類に基づき胚の状態の良否を判定するための指標として利用することができる。
　また、各細胞をより詳細に評価するために、細胞の形状に対応する指標値を導入することが可能である。個々の細胞は、独立した状態であれば概ね球形を示すと考えられる。しかしながら、複数の細胞が集まった割球の集合体では、特に他の細胞と接する部分において変形が生じる。また、胚培養の現場では細胞の表面が滑らかである方がよいとの知見もあり、このような形状に関する特徴を定量的に表現することが望まれる。
　そこで、この実施形態では、１つの細胞として特定された個々の領域が楕円体に近似される（ステップＳ１０８）。これにより、各細胞の表面と近似楕円体の表面との乖離の程度を定量化して表すことで、上記要求に応えられるようにしている。具体的には次の通りである。なお、ここでは互いに接した状態の２つの細胞について考える。しかしながら、より多くの細胞が互いに接している割球の集合体についても同様に考えることが可能である。
　図８Ａないし図８Ｃは領域分割による個々の細胞の分離を模式的に示す図である。また、図９Ａないし図９Ｄは細胞の楕円体近似を模式的に示す図である。図８Ａに示すように、領域分割前の状態の二値化画像では２つの細胞Ｃ１、Ｃ２が一体のオブジェクトとなっている。Watershedアルゴリズムによる領域分割処理を実行することで、図８Ｂおよび図８Ｃに示すように、２つの細胞Ｃ１、Ｃ２の境界が画定される。これにより両者は三次元空間において独立したオブジェクトとして分離される。図８Ｃにおいて、実線は他の細胞と接することなく露出した細胞の表面（以下「露出表面」という）の輪郭を示し、点線は領域分割処理によって画定された細胞間の境界を示す。
　こうして分割された細胞は全体として球形または楕円体をなしていると考えられる。ただし、他の細胞と接していた部分では大きく変形している場合がある。そこで、分割された領域の表面のうち、他の細胞とは接していなかった露出表面の情報のみを使って楕円体近似を行う。こうすることにより、細胞の変形に起因する誤差を抑えることができる。
　具体的には、まず露出表面の抽出を行う。すなわち、分割された細胞の領域の二値化画像を１画素分膨張処理し、分割前の二値化画像における割球の集合体の周囲領域との間で、画素ごとの論理積に相当する画像を作成する。このようにすると、１画素分の厚さを有する露出表面の像が得られる。分割処理により作成された境界面は元の割球の集合体の内部に含まれていたものであるため、論理積の結果には現れない。
　こうして得られた露出表面の形状から近似楕円体を求める。具体的には、露出表面に含まれる各点の空間座標（ｘ，ｙ，ｚ）と、楕円体の一般式である下式：
　　Ａｘ2＋Ｂｙ2＋Ｃｚ2＋Ｄｘｙ＋Ｅｘｚ＋Ｆｙｚ＋Ｇｘ＋Ｈｙ＋Ｉｚ＋Ｊ＝０
から、例えば最小二乗法によって各定数Ａ～Ｊを求める。こうすることにより、当該領域を近似する楕円体を特定することができる。図９Ａに破線で示すように、細胞Ｃ１に対応する近似楕円体Ｅａ１は、実線で示される細胞Ｃ１の露出表面の形状を近似したものであり、点線で示される、領域分割によって作成された境界とは必ずしも一致しない。すなわち、近似楕円体Ｅａ１は、他の細胞との接触による細胞Ｃ１の変形による影響を受けない。細胞Ｃ２に対応する近似楕円体Ｅａ２についても同様である。
　このように、個々の細胞を楕円体によって近似的に表すことにより、その形状を表すパラメータによって、細胞の形状を定量的に表すことができる。パラメータとしては、例えば長軸および短軸の長さ、体積、扁平率、長軸の方向等を用いることができる。これにより、個別の細胞の形状の評価および複数の細胞間での形状の比較等を容易に、かつ定量的に行うことが可能となる。
　細胞Ｃ１の露出表面においても、凹凸のある表面と近似楕円体Ｅａ１の表面とは完全には一致しない。ただし、細胞の状態が良好で表面が滑らかであれば、近似楕円体の表面と露出表面との乖離は小さいと考えられる。すなわち、細胞の露出表面と近似楕円体の表面との乖離の程度を表す数値は、細胞の状態を定量的に指標する値となり得る。
　このような指標値を導出するために、細胞表面の各点の座標を極座標により表すこととする。撮像対象物に対しｘｙ方向へのビーム走査で取得される原画像データ（ボクセルデータ）では、ｘｙｚ座標系で三次元空間内の各位置が表されている。しかしながら、細胞の概略形状は楕円体で表すことができるから、細胞の形状をより端的に表現するためには、細胞の中心を原点とする球座標で各位置を表す方が好ましい。そこで、図９Ｂに示すように、ｘｙｚ直交座標系から、動径ｒおよび２つの偏角θ、φを座標変数とするｒθφ球座標系への座標変換を行う。
　よく知られているように、直交座標系における点Ｐの座標（ｘ，ｙ，ｚ）と球座標系における点Ｐの座標（ｒ，θ，φ）との間には、原点Ｏが共通であれば下式：
　　ｘ＝ｒ・ｓｉｎθ・ｃｏｓφ
　　ｙ＝ｒ・ｓｉｎθ・ｓｉｎφ
　　ｚ＝ｒ・ｃｏｓθ
の関係があり、相互に変換可能である。
　具体的には、先に求められた近似楕円体の重心を球座標系の原点Ｏとする。この原点Ｏは、原画像データにおけるｘｙｚ直交座標系の原点と一致する必要はない。そして、適宜の変換処理により直交座標から球座標への座標変換を行う（ステップＳ１０９）。このように座標変換を行うことにより、ｘｙｚ直交座標系で特定される三次元空間内の各位置が、ｒθφ球座標系で表されることになる。楕円体の重心は、例えば楕円体の長軸または短軸を表す線分の中点とすることができる。
　このような座標変換により、細胞Ｃ１（Ｃ２）の露出表面および近似楕円体Ｅａ１（Ｅａ２）の表面の各点の座標は極座標により表される。図９Ｃおよび図９Ｄに示すように、１つの動径ｒ方向における細胞Ｃ１（Ｃ２）の露出表面と近似楕円体Ｅａ１（Ｅａ２）の表面との距離Δｒを求めると、この距離Δｒは当該方向における細胞と楕円体との乖離度を表す値となる。つまり、距離Δｒが小さければ細胞の表面と楕円体の表面との乖離は小さく、大きいほど乖離が顕著であると言える。
　細胞Ｃ１全体としての近似楕円体Ｅ１に対する乖離度を表すには、例えば露出表面の全体について動径方向ごとに上記と同様に距離Δｒを求め、その最大値、最大値と最小値との差、絶対値の積分値等の少なくとも１つを指標値とすることができる（ステップＳ１１０）。位置を球座標で表したことにより、このような細胞とその近似楕円体との間の定量的な比較演算を簡単に行うことが可能になっている。これらの指標値は、細胞の歪みの程度を表したものといえる。このような指標値が小さいほど、細胞はその表面が滑らかで歪みが小さく、良好な細胞であると推定することができる。
　上記のように、個々の構造体の形状面での評価は、ＯＣＴ撮像により得られる三次元像を二値化した二値化画像に基づいて行うことが可能である。一方、個々の構造体の内部における充実の程度やその均一性を評価するためには、構造体の密度の違いに起因する三次元オブジェクトの濃淡に関わる情報が必要である。
　ここで、上記のように個々の細胞およびその他の構造体が個別の領域に分割された二値化画像をマスクとして用い、原画像データで表される濃淡を含む三次元像に対して作用させる。そうすると、原画像に含まれる三次元オブジェクトを、領域ごとに個別に切り出すことが可能である。すなわち、原画像と、領域分割された二値化画像との画素ごとの論理積に相当する画像を作成することで（ステップＳ１１１）、原画像を領域分割した画像が得られる。また、原画像中の三次元オブジェクトを、上記した領域分割により画定された境界で分割することによっても同様の結果が得られる。
　これにより、分割された個々の領域においては、その濃淡は当該領域内の各画素の画素値（輝度値）によって表される。この情報を用いて、その内部構造を個別に評価することが可能となる。例えば、分割された領域内の平均濃度は当該領域内の物質の平均密度を表す指標値であり、例えば核を含む細胞とこれを含まないフラグメント等との識別に役立つ情報となる。また例えば分割された各領域につき、その内部の輝度値に基づき濃度ヒストグラムを作成する（ステップＳ１１２）。こうすれば、内部構造を推定するのに有用な情報を得ることが可能になる。１つの領域全体について濃度ヒストグラムを作成すれば、当該領域全体の性質を表すことができる。また、１つの領域のうち一部について濃度ヒストグラムを作成することで、領域内部の構造をより詳細に示すことが可能となる。
　図１０Ａおよび図１０Ｂは濃度ヒストグラムの例を示す図である。図１０Ａは、領域分割により画定された細胞Ｃ１の内部の２つの領域Ｒ１，Ｒ２についての濃度ヒストグラムの例を示す。領域Ｒ１についての濃度ヒストグラムでは、比較的低い濃度域に幅の狭いピークが現れている。この領域が、物質密度の低い、例えば空洞であると推定することができる。一方、領域Ｒ２からは比較的高い濃度域に幅の広いピークが現れており、物質が高濃度かつ局在して含まれていることが示される。例えば細胞内の小器官が集中している領域や細胞分裂のための動きが活発な領域において、このような傾向が見られる。
　また、細胞内の核の領域とその他の領域との間でも、濃度ヒストグラムに大きな違いが生じると考えられる。このように、濃度ヒストグラムを求めることにより、例えばピークの位置、高さ、幅など細胞の内部構造に関する手がかりとなる情報をユーザに提供することが可能になる。
　このような濃度ヒストグラムを求める手法は、胚盤胞の評価にも有用なものである。すなわち、図１０Ｂに示すように、胚盤胞では薄い栄養外胚葉Ｔの内部に内細胞塊Ｉが形成されている。状態の良い胚では内細胞塊Ｉが多数の細胞で密に構成されているのに対し、そうでない胚の内細胞塊では、細胞数が少ないあるいは密度が低い疎な状態である場合がある。したがって、領域分割により特定された内細胞塊Ｉについて濃度ヒストグラムを求めたとき、状態の良い胚では、図１０Ｂ下部のグラフにおいて実線で示すようにピークの広がりが小さい。これに対し、状態の良くない胚では、破線で示すようにピークの広がりが大きくなる。このように、内細胞塊Ｉについても、ピークの位置、高さ、幅などを指標する情報を、胚の評価に用いることができる。例えばこれらの指標値をGardner分類に適用することで、従来は主観的判断に頼っていた分類を定量的に行うことが可能となる。
　以上のように、この実施形態では、胚をＯＣＴ撮像することにより得られる三次元像から、細胞、フラグメント等の各種の構造体に対応する領域を個別に特定することが可能である。こうして特定された領域ごとにその特徴を定量的に示す指標値を算出することで、各領域の評価に有用な情報が客観的な数値として得られる。このような指標値を例えば表示部３５２に表示する等の態様でユーザに提供することで、この実施形態では、ユーザによる胚の評価作業を効果的に支援することが可能である。
　また、三次元像に含まれる構造体としては、割球、細胞、フラグメント等種々の種類のものが含まれ得る。しかしながら、こうして求められる指標値に対して予め適宜の閾値を設定しておけば、構造体がどの種類のものかを自動的に判定する画像判定方法として実施することが可能となる。
　なお、本発明は上記した実施形態に限定されるものではなく、その趣旨を逸脱しない限りにおいて上述したもの以外に種々の変更を行うことが可能である。例えば、上記実施形態の画像処理装置１は、試料ＳｐをＯＣＴ撮像する機能、および撮像データから出力画像を作成し出力する機能を有するものである。しかしながら、本発明の画像処理方法は、自身は撮像機能を持たず、撮像機能を有する他の装置での撮像により得られた撮像データを取得したコンピュータ装置によって実行することも可能である。これを可能とするために、図５の各処理ステップのうちステップＳ１０２ないしＳ１１２をコンピュータ装置に実行させるためのソフトウェアプログラムとして、本発明が実施されてもよい。
　このようなプログラムの配布は、例えばインターネット等の電気通信回線を介してダウンロードする形式によって行うことが可能である。また、当該プログラムを記録したコンピュータ読み取り可能な記録媒体を配布することによっても可能である。また、既存のＯＣＴ撮像装置にインターフェースを介してこのプログラムを読み込ませることで、当該撮像装置により本発明を実施することも可能となる。
　また例えば、上記実施形態では、割球の集合体を構成する細胞を個別に評価する目的に供するために、Local Thickness演算で区分された領域のうちフラグメント等に該当する第２領域を消去する処理が実行されている。しかしながら、例えばフラグメントを詳しく観察する目的のために、第２領域を残して第１領域を消去することも可能である。また、第１領域と第２領域とを個別に表す画像が作成されてもよい。
　また例えば、上記した実施形態における画像処理（図５）では、胚の三次元像を領域分割するための処理（ステップＳ１０４～Ｓ１０６）、分割された領域のサイズ算出（ステップＳ１０７）、個々の細胞を楕円体近似する処理（ステップＳ１０８）、細胞の表面を極座標により表現する処理（ステップＳ１０９）、細胞と近似楕円体との乖離度を求める処理（ステップＳ１１０）、領域内の濃度分布を求める処理（ステップＳ１１１～Ｓ１１２）がこの順番で実行される。しかしながら、これらの全てが実行されることは必須でなく、目的に応じ必要な処理が適宜選択されて実行されてもよい。また、実行順序も上記に限定されるものではない。
　また、上記実施形態は、ＯＣＴ撮像により得られる三次元像から、ユーザによる評価作業に有用な各種の定量的情報を提示するものである。ここで、ユーザによるそれらの情報の利用形態については、上記に例示したものに限定されず任意である。
　以上、具体的な実施形態を例示して説明してきたように、本発明に係る画像処理方法においては、Local Thickness演算により三次元像内の各点に与えられる階調値が指標値とされてもよい。一般ユーザが利用可能な態様で公開されているLocal Thickness演算アルゴリズムでは、オブジェクトの大きさが階調値によって表現される。このような階調値を、上記した処理における指標値として好適に適用することが可能である。
　また例えば、領域分割する工程は、複数の領域の間における境界を画定し、原画像データが表す三次元像を画定された境界で複数に分割するように構成されてもよい。これらの境界で三次元像中のオブジェクトを分割することで、割球の個々の細胞およびその他の構造物を互いに分離し個別に評価することが可能となる。
　また例えば、分割された領域について、当該領域の表面のうち他の領域と接していない露出表面が楕円体近似されるように構成されてもよい。細胞は不定形であるが概ね楕円体形状を有している、そのため、細胞を楕円体近似することで、その形状を具体的な数値として表すことが可能になる。この場合、他の細胞と接している部分では細胞の表面が変形している。そのため、そのような変形の少ない露出表面を楕円体近似することで、表面の変形に起因する誤差を低減することができる。
　この場合、例えば、楕円体近似により求められた近似楕円体の表面と露出表面との乖離度を指標する指標値がさらに算出されてもよい。近似楕円体は細胞の理想的な形状を示すものと言えるので、これからの乖離の大きさは、細胞の歪みの大きさや表面の凹凸の大きさを示すものである。このような形状に関する特徴を定量的に示す指標値があれば、ユーザによる評価作業に有用なものとなる。
　さらに、例えば、露出表面および近似楕円体の表面の各位置の座標を近似楕円体の重心を原点とする極座標によって表し、同一動径上にある露出表面上の点と近似楕円体の表面上の点との距離に基づき指標値が求められてもよい。楕円体の重心から見た細胞と近似楕円体との距離は、それらの間の乖離度を端的に表す情報として利用可能である。
　また例えば、分割された領域内の画素各々の画素値のヒストグラムを算出するようにしてもよい。このような構成によれば、細胞の内部構造を定量化した情報を作成することが可能である。
　また、この発明に係る画像処理方法は、上記方法により胚の三次元像を複数の領域に分割し、そのうち指標値が第１閾値より大きい第１領域を割球の集合体と判定するように構成することができる。胚において正常な割球は他の構造体に比べて大きな構造を有していると考えられる。したがって、比較的大きな構造と推定される第１領域を割球の集合体と判定することは極めて合理的である。
　この場合、第１領域からさらに分割された１つの領域を１つの細胞と判定してもよい。第１領域を複数の細胞からなる割球と見なしたとき、それをさらに分割した個々の領域は、それぞれ１つの細胞に対応すると考えるのが妥当である。
　また例えば、指標値が第１閾値より小さい第２領域に区分された領域のうち球形度が第２閾値より大きい領域をフラグメントと判定してもよい。フラグメントは細胞に比べて構造が小さいが、一定の球形度を有している。したがって、第２領域に区分される領域のうち比較的高い球形度の領域をフラグメントと見なすことができる。こうすることで、フラグメントとその他の構造体とを区別することが可能になる。
　以上、特定の実施例に沿って発明を説明したが、この説明は限定的な意味で解釈されることを意図したものではない。発明の説明を参照すれば、本発明のその他の実施形態と同様に、開示された実施形態の様々な変形例が、この技術に精通した者に明らかとなるであろう。故に、添付の特許請求の範囲は、発明の真の範囲を逸脱しない範囲内で、当該変形例または実施形態を含むものと考えられる。
　この発明は、培養された胚の状態を評価する作業を支援する目的に適用することが可能であり、例えば不妊治療における体外受精、人工授精の成功確率を高めるために利用することができる。
　１　画像処理装置
　１０　保持部
　２０　撮像ユニット
　２１　光源
　２２　ビームスプリッタ
　２４　基準ミラー
　２６　光検出器
　３０　制御ユニット
　３３　信号処理部
　３４　３Ｄ復元部
　３５２　表示部
　Ｓｐ　試料

Claims (12)

  1. 　培養された胚を光干渉断層撮像して得られる、前記胚の三次元像に対応する原画像データを取得する工程と、
    　前記原画像データに基づき、前記三次元像を複数の領域に領域分割する工程と
    を備え、
    　前記領域分割する工程では、
    　前記三次元像に対しLocal Thickness演算を実行して前記三次元像に含まれるオブジェクトの大きさを指標する指標値を求め、該指標値が所定の第１閾値より大きい領域と小さい領域とに前記三次元像を分割し、
    　分割された領域の各々を、Watershedアルゴリズムにより分割する、
    画像処理方法。
  2. 　前記Local Thickness演算により前記三次元像内の各点に与えられる階調値が前記指標値とされる請求項１に記載の画像処理方法。
  3. 　前記領域分割する工程では、前記複数の領域の間における境界を画定し、前記原画像データが表す前記三次元像を画定された境界で複数に分割する請求項１または２に記載の画像処理方法。
  4. 　分割された前記領域について、当該領域の表面のうち他の領域と接していない露出表面が楕円体近似される請求項１ないし３のいずれかに記載の画像処理方法。
  5. 　前記楕円体近似により求められた近似楕円体の表面と、前記露出表面との乖離度を指標する指標値を算出する請求項４に記載の画像処理方法。
  6. 　前記露出表面および前記近似楕円体の表面の各位置の座標を、前記近似楕円体の重心を原点とする極座標によって表し、同一動径上にある前記露出表面上の点と前記近似楕円体の表面上の点との距離に基づき、前記指標値を求める請求項５に記載の画像処理方法。
  7. 　分割された前記領域内の画素各々の画素値のヒストグラムを算出する請求項１ないし６のいずれかに記載の画像処理方法。
  8. 　請求項１ないし７のいずれかに記載の画像処理方法により、前記胚の三次元像を複数の領域に分割し、前記複数の領域のうち前記指標値が前記第１閾値より大きい第１領域を割球の集合体と判定する画像処理方法。
  9. 　前記第１領域から分割された１つの領域を１つの細胞と判定する請求項８に記載の画像処理方法。
  10. 　前記指標値が前記第１閾値より小さい第２領域から分割された領域のうち、球形度が第２閾値より大きい領域をフラグメントと判定する請求項８または９に記載の画像処理方法。
  11. 　培養された胚を光干渉断層撮像して得られる、前記胚の三次元像に対応する原画像データを取得する工程と、
    　前記原画像データに基づき、前記三次元像を複数の領域に領域分割する工程と
    を備え、
    　前記領域分割する工程では、
    　前記三次元像に対しLocal Thickness演算を実行して前記三次元像に含まれるオブジェクトの大きさを指標する指標値を求め、該指標値が所定の第１閾値より大きい領域と小さい領域とに前記三次元像を分割し、
    　分割された領域の各々を、Watershedアルゴリズムにより分割する、
    画像処理を、コンピュータに実行させるためのプログラム。
  12. 　請求項１１に記載のプログラムを記録した、コンピュータ読み取り可能な記録媒体。
PCT/JP2019/003106 2018-01-31 2019-01-30 画像処理方法、プログラムおよび記録媒体 WO2019151302A1 (ja)

Priority Applications (3)

Application Number Priority Date Filing Date Title
US16/966,488 US11710237B2 (en) 2018-01-31 2019-01-30 Image processing method and recording medium
EP19747001.6A EP3748580B1 (en) 2018-01-31 2019-01-30 Image processing method, program and recording medium
CN201980009540.7A CN111630564B (zh) 2018-01-31 2019-01-30 图像处理方法、程序及记录介质

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2018-015378 2018-01-31
JP2018015378A JP6996709B2 (ja) 2018-01-31 2018-01-31 画像処理方法、画像判定方法、プログラムおよび記録媒体

Publications (1)

Publication Number Publication Date
WO2019151302A1 true WO2019151302A1 (ja) 2019-08-08

Family

ID=67478725

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2019/003106 WO2019151302A1 (ja) 2018-01-31 2019-01-30 画像処理方法、プログラムおよび記録媒体

Country Status (4)

Country Link
US (1) US11710237B2 (ja)
EP (1) EP3748580B1 (ja)
JP (1) JP6996709B2 (ja)
WO (1) WO2019151302A1 (ja)

Families Citing this family (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP7382289B2 (ja) * 2020-06-19 2023-11-16 株式会社Ｓｃｒｅｅｎホールディングス 画像処理方法、プログラムおよび記録媒体
KR102530717B1 (ko) * 2020-06-19 2023-05-09 가부시키가이샤 스크린 홀딩스 화상 처리 방법, 프로그램 및 기록 매체
EP3949830A1 (en) * 2020-08-05 2022-02-09 Leica Microsystems Inc. Control system for an oct imaging system, oct imaging system and method for imaging
JP2022143662A (ja) 2021-03-18 2022-10-03 株式会社Ｓｃｒｅｅｎホールディングス 受精卵の発生ステージ判定方法、プログラム、記録媒体、撮像方法および撮像装置
WO2022226647A1 (en) * 2021-04-27 2022-11-03 The Governing Council Of The University Of Toronto Systems and methods for detection and 3d imaging of cell junctions of blastocysts and other biological structures
JP2023018827A (ja) * 2021-07-28 2023-02-09 株式会社Ｓｃｒｅｅｎホールディングス 画像処理方法、プログラムおよび記録媒体
EP4145389A1 (en) 2021-09-01 2023-03-08 SCREEN Holdings Co., Ltd. Image processing method, fertilized egg evaluation method, computer program and recording medium
JP2023035518A (ja) 2021-09-01 2023-03-13 株式会社Ｓｃｒｅｅｎホールディングス 画像処理方法、コンピュータープログラムおよび記録媒体
JP2023125282A (ja) 2022-02-28 2023-09-07 株式会社Ｓｃｒｅｅｎホールディングス 解析方法および解析装置
WO2024062953A1 (ja) * 2022-09-22 2024-03-28 株式会社レゾナック 画像処理方法及び画像処理プログラム
JP2024046836A (ja) * 2022-09-26 2024-04-05 株式会社Ｓｃｒｅｅｎホールディングス 分類方法およびコンピュータプログラム
CN115511831B (zh) * 2022-09-27 2023-04-04 佳木斯大学 一种组织胚胎病理切片的数据分析处理系统及方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2006338191A (ja) * 2005-05-31 2006-12-14 Olympus Corp 画像処理装置及び領域分割プログラム
JP2014016666A (ja) * 2012-07-05 2014-01-30 Olympus Corp 細胞分裂過程追跡装置、及び細胞分裂過程追跡プログラム
JP2015130806A (ja) * 2014-01-09 2015-07-23 大日本印刷株式会社 成育情報管理システム及び成育情報管理プログラム
JP2016146132A (ja) * 2015-02-09 2016-08-12 学校法人日本大学 形状特徴抽出方法、形状特徴抽出処理装置、形状記述方法及び形状分類方法
JP2017521067A (ja) 2014-07-01 2017-08-03 アンスティチュ ナショナル ドゥ ラ サンテ エ ドゥ ラ ルシェルシュ メディカル 胚の品質の三次元再構成及び決定のための方法

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8515143B2 (en) * 2009-01-09 2013-08-20 Dai Nippon Printing Co., Ltd. Embryo quality evaluation assistance system, embryo quality evaluation assistance apparatus and embryo quality evaluation assistance method
JP2012531584A (ja) * 2009-06-25 2012-12-10 フェイズ ホログラフィック イメージング ペーホーイー アーベー デジタルホログラフィック撮像による卵又は胚の分析
GB201810634D0 (en) * 2018-06-28 2018-08-15 Vitrolife As Methods and apparatus for assessing embryo development

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2006338191A (ja) * 2005-05-31 2006-12-14 Olympus Corp 画像処理装置及び領域分割プログラム
JP2014016666A (ja) * 2012-07-05 2014-01-30 Olympus Corp 細胞分裂過程追跡装置、及び細胞分裂過程追跡プログラム
JP2015130806A (ja) * 2014-01-09 2015-07-23 大日本印刷株式会社 成育情報管理システム及び成育情報管理プログラム
JP2017521067A (ja) 2014-07-01 2017-08-03 アンスティチュ ナショナル ドゥ ラ サンテ エ ドゥ ラ ルシェルシュ メディカル 胚の品質の三次元再構成及び決定のための方法
JP2016146132A (ja) * 2015-02-09 2016-08-12 学校法人日本大学 形状特徴抽出方法、形状特徴抽出処理装置、形状記述方法及び形状分類方法

Non-Patent Citations (1)

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

Also Published As

Publication number Publication date
EP3748580A4 (en) 2021-10-27
EP3748580A1 (en) 2020-12-09
US11710237B2 (en) 2023-07-25
US20210056704A1 (en) 2021-02-25
JP2019133429A (ja) 2019-08-08
JP6996709B2 (ja) 2022-01-17
CN111630564A (zh) 2020-09-04
EP3748580B1 (en) 2024-04-17

Similar Documents

Publication Publication Date Title
WO2019151302A1 (ja) 画像処理方法、プログラムおよび記録媒体
JP6867629B2 (ja) 画像処理方法、プログラムおよび記録媒体
JP4847530B2 (ja) 病変組織のより有意義な解剖学的特徴を明らかにする生体画像データの処理方法
KR102580984B1 (ko) 화상 처리 방법, 프로그램 및 기록 매체
US20230063443A1 (en) Image processing method and recording medium
US20220301155A1 (en) Fertilized egg development stage determination method, recording medium, imaging method and imaging apparatus
JP7382289B2 (ja) 画像処理方法、プログラムおよび記録媒体
Taieb et al. Classification of tissue biopsies by raman spectroscopy guided by quantitative phase imaging and its application to bladder cancer
JP6895297B2 (ja) 細胞塊の評価方法及び細胞塊の状態解析装置
JP2023066991A (ja) 画像処理方法、受精卵の評価方法、コンピュータープログラムおよび記録媒体
JP2023035517A (ja) 画像処理方法、受精卵の評価方法、コンピュータープログラムおよび記録媒体
JP2022143660A (ja) 画像処理方法、プログラムおよび記録媒体、ならびに画像処理装置
CN111630564B (zh) 图像处理方法、程序及记录介质
JP7382290B2 (ja) 画像処理方法、プログラムおよび記録媒体
US20230061181A1 (en) Image processing method, fertilized egg evaluation method and recording medium
Sawyer Multimodal Optical Imaging for Tissue Characterization and Disease Diagnosis
Weng Integrating Coherent Anti-Stokes Raman Scattering Imaging and Deep Learning Analytics for High Precision, Real Time, Label Free Cancer Diagnosis

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

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

ENP Entry into the national phase

Ref document number: 2019747001

Country of ref document: EP

Effective date: 20200831