JP2022001984A - Image processing method, program, and storage medium - Google Patents

Abstract

Kind Code: A1 A quantitative information useful for analyzing the surface structure of a cell mass is obtained from OCT image data of a cell mass having cavities inside. Kind Code: A1 An image processing method according to the present invention comprises a step of acquiring three-dimensional image data representing a three-dimensional image of a cell mass by optical coherence tomography (step S111); obtaining a thickness distribution representing the thickness of the cell mass between the outer surface facing the external space and the inner surface facing the cavity of the cell mass for each position (step S116); and a step of detecting a protruding portion having a thickness protruding from the surrounding area based on the acquired result (step S118). [Selection drawing] Fig. 3

Description

The present invention relates to image processing for imaging a cell mass, particularly a cell mass having a cavity inside, and analyzing its structure.

For example, in assisted reproductive technology for the purpose of infertility treatment, embryos (fertilized eggs) fertilized in vitro and cultured for a certain period of time are returned to the body. However, its pregnancy success rate (in assisted reproductive technology) is not always high, and the mental and financial burden on the patient is also great. In order to solve this problem, a method for accurately determining the state of the embryo to be cultured is being sought.

For example, in an embryo that has progressed to the blastocyst stage, a cavity structure called a blastocoel is formed inside, and a layer of cells called a vegetative ectoderm surrounds the cavity structure. The number of cells constituting the vegetative ectoderm is an index for evaluating the state of the embryo, and therefore a technique for non-invasively counting the number is required.

In addition, for example, organoids in the intestinal epithelium are known to constitute a cell population having an internal cavity. When this is used as an evaluation model for the cell membrane transport system, it is ideal that a cavity covered with one layer of cell layer is formed, but in reality, a multi-layered cell layer is formed around the cavity. There are many. Therefore, there is a need for a technique for quantitatively evaluating the structure of the cell layer on the surface of the intestinal epithelial organoid, specifically, how the cells form the cell layer covering the cavity.

As a technique expected to meet such a demand, the applicant of the present application has previously disclosed Patent Documents 1 and 2. Patent Document 1 describes trophectoderm and inner cells from a three-dimensional image of an embryo (fertilized egg) imaged by a non-invasive tomographic imaging technique such as optical coherence tomography (OCT). A method for identifying and dividing a mass is described. Further, in Patent Document 2, the embryo imaged by OCT is represented by polar coordinates with the position of the center of gravity as the origin, and the reflected light intensity in each radial direction is visualized by a two-dimensional map replaced with a luminance value for observation. Techniques for assisting the embryo evaluation work by the person are described.

Japanese Unexamined Patent Publication No. 2019-133249 Japanese Unexamined Patent Publication No. 2019-132710

The above-mentioned conventional technique presents the three-dimensional structure of the observation object to the observer in an easy-to-understand manner, but it has not been possible to automatically extract quantitative information effective for supporting observation and evaluation work. Not in. In this respect, there is room for improvement in the above-mentioned prior art. More specifically, there is a need for a technique capable of obtaining quantitative information on the surface structure of a cell mass having an inner cavity such as a blastocyst stage embryo or an intestinal epithelial organoid.

The present invention has been made in view of the above problems, and is useful for analyzing the surface structure of a cell mass from the image data obtained by OCT imaging in an image process in which an image of a cell mass having a cavity inside is a processing target. The purpose is to provide a technology that can request quantitative quantitative information.

One aspect of the present invention is an image processing method for processing an image of a cell mass having a cavity inside, and in order to achieve the above object, the cell mass is imaged by optical interference tomography and the cell mass is imaged. The cell between the outer surface of the cell mass facing the outer space and the inner surface facing the cavity based on the step of acquiring the three-dimensional image data representing the three-dimensional image of the cell mass. It is provided with a step of acquiring a thickness distribution representing the thickness of the mass for each position and a step of detecting a protrusion portion having a thickness protruding compared to the surrounding region based on the acquisition result of the thickness distribution. There is.

Further, another aspect of the present invention is a program for causing a computer to execute each of the above steps. Further, another aspect of the present invention is a computer-readable recording medium that stores the above program.

Due to the internal structure of the cell such as the cell nucleus, each cell has a structure in which the central portion is swollen as compared with the peripheral portion. Therefore, in the cell layer formed by connecting cells in the lateral direction to each other, unevenness due to the shape of each cell occurs. In particular, in a cell layer composed of one layer of cells so as to cover the internal cavity, such as an embryo in the blastocyst stage (fertilized egg), a protrusion corresponding to each cell appears in the layer. become. Therefore, the number of protrusions can be an index showing the position and number of cells constituting the layer. In addition, the spacing between protrusions can be an indicator of the size of individual cells.

In the above invention, the distribution of the thickness of the cell mass between the outer surface of the cell mass and the inner surface facing the inner cavity can be obtained from the three-dimensional image data of the cell mass imaged by OCT. This thickness distribution can represent the thickness of the cell layer covering the cavity by position. The thickness of the cell layer is not always uniform, for example, in a cell layer composed of one layer of cells, it is thicker in the central part of each cell and thinner in the peripheral part. Therefore, if a protrusion site having a thickness larger than that of the surroundings is detected from the obtained thickness distribution, it can be said that the site is highly probable to indicate the position of one cell.

In this way, by obtaining the thickness distribution of the cell layer covering the internal cavity and detecting the protrusions contained therein, information that quantitatively represents the position, number, size, etc. of the cells constituting the cell layer. Can be obtained.

As described above, according to the present invention, from the image data obtained by OCT imaging a cell mass having a cavity inside, the position and number of individual cells useful for analysis of the surface structure of the cell mass are related. It is possible to obtain quantitative information.

It is a figure which shows the structural example of the image processing apparatus which executes the image processing method of this invention. It is a figure which shows typically the structure of the embryo as a sample in this embodiment. It is a flowchart which shows the image processing in this embodiment. It is a figure which shows the relationship between the focal depth and the area of the zona pellucida extracted. It is a figure which shows the correspondence relation between the XYZ Cartesian coordinate system and the polar coordinate system. It is a figure which illustrates the method of determining the thickness of the vegetative ectoderm. It is a flowchart which shows the process of determining the thickness of the vegetative ectoderm. It is a figure which shows an example of a two-dimensional mapping method. It is a figure which shows another example of a 2D mapping method. It is a figure which shows the example of the 3D mapping method. It is a figure which shows the principle of a peak detection processing. It is a flowchart which shows the peak detection processing. It is a figure which shows the example of the peak position detected on a two-dimensional map. It is a figure which shows the example of the peak position detected on a three-dimensional map.

FIG. 1 is a diagram showing a configuration example of an image processing apparatus suitable as an execution subject of the image processing method according to the present invention. This image processing apparatus 1 performs a tomographic image of a sample carried in a liquid, for example, an embryo (fertilized egg) cultured in a culture solution, and image-processes the obtained tomographic image to obtain a cross section of a sample. Create a cross-sectional image showing the structure. In addition, a stereoscopic image of the sample is created from a plurality of cross-sectional images. In order to show the directions in each of the following figures in a unified manner, the XYZ orthogonal coordinate axes are set as shown in FIG. Here, the XY plane represents the horizontal plane. Further, the Z axis represents the vertical axis, and more specifically, the (−Z) direction represents the vertical downward direction.

The image processing device 1 includes a holding unit 10. The holding unit 10 holds the sample container 11 containing the sample S to be an image pickup object in a horizontal posture. The sample container 11 is, for example, a flat-bottomed, shallow dish-shaped container called a dish. Alternatively, for example, the sample container 11 is a well plate in which a large number of recesses (wells) W capable of supporting a liquid are formed on the upper surface of the plate-shaped member. A medium M such as a culture solution is injected into the sample container 11, and the sample S is supported inside the medium M.

The imaging unit 20 is arranged below the sample container 11 held by the holding unit 10. An optical coherence tomography (OCT) device capable of non-contact and non-destructive (non-invasive) imaging of a tomographic image of an image to be imaged is used for the image pickup unit 20. As will be described in detail later, the image pickup unit 20 which is an OCT device includes a light source 21 that generates illumination light for an object to be imaged, an optical fiber coupler 22, an object optical system 23, a reference optical system 24, and a spectroscope 25. And an optical detector 26.

The imaging unit 20 further includes a microscope imaging unit 28 for performing optical microscope imaging. More specifically, the microscope image pickup unit 28 includes an image pickup optical system 281 and an image pickup element 282. The imaging optical system 281 includes an objective lens whose focus is on the sample S in the sample container 11. Further, as the image pickup device 282, for example, a CCD image pickup device, a CMOS sensor, or the like can be used. As the microscope imaging unit 28, one capable of bright-field imaging or phase-difference imaging is preferable. The object optical system 23 and the microscope imaging unit 28 are supported by a support member (not shown) that can move in the horizontal direction, and the position in the horizontal direction can be changed.

The image processing device 1 further includes a control unit 30 that controls the operation of the device, and a drive unit 40 that drives a movable portion of the image pickup unit 20. The control unit 30 includes a CPU (Central Processing Unit) 31, an A / D converter 32, a signal processing unit 33, an image pickup control unit 34, an interface (IF) unit 35, an image memory 36, and a memory 37.

The CPU 31 controls the operation of the entire device by executing a predetermined control program, and the control program executed by the CPU 31 and the data generated during the processing are stored in the memory 37. The A / D converter 32 converts the signal output from the photodetector 26 of the image pickup unit 20 and the image pickup element 282 according to the amount of received light into digital data. The signal processing unit 33 performs signal processing described later based on the digital data output from the A / D converter 32 to create an image of the object to be imaged. The image data thus created is appropriately stored and stored in the image memory 36.

The image pickup control unit 34 controls each part of the image pickup unit 20 to execute the image pickup process. Specifically, the imaging control unit 34 selectively positions the object optical system 23 for OCT imaging and the microscope imaging unit 28 for optical microscope imaging at a predetermined imaging position within the imaging field of view of the imaging sample S. .. When the object optical system 23 is positioned at the imaging position, the imaging control unit 34 causes the imaging unit 20 to execute an OCT imaging process described later to acquire three-dimensional image data representing the three-dimensional structure of the sample S. On the other hand, when the microscope image pickup unit 28 is positioned at the image pickup position, the image pickup control unit 34 optically obtains two-dimensional image data representing a planar image of the sample S imaged on the light receiving surface of the image pickup element 282 by the image pickup optical system 281. Have the image pickup unit 28 acquire the image.

The interface unit 35 is responsible for communication between the image processing device 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 operation input from the user and notifying the user of various information. For this purpose, the interface unit 35 has an input device 351 such as a keyboard, a mouse, and a touch panel capable of receiving operation inputs related to device function selection and operating condition setting, and a fault created by the signal processing unit 33. It is connected to a display unit 352 made of, for example, a liquid crystal display, which displays various processing results such as an image and a stereoscopic image created by the 3D restoration unit 34.

In the image pickup unit 20, a low coherence light beam containing a wide-band wavelength component is emitted from a light source 21 having a light emitting element such as a light emitting diode or a superluminescent diode (SLD). For the purpose of imaging a sample such as cells, it is preferable to use, for example, near infrared rays in order to allow the incident light to reach the inside of the sample.

The light source 21 is connected to an optical fiber 221 which is one of the optical fibers constituting the optical fiber coupler 22, and the low coherence light emitted from the light source 21 is transferred to the two optical fibers 222 and 224 by the optical fiber coupler 22. It is branched into the light of. The optical fiber 222 constitutes an object-based optical path. More specifically, the light emitted from the end of the optical fiber 222 is incident on the object optical system 23.

The object optical system 23 includes a collimator lens 231 and an objective lens 232. The light emitted from the end of the optical fiber 222 enters the objective lens 232 via the collimator lens 231. The objective lens 232 has a function of converging the light (observation light) from the light source 21 on the sample and a function of condensing the reflected light emitted from the sample and directing it toward the optical fiber coupler 22. Although a single objective lens 232 is shown in the figure, a plurality of optical elements may be combined. The reflected light from the object to be imaged enters the optical fiber 222 as signal light via the objective lens 232 and the collimator lens 231. The optical axis of the objective lens 232 is orthogonal to the bottom surface of the sample container 11, and in this example, the optical axis direction coincides with the vertical axis direction.

The CPU 31 gives a control command to the image pickup control unit 34, and the image pickup control unit 34 causes the image pickup unit 20 to move in a predetermined direction in response to the control command. More specifically, the image pickup control unit 34 moves the image pickup unit 20 in the horizontal direction (XY direction) and the vertical direction (Z direction). The horizontal movement of the image pickup unit 20 changes the image pickup range in the horizontal direction. Further, due to the vertical movement of the image pickup unit 20, the focal position of the objective lens 232 in the optical axis direction changes with respect to the sample S which is the object to be imaged.

A part of the light incident on the optical fiber coupler 22 from the light source 21 is incident on the reference optical system 24 via the optical fiber 224. The reference optical system 24 includes a collimator lens 241 and a reference mirror 243, which together with the optical fiber 224 form a reference system optical path. Specifically, the light emitted from the end of the optical fiber 224 is incident on the reference mirror 243 via the collimator lens 241. The light reflected by the reference mirror 243 is incident on the optical fiber 224 as reference light.

The reference mirror 243 is supported by an advance / retreat member (not shown) operated by a control command from the image pickup control unit 34, and can move forward / backward in the Y direction. The optical path length of the reference light reflected by the reference mirror 243 is adjusted by moving the reference mirror 243 in the Y direction, that is, in the approaching and separating directions with respect to the collimator lens 241.

The reflected light (signal light) reflected by the surface or internal reflection surface of the sample and the reference light reflected by the reference mirror 243 are mixed by the optical fiber coupler 22 and incident on the optical detector 26 via the optical fiber 226. do. At this time, interference due to the phase difference occurs between the signal light and the reference light, but the spectral spectrum of the interference light differs 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 detecting the amount of light by splitting the interference light for each wavelength and performing a Fourier transform on the detected interference signal. The OCT imaging technique based on such a principle is called a Fourier domain OCT (FD-OCT).

In the image pickup unit 20 of this embodiment, the spectroscope 25 is provided on the optical path of the interference light from the optical fiber 226 to the photodetector 26. As the spectroscope 25, for example, one using a prism, one using a diffraction grating, or the like can be used. The interference light is separated by the spectroscope 25 for each wavelength component and received by the photodetector 26.

By Fourier-converting the interference signal output from the photodetector 26 according to the interference light detected by the photodetector 26, the reflected light intensity in the depth direction at the incident position of the illumination light in the sample, that is, in the Z direction. The distribution is required. By scanning the light beam incident on the sample container 11 in the X direction, the reflected light intensity distribution in a plane parallel to the XZ plane can be obtained, and from the result, a tomographic image of the sample having the plane as a cross section can be created. can. Since the principle is well known, detailed description thereof will be omitted.

Further, by taking a tomographic image each time while changing the beam incident position in the Y direction in multiple stages, it is possible to obtain a large number of tomographic images obtained by taking a tomographic image of the sample in a cross section parallel to the XZ plane. If the scanning pitch in the Y direction is reduced, image data having a resolution sufficient to grasp the three-dimensional structure of the sample can be obtained. From these tomographic image data, it is possible to create three-dimensional image data (so-called voxel data) corresponding to the three-dimensional image of the sample.

As described above, the image processing apparatus 1 has a function of acquiring an image of the sample S supported on the sample container 11 together with the medium M. As the image, it is possible to acquire two-dimensional image data obtained by optical microscope imaging, tomographic image data obtained by OCT imaging, and three-dimensional image data based on the tomographic image data.

Hereinafter, an embodiment of the image processing method according to the present invention, which can be executed by using the image processing apparatus 1 configured as described above, will be described. In the image processing in this embodiment, a fertilized egg in the blastocyst stage (hereinafter, simply referred to as “embryo”) is imaged as sample S, and from the image, the zona pellucida, which is the main structure constituting the embryo, and ectoderm It has the content of individually extracting the regions corresponding to the germ layer and the inner cell mass.

Based on the data obtained by this, it becomes possible to effectively support the evaluation work of the embryo by the user (specifically, a doctor or an embryo cultivator). For example, in the culture of embryos for the purpose of infertility treatment, the image processing method of the present embodiment can be applied for the purpose of obtaining findings for determining whether or not the culture is proceeding well.

FIG. 2 is a diagram schematically showing the structure of an embryo as a sample in this embodiment. As is already known, cleavage is started when an egg is fertilized, and a blastocyst is formed through a state called a morula. FIG. 2A schematically shows the internal structure of the embryo in the blastocyst stage. In the blastocyst stage, embryo E has a cavity called blastocoel B inside. More specifically, cells with advanced cleavage form a thin layer on the surface of the embryo to form a blastocoel T, and the internal space surrounded by the ectoderm T forms a blastocoel B. In addition, an inner cell mass I in which a large number of cells are densely formed is formed in a part of the internal space.

Further, a zona pellucida ZP is formed so as to cover the outer surface of the vegetative ectoderm T. The zona pellucida ZP is a membrane mainly composed of glycoprotein and having a substantially uniform thickness. On the other hand, the vegetative ectoderm T formed by a large number of cells has various thicknesses depending on the position and is distributed so as to adhere to the entire inner surface of the zona pellucida ZP. As shown in FIG. 2 (b), the vegetative ectoderm T is formed by arranging a large number of cells C along the inner surface of the zona pellucida ZP. Due to the difference in the shape and size of the individual cells C, the vegetative ectoderm T becomes a layer having a different thickness depending on the position. In addition, in FIG. 2, the variation in the thickness of the vegetative ectoderm T is emphasized more than it actually is.

These structures that make up the embryo E, namely the zona pellucida ZP, the vegetative ectoderm T, and the inner cell mass I, have important significance as areas of interest in embryo evaluation. Therefore, the technique of automatically extracting the region corresponding to these structures from the captured image has great significance in supporting the evaluation work of the embryo by the user. However, in the tomographic image or the three-dimensional image obtained by OCT imaging, there is almost no difference in the luminance information between these structures. For this reason, it is difficult to accurately separate these structures by simple separation based on the difference in brightness.

As described above, a technique for automatically separating zona pellucida ZP, vegetative ectoderm T, and inner cell mass I from OCT images has not been established so far. In view of such a problem, Patent Document 2 previously disclosed by the applicant of the present application makes it possible to separate the region corresponding to the vegetative ectoderm and the inner cell mass from the OCT image. However, even with this technique, it cannot be said that the zona pellucida and the vegetative ectoderm are clearly separated. Specifically, the region extracted as vegetative ectoderm by this technique may include the region corresponding to the zona pellucida.

The first object of the image processing in the present embodiment is to separate the transparent band ZP and the feeding ectoderm T by analyzing the OCT image using the information obtained from the two-dimensional image obtained by the optical microscope imaging. Finally, in the three-dimensional image, the regions occupied by the transparent band ZP, the trophectoderm T, and the inner cell mass I are individually identified. The second purpose is to form the number of cells constituting the vegetative ectoderm T thus separated.

In a bright-field image or a phase-difference image obtained by imaging with an optical microscope, it is relatively easy to distinguish the zona pellucida ZP from the vegetative ectoderm T due to the difference in brightness and texture. And in well-cultured blastocyst stage embryos, the zona pellucida ZP has a generally uniform thickness. In other words, it can be considered that the zona pellucida ZP occupies the region of the embryo from its surface to a certain depth. From these facts, the thickness of the zona pellucida ZP was identified from the two-dimensional image data obtained by optical microscope imaging, and the zona pellucida ZP was identified from the three-dimensional image data obtained by OCT imaging when viewed from the surface of embryo E. By specifying the region to the depth corresponding to the thickness, the region corresponding to the zona pellucida ZP can be extracted separately from the trophozoite T. Further, for the separation of the inner cell mass I, for example, the method described in Patent Document 2 can be applied. With these, the first purpose can be achieved.

The vegetative ectoderm T thus separated from the other structures has irregularities due to the shape of the individual cells C, as shown in FIG. 2 (b). More specifically, protrusions due to the shape of individual cells C are formed on the surface of the vegetative ectoderm T facing the blastocoel B. Therefore, by detecting such a protrusion site and counting the number thereof, the number of cells C constituting the vegetative ectoderm T can be known. In this way, the second purpose is achieved.

FIG. 3 is a flowchart showing image processing in the present embodiment. This process is realized by the CPU 31 executing a control program prepared in advance and causing each part of the device to perform a predetermined operation. When the sample container 11 containing the embryo to be evaluated is set in the holding unit 10 (step S101), the optical microscope imaging and the OCT imaging by the imaging unit 20 are executed with the embryo as the imaging target.

In FIG. 3, steps S102 to S105 represent optical microscope imaging and data processing using the resulting two-dimensional image data. On the other hand, steps S111 and S112 represent OCT imaging and image processing using the three-dimensional image data obtained thereby. These two processes may be performed one after the other in chronological order, and in this case, it is arbitrary which process is performed first. Further, the processing time may be shortened by performing these processings in parallel.

First, the outline of the process of steps S102 to S105 will be described. In this process, the embryo is imaged with an optical microscope to acquire two-dimensional image data of the embryo, and the thickness of the zona pellucida is calculated based on the image data. Specifically, the microscope imaging unit 28 is positioned at the imaging position, the focal position is changed and set in multiple stages in the depth direction (Z direction), and imaging is performed each time. As a result, a plurality of image sets having different focal depths, so-called Z-stack images, are acquired (step S102).

A region corresponding to the zona pellucida ZP is extracted from each of these images (step S103). Then, one of the plurality of images that is most well focused on the zona pellucida ZP is selected (step S104), and the thickness of the zona pellucida ZP is calculated based on the selected images (step S105). These processes, that is, the process of extracting the region corresponding to the transparent band ZP (step S103), the process of selecting the most focused image (step S104), and the process of calculating the thickness of the transparent band ZP (step S105). Details will be described later. By the treatment up to this point, the thickness of the zona pellucida ZP in the embryo E to be evaluated is known.

On the other hand, in step S111, the object optical system 23 is positioned at the imaging position, and the embryo E is tomographically imaged while scanning the imaging position, so that three-dimensional image data of the embryo E is acquired. By binarizing the brightness of each part shown in the three-dimensional image data with a predetermined threshold value, the three-dimensional image of the embryo E is occupied by some structure having a relatively high density and becomes high brightness (hereinafter referred to as “high brightness”). , "Structural region") and a region with lower density and lower brightness. For example, medium M has low brightness in the OCT image.

As described above, the structure occupying the structural region may include zona pellucida ZP, vegetative ectoderm T, and inner cell mass I. In steps S113 and S114, they are separated from each other. Specifically, in step S113, the region corresponding to the zona pellucida ZP is separated from other structural regions. Then, in step S114, the vegetative ectoderm T and the inner cell mass I are separated.

By the treatments up to this point (steps S101 to S114), the first object of identifying the regions occupied by the zona pellucida ZP, the vegetative ectoderm T, and the inner cell mass I from the three-dimensional image of the embryo E is achieved. Subsequent treatments in steps S115-S118 are for the second purpose of counting the number of cells constituting the isolated vegetative ectoderm T. This process will be described in detail later.

In the following, each elemental technique for executing each step (steps S103 to S105, S111) of the above processing will be described in order. Since the process for acquiring the Z-stack image in step S102 is well known, the description thereof will be omitted. Further, as for the process of step S114, the technique described in Patent Document 2 can be applied, and thus the description thereof will be omitted.

In step S103, a region corresponding to the zona pellucida ZP is extracted from the two-dimensional image data obtained by imaging the embryo E with an optical microscope. This process can be performed using an appropriate image processing technique. For example, a pattern recognition technique for extracting a region having a specific feature from an image can be applied. Specifically, a classification model is constructed by supervised learning using the image of the zona pellucida acquired in advance as a teacher image, and the optical microscope image of the embryo E to be evaluated is divided into regions using this to obtain the image from the image. The area corresponding to the zona pellucida ZP can be extracted.

In this embodiment, a known semantic segmentation method is used as an example of the region division process. The semantic segmentation method is a technique for labeling each pixel in an image using a classification model pre-constructed by a deep learning algorithm. In this embodiment, this method can be used as follows.

First, an optical microscope image of an embryo in which the zona pellucida is captured with good image quality is prepared, and each pixel in the region corresponding to the zona pellucida in the image is labeled to that effect. Then, a classification model is constructed by executing deep learning using the original microscope image as input data and the label image as correct answer data. By giving an unknown image as input data to the classification model constructed in advance in this way, it is possible to obtain an output image in which the region corresponding to the transparent band of the image is labeled to that effect. By extracting such an area from the output image, it is possible to extract the transparent band as a result.

An example of a concrete method for constructing a classification model will be described. This process can be executed by various computer devices having a function of displaying an image and a function of accepting an operation input from a user. For example, it can be executed by the image processing device 1 or a general-purpose computer device such as a personal computer.

First, an optical microscope image of the embryo, pre-photographed in focus (in focus) on the zona pellucida, is displayed. In the image processing device 1, the display unit 352 can display the image. For the image displayed in this way, the instruction input from the user for designating the area corresponding to the in-focus transparent band is accepted. In this case, it is desirable that the user is an expert who has sufficient knowledge about the image of the embryo. Further, when the image processing device 1 is used, the teaching input can be received via the input device 351.

The area designated as the zona pellucida is given a label to that effect. By performing deep learning with the image labeled in this way as the correct answer data and the original image as the input data, a classification model for extracting the transparent band from the image is constructed. If necessary, a label other than the zona pellucida may be used.

This classification model is constructed by using an image in focus on the zona pellucida as an input image. Therefore, when this is applied to an unknown test image and the semantic segmentation method is executed, a region having a strong characteristic as a transparent band is extracted from the test image. Once the zona pellucida is extracted in this way, it is possible to determine its thickness. In order to obtain the thickness accurately, it is desirable to focus on the zona pellucida in the widest possible area in the image. That is, it can be said that it is desirable that the thickness of the zona pellucida is calculated based on an image in which the area of the region extracted by the semantic segmentation method is as large as possible. The area can be expressed by the number of pixels.

On the other hand, since the embryo E to be evaluated has a three-dimensional structure, the in-focus state to the zona pellucida ZP may not always be good in the image taken by appropriately determining the focal depth. Therefore, from the Z stack images acquired at different focal positions in the depth direction (Z direction), one image having the largest area of the region extracted as corresponding to the zona pellucida ZP is selected, and this image is used. Is used to calculate the thickness of the zona pellucida ZP.

In step S104, the image most focused on the zona pellucida ZP is selected from the Z stack images. In the semantic segmentation method in this embodiment, since the region of the zona pellucida in focus is extracted from the image, it is highly probable that the image having the largest area of this region is the image most focused on the zona pellucida ZP. It can be said that. Images that meet such conditions may be selected.

FIG. 4 is a diagram showing the relationship between the focal depth and the area of the extracted transparent band. In the upper part of FIG. 4A, an optical microscope image (Z stack image) captured at various focal positions in the depth direction and a region (corresponding to the zona pellucida ZP) extracted from these images by the semantic segmentation method. The relationship with the region) is schematically shown. Further, the graph at the bottom of FIG. 4A exemplifies the relationship between the focal position at the time of imaging and the area (number of pixels) of the extracted region. As shown in these, in the image Ia in which the Z-stack image is best focused and the zona pellucida ZP is clear, the area of the extracted region is maximized. In other words, the focal position corresponding to the image having the largest extracted area in the Z stack image can be regarded as the in-focus position. Therefore, this image Ia may be selected as the image that is the basis for calculating the thickness of the zona pellucida ZP.

However, for example, the image may be blurred due to vibration during imaging, and as a result, the apparent area of the zona pellucida may appear in the image larger than it actually is. FIG. 4B shows such a situation. For example, in the image Ib, it is assumed that the image is blurred due to vibration. Then, the area extracted as corresponding to the transparent band becomes apparently large, and this position may be erroneously determined as the in-focus position.

In this embodiment, this problem is solved by using the luminance difference between the pixels sandwiching the peripheral edge of the extracted region. That is, in a well-focused image, the boundary between the area corresponding to the zona pellucida and the surrounding area is clear, and therefore it is considered that there is a clear contrast between the brightness of those areas. On the other hand, in an out-of-focus image, the boundaries between these areas are not clear and therefore the contrast is not clear.

For this reason, instead of simply evaluating by the size of the area (number of pixels) of the extracted area, it is necessary to introduce an evaluation value that also reflects the magnitude of the brightness change in the edge portion of the area, that is, the sharpness. Therefore, it is considered that the above-mentioned erroneous determination can be reduced. There are various methods for quantifying the amount of such edge change, and they can be appropriately selected and applied.

In this embodiment, as an example thereof, a value obtained by multiplying the area of the extracted region by a coefficient reflecting the magnitude of the change in luminance at the edge portion is used as an evaluation value. As this coefficient, for example, a value obtained by squared the difference in luminance between pixels sandwiching an edge can be used. More specifically, the difference between the average value of the brightness of all the pixels adjacent to the edge in the extracted area and the average value of the activation of all the pixels adjacent to the edge outside the edge is obtained, and the difference is obtained. Is squared to obtain the above coefficient.

By doing so, as shown in the graph at the bottom of FIG. 4B, it is possible to reduce the risk that an increase in the extraction area due to the influence of the deviation in the image Ib causes an erroneous determination of the in-focus position. Although the luminance difference is squared here in order to make the coefficient a positive value, the absolute value of the luminance difference may be used as the coefficient instead.

The process for selecting an image focused on the transparent band can be configured as follows. This process corresponds to step S104 in FIG. For each image constituting the Z-stack image, a region corresponding to the transparent band is extracted in step S103. In order to obtain the area of the region extracted in this way, the number of pixels belonging to the region is counted for each image. Then, the average brightness of the adjacent pixels inside and outside the edge of the extracted region is obtained, and the difference is calculated. Based on these values, an evaluation value indicating the degree of focusing of each image is calculated. Specifically, the evaluation value is calculated by multiplying the number of pixels in the extraction region by a coefficient represented by the square of the brightness difference inside and outside the edge. The image having the maximum evaluation value obtained in this way is selected as the image most focused on the zona pellucida.

Next, the processing content of step S105, which calculates the thickness of the zona pellucida ZP from one optical microscope image selected from the Z stack image, will be described. For the selected image, the region corresponding to the zona pellucida ZP is extracted in step S103. In a well-cultured embryo E, it is considered that an annular region having a substantially constant width is extracted as a region corresponding to the zona pellucida ZP. In the image focused on the zona pellucida ZP, the width of the annulus, that is, the distance between the inner and outer edges of the annulus, corresponds to the thickness of the zona pellucida ZP.

Various methods can be considered for finding the width of the annulus, but a simple method is to use the Distance Transform function provided in the OpenCV (Open Source Computer Vision) library. By applying the Distance Transform function with one pixel on the inner edge of the annulus as the pixel of interest, the distance from that pixel to the nearest pixel on the outer edge of the annulus can be specified. This distance represents the width of the annulus at that position, that is, the thickness of the zona pellucida ZP. On the contrary, it is equivalent to use the pixel on the outer edge of the annulus as the pixel of interest and find the shortest distance from this to the inner edge. In this way, the average value of the widths obtained at each position on the circle can be used as the average value of the thickness of the transparent band ZP. In the following, the average value of this thickness will be represented by the symbol Ta.

The detailed processing contents of steps S103 to S105 have been described above. Next, the transparent band separation process in step S113 will be described. Here, in the three-dimensional image of embryo E obtained by OCT imaging, the range from the surface to the depth Ta is regarded as the zona pellucida ZP. Therefore, by extracting only this range from the three-dimensional image, only the structure corresponding to the zona pellucida ZP can be extracted. On the other hand, by erasing the structures in this range, the structures other than the zona pellucida ZP, that is, the vegetative ectoderm T and the inner cell mass I can be taken out.

For a three-dimensional image in which the region corresponding to the zona pellucida ZP is erased, as described in Patent Document 2, the remaining structure is subjected to region division processing using, for example, a Local Thickness operation. Thereby, the region corresponding to the vegetative ectoderm T and the region corresponding to the inner cell mass I can be separated. In this way, the regions occupied by the zona pellucida ZP, the vegetative ectoderm T, and the inner cell mass I can be individually extracted from the three-dimensional image of the embryo E.

Next, a process for achieving the second object of the present embodiment, that is, counting the number of cells constituting the vegetative ectoderm T (steps S115 to S119) will be described. The general flow of this process is to calculate the thickness distribution at each position of the vegetative ectoderm T, detect the peak position in the obtained thickness distribution, and assume that the number of significant peaks corresponds to individual cells. Is to count. Hereinafter, each step will be described in order.

In step S115, the three-dimensional image data represented by the XYZ Cartesian coordinate system is converted into a polar coordinate representation. Since the embryo E is generally spherical and the inside is hollow, the vegetative ectoderm T has a shape close to that of a spherical shell. In order to express such a structure more simply, it is preferable to express each position in polar coordinates (spherical coordinates) with the center of embryo E as the origin. Therefore, as shown in FIG. 5, coordinate conversion is performed from the XYZ Cartesian coordinate system to the rθφ polar coordinate system using the moving diameter r and the two declinations θ and φ as coordinate variables.

FIG. 5 is a diagram showing the correspondence between the XYZ Cartesian coordinate system and the polar coordinate system. As is well known, the origin O is common between the coordinates of the point P in the Cartesian coordinate system (x, y, z) and the coordinates of the point P in the spherical coordinate system (r, θ, φ). If there is, the following formula:
x = r ・ sinθ ・ cosφ
y = r ・ sinθ ・ sinφ
z = r · cosθ
There is a relationship.

Specifically, the center position of the embryo E is specified from the three-dimensional image data, and this center position 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 Cartesian coordinate system in the original signal data. Then, coordinate conversion from orthogonal coordinates to spherical coordinates is performed by an appropriate conversion process. By performing the coordinate transformation in this way, each position in the three-dimensional space specified by the XYZ Cartesian coordinate system can be expressed in polar coordinates.

The "center of the embryo" can be obtained, for example, as follows based on the three-dimensional image data. When the surface of the three-dimensional image of the embryo E represented by the three-dimensional image data can be regarded as a spherical surface, the center of gravity of the sphere in the image can be the center of the embryo. In three-dimensional image processing, a method for obtaining the center of gravity of a solid object is known, and such a method can be applied. Further, for example, a spherical surface or a spheroidal surface that approximately represents the surface of the embryo may be specified, and the center of the approximate surface may be the center of the embryo.

In step S116, the thickness distribution of the vegetative ectoderm T can be obtained from the three-dimensional image data thus transformed in polar coordinates. Specifically, the thickness distribution of the vegetative ectoderm T can be obtained by calculating the thickness of the vegetative ectoderm T in one radial direction in the polar coordinate space and performing this in various radial directions. can.

FIG. 6 is a diagram illustrating how to determine the thickness of vegetative ectoderm. In the example shown in FIG. 6A, considering one radius r1 extending from the origin O, the inner surface of the vegetative ectoderm T in the direction of the radius r1, that is, the surface Si facing the blastocoel B, and the outer surface, that is, the zona pellucida. The distance T1 from the surface So facing the ZP is defined as the thickness of the vegetative ectoderm T in the radial direction.

Further, in the example shown in FIG. 6B, attention is paid to the point P1 at which the radius r1 intersects the inner surface Si of the vegetative ectoderm T, and the shortest distance between the point P1 and the outer surface So of the vegetative ectoderm T is set. It is the thickness of the vegetative ectoderm T in the radius r1 direction. Alternatively, focusing on the point P2 where the radius r1 intersects the outer surface So of the vegetative ectoderm T, the shortest distance between the point P2 and the inner surface Si may be the thickness of the vegetative ectoderm T in the direction of the radius r1. good. For the calculation of the shortest distance, it is possible to use the Distance Transform function in the OpenCV library in the same manner as the above-mentioned zona pellucida ZP thickness calculation process.

FIG. 7 is a flowchart showing a process for determining the thickness of the vegetative ectoderm. First, one radial direction is selected (step S201). Specifically, one radial direction is specified by appropriately temporarily setting the values of the declinations θ and φ. Next, in the polar coordinate space, the position of the inner surface Si of the vegetative ectoderm T in the radial direction (the position of the point P1 shown in FIG. 6 (b)) or the position of the outer surface So (shown in FIG. 6 (b)). The pixel corresponding to the position of the point P2) is specified as the pixel of interest (step S202). For later processing, it is desirable to store the distance from the origin O to the pixel of interest in association with the radial direction specified by the set of declinations θ and φ. This information is expressed by the function R (θ, φ).

Subsequently, the distance between the inner surface Si and the outer surface So of the vegetative ectoderm T in the radial direction, that is, the thickness of the vegetative ectoderm T is determined (step S203). The thickness can be determined by the method shown in either FIG. 6 (a) or FIG. 6 (b).

The thickness of the vegetative ectoderm T thus obtained is associated with one radial direction and can therefore be expressed as a function of the declinations θ and φ. In the following, the thickness of the vegetative ectoderm T in one radial direction specified by the set of deviation angles θ and φ will be expressed as T (θ, φ). The obtained thickness T (θ, φ) is stored and stored in the memory 37 (step S204). By repeating the above treatment in various radial directions (step S205), a thickness distribution representing the thickness of the vegetative ectoderm T in each direction can be obtained.

In step S117, a map showing the spatial distribution of the thickness of the vegetative ectoderm T thus obtained is created. As the thickness distribution map, a two-dimensional map and a three-dimensional map can be considered. Further, as a two-dimensional map, a method of storing all directions in one map and a method of dividing the polar coordinate space into several spaces and storing each in one map can be considered. Hereinafter, specific examples of creating these mapping methods will be described.

FIG. 8 is a diagram showing an example of a two-dimensional mapping method. To map each direction expressed in polar coordinates on a two-dimensional plane, for example, a pseudo-cylindrical projection can be used. That is, as shown in FIG. 8A, the declination θ and φ correspond to latitude and longitude, respectively, and the same projection as when expressing the topography of the earth's surface, which is a substantially sphere, on a plane map is used. be able to. Several such drawing methods are known and can be appropriately selected according to the purpose.

Here, the Sinusoidal projection, which is a map projection that can maintain the area in an entity on a map by treating the argument θ as latitude and the argument φ as longitude, is used. At each point on the map specified by the combination of declinations θ and φ, pixels having a luminance value corresponding to the thickness T (θ, φ) of the vegetative ectoderm T calculated in each radial direction are arranged. Will be done. By expressing the thickness in each direction on the map in this way, the thickness distribution can be visualized. FIG. 8B shows a map showing the thickness of the vegetative ectoderm T at each point by luminance, using isoluminance lines instead of directly expressing by luminance. It can be said that the thickness in each direction is simply represented by contour lines. In addition, the thickness of the vegetative ectoderm T can be expressed by color or shade.

In FIG. 8B, the hatched upper right region represents the region where the inner cell mass I exists. The vegetative ectoderm T also exists between the inner cell mass I and the zona pellucida ZP, but the shape of the surface in contact with the inner cell mass I does not necessarily represent the shape of an individual cell. Therefore, in this embodiment, the region where the inner cell mass I exists is excluded from the derivation target of the thickness distribution.

As is a problem in the case of cartography, when the surface of a substantially sphere is represented by a two-dimensional map, all the information cannot be accurately represented. For example, in the case of the Sanson projection, which is a kind of equal area pseudo-cylindrical projection, the area of the ground surface can be accurately expressed, but the distance and direction cannot always be expressed correctly. In the case of the mapping of the present embodiment, the distortion becomes large especially at the peripheral edge of the map plane.

FIG. 9 is a diagram showing another example of the two-dimensional mapping method. In FIG. 9A, a three-dimensional image of the vegetative ectoderm T, which is represented as a substantially spherical surface in the polar coordinate space, is divided into four by two planes that pass through the origin O and are perpendicular to each other. Then, each of the divided sections Ta to Td is represented on one map. At this time, in the map corresponding to the section Ta, the intersection Ca of the straight line and the spherical surface passing through the center of the divided spherical surface, more specifically, the origin O and intersecting each of the two planes at an angle of 45 degrees. However, it should be roughly centered on the map. The same applies to the other compartments Tb, Tc, and Td.

In this way, as shown in FIG. 9B, the vegetative ectoderm T divided into four in the polar coordinate space is represented using only the central part of the map plane. Therefore, although the map has four faces, it is possible to accurately express the thickness distribution of the whole vegetative ectoderm T in each map without using the peripheral portion where the distortion becomes large.

In reality, each section is partially overlapped when considering the division of the spherical surface corresponding to the vegetative ectoderm T, and when mapping, as shown by a dotted line in FIG. 9 (b). It is more preferable that the thickness distribution is shown not only in the compartment Ta but also in a certain range outside the compartment Ta. Although divided for convenience of mapping, the actual vegetative ectoderm T is continuous to the outside of each compartment, and by mapping to a certain range outside the compartment, the thickness at the outer edge of each compartment This is because the continuity of the distribution can be expressed.

In this sense, instead of performing division, it is conceivable to create a plurality of two-dimensional maps with different central positions of the map and use only the central part thereof. Further, the number of divisions is not limited to 4 above and is arbitrary.

The above-mentioned two-dimensional mapping method is suitable for an application such as a bird's-eye view of the entire thickness distribution on a screen or a printed paper surface, but it is inevitable that the accuracy will decrease at the peripheral portion thereof. Therefore, the 3D map is superior in terms of accuracy. In a three-dimensional map, it is possible to ensure the same degree of accuracy in all directions in the polar coordinate space.

When displaying a three-dimensional map on the screen, for example, it is convenient for the user to temporarily display the map projected in one line-of-sight direction and change the line-of-sight direction according to the user operation. Can be improved. For example, if the animation display changes the line-of-sight direction in real time in conjunction with the user operation, the user can perform the evaluation work as if he / she is directly observing the embryo E in front of the eye.

As a three-dimensional mapping method for the ectoderm T, a method of giving information indicating the thickness distribution in each direction to an approximate spherical surface (or a spheroid) representing the approximate shape of the ectoderm T, and actual nutrition. There is a method of giving information on the thickness distribution to the curved surface corresponding to the three-dimensional shape of the ectoderm T. Further, in the latter case, there is an option of making the shape of the curved surface correspond to the inner surface Si of the vegetative ectoderm T or to correspond to the outer surface So.

The outer surface So of the vegetative ectoderm T is in contact with the zona pellucida ZP and therefore its surface shape is relatively smooth. In this sense, the curved surface corresponding to the outer surface So of the vegetative ectoderm T is substantially not much different from the spherical surface. On the other hand, the inner surface Si of the vegetative ectoderm T has irregularities corresponding to the shapes of the individual cells constituting the vegetative ectoderm T. Therefore, it is easier to imagine the actual uneven shape by mapping to the curved surface corresponding to the inner surface Si of the vegetative ectoderm T.

As a method of visualizing the thickness on the map surface, as in the case of two-dimensional mapping, a method of converting the thickness into luminance and displaying it, a method of expressing contour lines (contour lines), color coding, shading, etc. Is applicable. In the following, a case where mapping is performed on a spherical surface and a case where mapping is performed on a curved surface corresponding to the inner surface Si of the vegetative ectoderm T will be illustrated. Contour lines shall be used as the mapping method.

FIG. 10 is a diagram showing an example of a three-dimensional mapping method. FIG. 10A shows an example of mapping on a spherical surface, and FIG. 10B shows an example of mapping on a curved surface corresponding to the inner surface Si of the vegetative ectoderm T. Here, the three-dimensional map is shown in a state of being projected on a two-dimensional plane, but in reality, mapping is performed on a three-dimensional curved surface. In the case of display output on the screen, it is possible to display in various different line-of-sight directions. Also in these maps, the region indicated by hatching corresponds to the inner cell mass I.

For example, in the three-dimensional map of FIG. 10A, pixels of brightness corresponding to the thickness calculation result in the direction of the radius passing through the position are arranged at each position on the spherical surface centered on the origin O in the polar coordinate space. It can be created by doing. That is, pixels having a luminance value corresponding to the thickness T (θ, φ) may be arranged at the intersection of the radius represented by the two deviation angles (θ, φ) and the spherical surface.

The 3D map solves the distortion problem that occurs in the 2D map. In the two-dimensional map, apart from the problem of distortion, there is a problem that the length of the radius, that is, the information of the distance from the origin O is not reflected in the map. In particular, when the shape of the sample S deviates greatly from the spherical surface, the unevenness displayed on the map may not reflect the actual shape well. In the case of a three-dimensional map, it is possible to reduce such a deviation by, for example, using a spheroid that is closer to the shape of the sample as an approximate curved surface.

In particular, as shown in FIG. 10B, by mapping to a curved surface that reflects the actual shape, it is possible to create a map that is more realistic. The information required for mapping is the position (radius length) of the inner surface Si (or outer surface So) of the vegetative ectoderm T in each radial direction and the thickness of the vegetative ectoderm T in that direction. be. These pieces of information have already been obtained as R (θ, φ) and T (θ, φ), respectively, in the step of calculating the thickness of the vegetative ectoderm T (step S116 in FIG. 3). Therefore, by arranging pixels having a brightness value corresponding to the thickness T (θ, φ) at positions represented by coordinates (R (θ, φ), θ, φ) in the polar coordinate space, a three-dimensional map can be obtained. Is created.

When the three-dimensional map shown in FIG. 10B is projected and displayed on the two-dimensional display screen, it is necessary to clearly display the uneven shape of the curved surface of the map and the thickness of the vegetative ectoderm T at each position. .. In particular, the shape of a curved surface other than the outer peripheral portion of the projected two-dimensional figure may be difficult to read from the displayed image. In order to avoid this problem, it is necessary to express the shape of the curved surface and the thickness at each point in different ways. For example, the shape of the curved surface may be represented by contour lines, while the thickness of each portion may be represented by shading or color coding.

Returning to FIG. 3, a method of counting the number of cells constituting the vegetative ectoderm T from the thickness distribution thus mapped will be continued. The surface irregularities of the vegetative ectoderm T visualized as described above are due to the shape of individual cells. From this, for example, when the thickness is expressed by brightness, the region with higher brightness than the surroundings scattered on the map corresponds to the protrusion site corresponding to each cell in the vegetative ectoderm T. It is thought that it is doing. Therefore, by extracting such a region thicker than the surroundings, that is, a high-intensity region, it is possible to estimate where the individual cells are, and it is possible to count the number of cells constituting the vegetative ectoderm T. ..

From this, it is highly possible that the peak appearing in the thickness profile of the vegetative ectoderm T indicates the protrusion site corresponding to the individual cells constituting the vegetative ectoderm T. Therefore, in step S118, the protrusion portion is indirectly detected by finding a significant peak in the thickness profile.

FIG. 11 is a diagram showing the principle of peak detection processing. Further, FIG. 12 is a flowchart showing the peak detection process. As shown by the solid line in FIG. 11A, when the relationship (thickness profile) between the position in the direction and the brightness (that is, the thickness) is plotted along a certain direction on the map, the fluctuation of the brightness is changed for each position. There is, and a high peak appears with the brightness protruding from the surroundings. These are candidates for peaks that correspond to cell locations.

First, the region having higher brightness than the surroundings is extracted in this way. Specifically, the maximum value filtering process is executed for the thickness profile shown by the solid line (step S301). By appropriately setting the window size and executing the maximum value filtering process, the peak width of the profile is expanded as shown by the dotted line in FIG. 11 (a). Due to the maximum value processing, the luminance value increases at positions other than the peak, but the luminance value does not change at the peak position. Therefore, the position where the luminance value does not change (indicated by the broken line in the figure) is extracted as a peak candidate (step S302).

FIG. 11B shows the extracted peak candidates. If the brightness at the peak is extremely high or extremely low, the peak may be different from the detection target such as image noise or a peak that does not correspond to cells. Therefore, threshold values Lth1 and Lth2 are set for the luminance value, and among the peak candidates, those whose luminance value exceeds the threshold Lth1 which defines the upper limit value and those whose luminance value is lower than the threshold value Lth2 which defines the lower limit value are set. Erase (step S303). In FIG. 12, this process is described as “upper limit / lower limit process”. In FIG. 11B, the black circles indicate the peak candidates eliminated by this process. The white circles indicate the peak candidates that remain without being erased.

For the purpose of the present embodiment of identifying the individual cells constituting the vegetative ectoderm T, the thresholds Lth1 and Lth2 correspond to the upper limit and the lower limit of the size appropriate for the cells constituting the vegetative ectoderm T. Can be set. As a result, extremely large or small peak candidates that cannot be regarded as one cell can be excluded, and an error in the counting result of the number of cells can be reduced. Although the upper limit value and the lower limit value are set for the peak luminance value here, either one may be used. Further, if it does not cause a large error, the upper limit / lower limit processing may be omitted.

On the other hand, there is also a reasonable cell size in the direction along the plane of the vegetative ectoderm T. In other words, if the peak candidates are adequately representative of the individual cells, the distance between those peaks should be within a given range. Therefore, the distance between the remaining peak candidates that are adjacent to each other is calculated (step S304), and if there are peak candidates that are extremely close to each other, one of them is regarded as noise and eliminated (step S305). It is realistic to eliminate smaller peaks. FIG. 11C shows an example of the result of eliminating peak candidates whose distance to the adjacent peak is less than the preset minimum distance Dmin.

In this way, the peak indicated by the white circle in FIG. 11 (c) finally remains. It is considered that these peak positions represent the positions of the protrusions corresponding to individual cells. That is, by the treatment up to this point, the protrusion site existing in the vegetative ectoderm T surrounding the blastocoel B has been detected.

In addition, although the principle was explained here using a one-dimensional thickness profile to make the concept easier to understand, the actual thickness profile uses two-dimensional and three-dimensional maps in the case of using a two-dimensional map. In the case, it is three-dimensional. Therefore, in the case of using the two-dimensional map, the maximum value processing in step S301 is two-dimensional filtering. A circular filter is used in this case to apply a constant window size in each direction. The distance between peaks is also the distance on the two-dimensional map plane.

The two-dimensional map does not reflect the information on the length of the radius, and the distortion increases as it gets closer to the peripheral edge. Therefore, strictly speaking, it cannot be said that the ideal filtering process is executed in the circular filter window. However, for the purpose of detecting the peak corresponding to each cell position, it is practically possible to obtain sufficient accuracy. In particular, in the method of dividing an object into a plurality of parts and mapping as shown in FIG. 9, since the peripheral portion of the map where the distortion becomes large is not used, it is possible to sufficiently suppress the error due to such a factor. By making the range shown by the dotted line in FIG. 9B at least larger than the size of the filter window, it is possible to avoid an error due to the filtering process.

Also, if the thickness profile is shown in a 3D map, the filter window will be a 3D sphere. That is, the maximum value filtering process is spherical filtering. The distance between peaks is also the distance within the three-dimensional space.

In the three-dimensional map (FIG. 10A) in which the thickness of the vegetative ectoderm T at each position is projected onto the spherical surface, the information on the length of the radius from the origin O to the pixel of interest is not reflected in the processing. This has the effect of simplifying the processing, but can also be an error factor in the filtering processing. On the other hand, in the three-dimensional map (FIG. 10 (b)) projected onto the curved surface corresponding to the surface shape of the vegetative ectoderm T, the position of the pixel of interest reflects the distance from the origin O, so that attention is paid. By spherical filter processing centered on pixels, it is possible to eliminate errors and appropriately detect peaks.

FIG. 13 is a diagram showing an example of peak positions detected on a two-dimensional map. As shown by the dotted line in FIG. 13 (a), among the regions where the thickness of the vegetative ectoderm T locally protrudes from the surroundings, those corresponding to the above conditions are detected as significant peaks (black circles). These peaks are considered to correspond to each of the cells constituting the vegetative ectoderm T. FIG. 13B is a plot of only the detected significant peaks on a two-dimensional map. In this way, when significant peaks corresponding to individual cells are detected, the number of protrusions is determined by counting the number of them (step S119), and the cells indirectly constituting the vegetative ectoderm T. You can know the number of.

It is difficult to know the distribution of cells constituting the vegetative ectoderm T in the region occupied by the inner cell mass I, but by assuming that the cells are distributed at the same density as other regions, It is possible to estimate the number of cells in the embryo E as a whole. That is, the cell density was obtained by dividing the counted number of cells by the surface area of the vegetative epiderm T excluding the portion in contact with the inner cell mass I, and the surface area of the vegetative ectodermal T including the portion in contact with the inner cell mass I was obtained. By multiplying by the cell density, it is possible to estimate the total number of cells constituting the whole vegetative epiderm T.

FIG. 14 is a diagram showing an example of peak positions detected on a three-dimensional map. In the three-dimensional map as well, among the regions where the thickness of the vegetative ectoderm T locally protrudes from the surroundings, those corresponding to the conditions are marked as significant peaks (black circles) corresponding to the cells. Detected. By counting this over the entire map surface, it is possible to estimate the number of cells constituting the vegetative ectoderm T. The area occupied by the inner cell mass I can be handled in the same manner as in the case of the two-dimensional map.

As described above, when evaluating the cells in the layer covering the internal cavity, the protrusions thicker than the surroundings are detected in the thickness profile obtained from the three-dimensional image, and this is regarded as corresponding to each cell. It is possible to facilitate the quantitative evaluation of cells. When the outer shape of the object is substantially spherical, the processing can be facilitated by expressing the image data in the polar coordinate space.

The present invention is not limited to the above-described embodiment, and various modifications can be made other than those described above as long as the present invention is not deviated from the gist thereof. For example, the image processing apparatus 1 of the above embodiment has a function of OCT imaging and optical microscope imaging of the sample S, and a function of creating and outputting an output image from the imaging data. However, the image processing method of the present invention does not have an image pickup function by itself, and can be executed by a computer device that has acquired the image pickup data obtained by image pickup by another device having an image pickup function. In order to make this possible, the present invention may be implemented as a software program for causing a computer device to execute steps other than step S101 among the processing steps of FIG.

Such a program can be distributed, for example, by downloading it via a telecommunication line such as the Internet, or by distributing a computer-readable recording medium on which the program is recorded. .. Further, by loading this program into an existing OCT imaging device via an interface, it is possible to carry out the present invention by the device.

Further, in the above embodiment, various processes are executed after converting the three-dimensional image data into polar coordinate data, but the same protrusion portion as described above can be detected without performing such coordinate conversion. be. For example, when using the Distance Transform function for thickness detection, it is not always necessary to set the center (center of gravity) of the embryo, so it is possible to obtain the thickness directly from the three-dimensional image data of the trophic ectoderm represented by the Cartesian coordinate system. It is possible. Then, when the three-dimensional map shown in FIG. 10B is used, the map can be created in the XYZ Cartesian coordinate space as long as the positions of the inner surface or the outer surface are specified. Therefore, it can be said that conversion to a polar coordinate system is not necessary if the processing is a combination of these.

Further, for example, in the above embodiment, the zona pellucida and the vegetative ectoderm are separated from the three-dimensional image of the embryo E, and then the thickness distribution in the vegetative ectoderm is obtained. However, for example, if the thickness of the zona pellucida can be regarded as constant, even if the thickness distribution of the zona pellucida and the vegetative ectoderm is obtained without separating the zona pellucida, the change in thickness that appears in the zona pellucida is substantially the vegetative ectoderm. It reflects the thickness of.

Further, in the above embodiment, the present invention is applied to the evaluation of embryos in the blastocyst stage. However, the scope of application of the present invention is not limited to this, and it can be applied to the observation and evaluation of various cell masses having an inner cavity. For example, organoids of intestinal epithelial cells have cavities similar to those described above, the surface of which is covered with a layer formed by the cells. The present invention can also be suitably applied to the evaluation of such a cell layer.

As described above, as described by exemplifying a specific embodiment, in the image processing method according to the present invention, for example, the thickness of the cell mass along the moving diameter with one point inside the cell mass as the origin is determined. It may be configured as desired. When the cell mass has an outer shape close to a sphere or a spheroid, processing can be facilitated by obtaining the thickness in the radial direction in this way.

Further, for example, polar coordinate data representing the position of each point of the cell mass as a position in the polar coordinate space may be created from the three-dimensional image data, and the thickness distribution of the cell mass may be obtained based on the polar coordinate data. In particular, it is preferable that the center of gravity of the cell mass is the origin of the polar coordinate space. With such a configuration, it is possible to facilitate data processing for a cell mass having a shape close to a sphere or a spheroid.

In this case, for example, a cell mass thickness profile in which the radial direction in the polar coordinate space is associated with the cell mass thickness obtained in the radial direction is created, and the protrusion site is detected based on the profile. You may. With such a configuration, it is possible to easily detect a protrusion portion whose thickness protrudes from the surroundings from the change in the thickness of the cell mass indicated by the thickness profile.

Specifically, for example, the maximum value filtering process can be executed on the thickness profile, and the position where the value does not change before and after that can be set as the position of the protrusion portion. Since the thickness of the cell mass is larger than the surroundings at the position of the protrusion, the thickness value does not change even by the maximum value filtering. On the other hand, in the region where there is a thicker portion around, the value after the maximum value filtering changes due to the influence of the thickness. By utilizing this, it is possible to detect a portion having a thickness larger than that of the surroundings.

For example, it is possible to represent a thickness profile by a three-dimensional map in which pixels having a luminance value corresponding to the thickness of the cell mass at the position corresponding to the inner surface or the outer surface of the cell mass are arranged. Further, for example, by a three-dimensional map in which pixels having a brightness value corresponding to the thickness of the cell mass in the radial direction passing through the position are arranged at each position on the surface of the approximate sphere that approximates the cell mass to a sphere in the polar coordinate space. It is possible to represent a thickness profile. In these cases, the maximum value filter processing can be executed as the spatial filter processing for the spherical region centered on the pixel of interest.

On the other hand, a pixel having a brightness value corresponding to the thickness of a cell mass in one radial direction specified in the polar coordinate space by the two deviation angles on a coordinate plane whose coordinate axes are two deviation angles represented by polar coordinate data. It is also possible to represent the thickness profile by a two-dimensional map with. In this case, the maximum value filter processing can be executed as the two-dimensional filter processing for the circular region centered on the pixel of interest on the two-dimensional map. As a method of expressing the thickness distribution in a three-dimensional space by a two-dimensional map, for example, a pseudo-cylindrical projection can be used.

Further, for example, the polar coordinate space may be divided into a plurality of parts by a plane passing through the origin, and a two-dimensional map may be individually created and the protrusion portion may be detected for each of the divided spaces. When the shape of an object in a three-dimensional space is represented by a two-dimensional map, it is inevitable that the object on the map will be distorted. By dividing the space and creating the map individually, it is possible to avoid a large distortion at a specific location, and it is possible to improve the detection accuracy of the thickness distribution.

The image processing method according to the present invention may further include a step of measuring the number of detected protrusions. By requesting quantitative information in this way, it is possible to provide more useful information to the user. For example, when the cell mass of the object is a blastocyst stage embryo or an intestinal epithelial cell organoid, the present invention can be applied to the counting of cells constituting the surface thereof.

The present invention is suitable for supporting the observation and evaluation work of a cell mass having an inner cavity, for example, supporting the work of evaluating the state of a cultured embryo (fertilized egg), and assisted reproductive technology. It can be used for the purpose of selecting a good embryo with a higher pregnancy success rate in.

1 Image processing device 10 Holding unit 20 Imaging unit 30 Control unit 33 Signal processing unit 34 Imaging control unit E Embryo (cell mass)
S Sample Si Inner surface So Outer surface

Claims (14)

It is an image processing method that targets an image of a cell mass having a cavity inside.
A step of acquiring a three-dimensional image data representing a three-dimensional image of the cell mass by imaging a light interference tomography of the cell mass, and a step of acquiring the three-dimensional image data.
Based on the three-dimensional image data, a step of acquiring a thickness distribution representing the thickness of the cell mass for each position between the outer surface of the cell mass facing the outer space and the inner surface facing the cavity. When,
An image processing method comprising a step of detecting a protrusion portion having a thickness protruding from the surrounding region based on the acquisition result of the thickness distribution. The image processing method according to claim 1, wherein the thickness of the cell mass is obtained along a radius with a point inside the cell mass as the origin. From the three-dimensional image data, polar coordinate data representing the position of each point of the cell mass as a position in the polar coordinate space is created.
The image processing method according to claim 1 or 2, wherein the thickness distribution of the cell mass is obtained based on the polar coordinate data. The image processing method according to claim 3, wherein the center of gravity of the cell mass is the origin of the polar coordinate space. The image according to claim 3 or 4, wherein the protrusion is detected based on the thickness profile of the cell mass in which the radial direction in the polar coordinate space is associated with the thickness of the cell mass obtained in the radial direction. Processing method. The image processing method according to claim 5, wherein the maximum value filter processing is performed on the thickness profile, and the position where the value does not change before and after the thickness profile is set as the position of the protrusion portion. A three-dimensional map representing the thickness profile is created by arranging pixels having a brightness value corresponding to the thickness of the cell mass at the position corresponding to the inner surface or the outer surface.
The image processing method according to claim 6, wherein the maximum value filter processing is executed as a spatial filter processing for a spherical region centered on a pixel of interest. The thickness is obtained by arranging pixels having a brightness value corresponding to the thickness of the cell mass in the radial direction passing through the position at each position on the surface of the approximate sphere that approximates the cell mass to a sphere in the polar coordinate space. Create a 3D map that represents the profile
The image processing method according to claim 6, wherein the maximum value filter processing is executed as a spatial filter processing for a spherical region centered on a pixel of interest. A coordinate plane having two declinations represented by the polar coordinate data as coordinate axes has a brightness value corresponding to the thickness of the cell mass in one radial direction specified in the polar coordinate space by the two declinations. By arranging the pixels, a two-dimensional map representing the thickness profile is created.
The image processing method according to claim 6, wherein the maximum value filter processing is executed as a two-dimensional filter processing for a circular region centered on a pixel of interest on the two-dimensional map. The image processing method according to claim 9, wherein the two-dimensional map is created by a pseudo-cylindrical projection. The image processing method according to claim 9 or 10, wherein the polar coordinate space is divided into a plurality of parts by a plane passing through the origin, and the creation of the two-dimensional map and the detection of the protrusion portion are individually performed for each of the divided spaces. .. The image processing method according to any one of claims 1 to 11, further comprising a step of measuring the number of detected protrusions. Based on the three-dimensional image data representing the three-dimensional image of the cell mass obtained by photo-interference tomography of the cell mass having a cavity inside, the outer surface of the cell mass facing the outer space and the inner surface facing the cavity. And the step of acquiring the thickness distribution representing the thickness of the cell mass for each position between
A program for causing a computer to perform a step of detecting a protrusion portion having a thickness protruding from the surrounding area based on the acquisition result of the thickness distribution. A computer-readable recording medium on which the program according to claim 13 is recorded.

Images