An Image Processing Protocol to Extract Variables Predictive of Human Embryo Fitness for Assisted Reproduction

Abstract

Despite the use of new techniques on embryo selection and the presence of equipment on the market, such as EmbryoScope^® and Geri^®, which help in the evaluation of embryo quality, there is still a subjectivity between the embryologist’s classifications, which are subjected to inter- and intra-observer variability, therefore compromising the successful implantation of the embryo. Nonetheless, with the acquisition of images through the time-lapse system, it is possible to perform digital processing of these images, providing a better analysis of the embryo, in addition to enabling the automatic analysis of a large volume of information. An image processing protocol was developed using well-established techniques to segment the image of blastocysts and extract variables of interest. A total of 33 variables were automatically generated by digital image processing, each one representing a different aspect of the embryo and describing a different characteristic of the blastocyst. These variables can be categorized into texture, gray-level average, gray-level standard deviation, modal value, relations, and light level. The automated and directed steps of the proposed processing protocol exclude spurious results, except when image quality (e.g., focus) prevents correct segmentation. The image processing protocol can segment human blastocyst images and automatically extract 33 variables that describe quantitative aspects of the blastocyst’s regions, with potential utility in embryo selection for assisted reproductive technology (ART).

1. Introduction

Digital image processing and artificial intelligence (AI) techniques are rapidly gaining acceptance in medical sciences, including reproductive medicine [1,2], as tools for objective identification of lesions on medical images and even for disease diagnosis and prognosis. In the field of assisted reproductive technology (ART), these techniques are currently under investigation for identifying embryological parameters predictive of fitness for freezing and implantation [3]. Furthermore, automated image processing is now being used to improve the standardization and accuracy of diagnostics and prediction [4].

In the field of ART, computer-aided analyses are especially suitable for evaluation and classification of gametes and embryos [5], and, when combined with other methods, to predict embryo fate for ART [6]. Considering the complexity of each step and the vast number of variables to be considered, ART is an ideal field to test the potential of AI techniques [7]. Indeed, several algorithms have been developed to classify embryos or predict reproductive success that apply classical statistical methods such as logistic regression or various AI techniques [8].

The raw data for embryo evaluation is often acquired using a time-lapse system (TLS), which yields dynamic information on embryo morphology [9]. Digital image processing of these TLS-acquired image series has improved the accuracy of human embryo evaluation [10,11]. For instance, selection based on TLS has resulted in a 20% increase in pregnancy rate using ART [12]. However, embryo evaluation is still largely subjective and based on the embryologist’s experience, resulting in inconsistencies that contribute to the variability in clinical outcome [13,14]. Therefore, automation of embryo evaluation, although challenging due to the complexity of embryological variants [15], may lead to more consistent ART success.

Digital image processing as applied to ART includes multiple computational techniques that can extract mathematical variables from either conventional or TLS-acquired embryo images. Furthermore, AI techniques can automate this parameter selection process to obtain objective metrics associated with embryo morphology and quality [16,17]. Currently, the most widely used embryo classification system in Brazil and Argentina is that proposed by Gardner and Schoolcraft [18] and recently detailed by Puga-Torres et al. [19], which evaluates embryo fitness according to expansion of the blastocyst (EX), quality of the inner cell mass (ICM), and quality of the trophectoderm (TE) [20]. Image processing can evaluate embryonic quality by combining metrics for each region of the blastocyst (including the zona pellucida and the segmented ICM and TE) as described for bovine embryos [21]. In conjunction with AI tools, such image processing methods may help embryologists select the best oocyte or embryo for freezing, in addition to helping to distinguish euploid from aneuploid embryos [8,22].

Several previous studies have used computational techniques to extract various quantitative parameters describing the blastocyst stage. Santos Filho et al. [23] developed a method that yielded segmented images of the ICM and TE, while Singh et al. [24] applied identification and segmentation algorithms to isolate the TE. Matos et al. [25] even used digital image processing techniques for semiautomated selection of mouse blastocyst quality.

The inclusion of more objective embryo features, such as blastocyst width and/or area, may further improve selection for successful pregnancy [26]. With this study, the authors, based on their previous experience with the image segmentation of bovine and murine embryos [21,25], aimed to provide an objective and comprehensive assessment of human blastocyst morphology based on images captured by two distinct sources, EmbryoScope^® and Geri^® incubators, using various image processing techniques from the widely available Matlab^® platform. In the future, the feature variables extracted from image processing should be tested and validated in the assisted reproduction research (i.e., as program’s inputs to predict the human blastocyst fitness).

2. Materials and Methods

2.1. Study Design

The study was approved by the Institutional Review Board (IRB) of IVI Foundation, whose code corresponds to 1709-VLC-094-MM. In this research project, encrypted data was used, i.e., it was not possible to associate identified persons because the link with all information that identifies the subjects has been destroyed. Consequently, it was estimated that it would neither be necessary nor feasible to obtain the informed consent of the subjects, as it is a study with anonymized data. In addition, there was not express opposition from the subjects for the use of the data for the intended purposes.

2.2. Image Dataset

All raw images were of blastocysts collected from patients receiving assisted reproduction treatment between 2017 and 2019 at IVI-RMA (Valencia, Spain) and grown in either the EmbryoScope^® time-lapse incubator from Vitrolife or the Geri^® time-lapse incubator from Genea Biomedx (Figure 1).

The original dataset comprised 353 images of embryos growing in an EmbryoScope^® incubator and 474 images of embryos growing in a Geri^® incubator, of which 33 images from the EmbryoScope^® incubator and 165 from the Geri^® incubator were excluded. The 33 EmbryoScope^® images were excluded due to poor focus (n = 16), display of a nonblastocyst stage embryo (n = 14), incomplete visualization of the entire blastocyst (n = 2), or insufficient light to recognize individual regions (n = 1; Supplementary File, Figure S1). The 165 Geri^® images were excluded due to display of a nonblastocyst stage embryo (n = 79), blurred inner cell mass (n = 82), or display of a hatched blastocyst (n = 4). Therefore, 320 EmbryoScope^® and 309 Geri^® images were included (629 images in total). Of these, 289 images from the EmbryoScope^® incubator (90.3%) and 215 images from the Geri^® incubator (69.6%) were correctly processed by the methods proposed in the following sections, while the other 125 images were not correctly segmented, mainly due to the low quality of the raw image (Supplementary File, Figure S2).

2.3. Image Processing and Variable Extraction

The images provided by each incubator have unique characteristics (Figure 1) requiring different approaches for processing and feature extraction. The principal differences in processing relevant to image segmentation are presented in

Table 1. Different segmentation steps required for images provided by Geri^® and EmbryoScope^® incubators, and the variables extracted by each step.

Segmentation Step		Embryo Scope®	Geri®	Variable Extraction
1	Conversion of the input image to grayscale	X	X	-
2	Edge detection using a Canny filter	X	X	-
3	Partial embryo isolation	X	X	-
4	Histogram adjustment (contrast, stretching, and tone)	X	X	-
5	Prewitt adjustment	X	X	-
6	Binary conversion	X	X	-
7	Isolation of the whole blastocyst	X	X	Measurement of the area and radius of the blastocyst by equations using the Matlab® Image Processing ToolboxTM
8	TE and blastocoel + ICM isolation	X	X	Variables describing TE and the blastocoel + ICM using the Matlab® Image Processing ToolboxTM, local binary pattern (LBP) algorithm, and gray-level cooccurrence matrix (GLCM)
9	Threshold adjustment	X	X	-
10	Segmentation based on the grayscale		X	-
11	Determination of the binary distance		X	-
12	ICM isolation	X	X	Variables describing the ICM using the Matlab® Image Processing ToolboxTM, LBP algorithm, and GLCM

ICM: inner cell mass; TE: trophectoderm.

.

Segmentation of embryo images required between 10 and 12 steps to completely isolate the blastocyst, ICM, and TE. The steps used for segmentation of images from each incubator and typical changes to the image following each step are presented in Figure 2 and Figure 3.

2.4. Grayscale Conversion and Image Resizing

To reduce the computational load associated with image processing, raw images were first converted to grayscale [27,28] and adjusted to a standard size using the Matlab^® Image Processing Toolbox™ grayscale conversion and resizing functions [29].

2.5. Canny Filter Application for Edge Detection

In this work, a Canny filter (CF) was applied to detect image edges (Figure 4). In addition to a CF, Santos [30] added a Gaussian filter to reduce noise for improved location of edges. Therefore, the Gaussian filter described by Equation (1) was also applied, where σ_x and σ_y represent the standard deviations of pixels along the x and y axes, respectively.

$$G\left(x,y\right)=\frac{1}{2 \pi {\sigma }_x {\sigma }_y} e^{\frac{-x^2}{2{\sigma }_x^2}+\frac{-y^2}{2{\sigma }_y^2}}$$

(1)

This filter application is necessary to obtain information for isolating the region of interest (ROI), specifically the image of the blastocyst without surrounding background. A CF application can be found in the Matlab^® Image Processing Toolbox™ [29].

2.6. Detection of Circles by the Hough Transform

The Hough transform (HT) algorithm is used to detect forms such as lines, circles, and ellipses in digital images [31]. According to Seifozzakerini et al. [32], this technique is useful for the recognition of straight lines or circular forms, and thus is appropriate for ROI detection since embryos present a circular form. Circle parameterization described by Hough [33] is presented in Equation (2), where (a, b) represents the central coordinate and r represents the radius.

f_c (x, y) = (x − a)² + (y − b)² − r² = 0

(2)

These functions yield parameters such as the circle center and radius length that in turn allow for isolation of the ROI (Figure 5). Hough transformation was then followed by a final contrast adjustment to enhance ROI boundaries. The Matlab^® Image Processing Toolbox™ has all the functions needed to apply the HT as well as the contrast adjustment [29].

2.7. Histogram Adjustment

Contrast is an essential parameter determining image resolution [34]. To further improve image quality, contrast was optimized by histogram adjustment using the Matlab^® Image Processing Toolbox™ [29]. However, different contrast adjustments are required for images obtained by Geri^® (Figure 6) and EmbryoScope^® (Figure 7) to maintain an adequate standard for the desired variable acquisition.

2.8. Maxima of Gradient Magnitude

In a grayscale image, gradient magnitude (GM) is defined as the root mean square of the directional gradients along two orthogonal directions [35]. In general, GM maxima can be detected at the boundaries (edges) of a region [36]. For extracting the GM, several operators can be used, such as the Prewitt method, which is efficient for detecting edges. According to Yang et al. [37], this method can detect edges in the vertical (Gy) and horizontal (Gx) directions and calculate the GM using a pair of 3 × 3 convolution masks (Figure 8). This technique is applied using the Matlab^® Image Processing Toolbox™.

Using this technique, GMs are obtained (Figure 9b), which can later be converted into maximum values (Figure 9c). This process allows brightness-invariant areas (i.e., regions of the image with similar characteristics) to be captured and grouped. Based on the maximum magnitude values, pixels within a specific range can be extracted.

2.9. Binary Image Construction

A binary image is composed of a matrix containing only 0s and 1s, where the 1s represent the object of interest [38]. To obtain a binarized image, pixel values within a defined range are converted to 1s, while all other pixels are set to 0 (Figure 10b). This process is important for applying the HT and allows total isolation of the blastocyst (Figure 10c), TE (Figure 10d), and blastocoel plus ICM (Figure 10e).

Segmentation of the image into blastocyst, TE, and blastocoel plus ICM allows the extraction of several variables. In this process, variation and similarity of texture, uniformity in gray-level distribution, proximity of the grayscale to the GLCM diagonal in ICM and blastocoel, mean gray level, and standard deviation of gray level are extracted. In addition, brightest and darkest regions, average brightness, most frequent light intensity (modal value), binary image sum, blastocyst radius, and area measurements are calculated.

The textures were examined by the local binary pattern (LBP) algorithm and by constructing a gray-level cooccurrence matrix (GLCM). The LBP algorithm is a widely applied approach to recognize and analyze texture in images [39]. In this study, the LBP technique was applied to determine the textures of the ICM, TE, and EX. In Matlab^®, the LBP algorithm is included in the Computer Vision Toolbox^TM [40]. Feature extraction using the LBP algorithm consists of two main steps (Figure 11), thresholding and encoding. In the first step, all neighboring pixels in each defined area (e.g., pattern) are compared to the value of the central pixel. Then, all values of the neighboring pixels are changed to a binary value (0 if lower than the value of the central pixel or 1 if larger). After this transformation, the encoding step involves conversion of the binary number obtained to a decimal value [39].

The GLCM algorithm was also applied to examine the textures of the ICM, TE, and the blastocoel plus ICM images obtained by segmentation. This statistical method, included in the Image Processing Toolbox^TM, essentially examines texture by analyzing the spatial relationships among pixels, specifically how often pairs of pixels with specific values and in a specified spatial relationship occur in an image. This analysis yields the GLCM, from which statistical measures are extracted [29]. Using this algorithm, it is possible to extract information such as variation of the texture level, texture similarities, uniformity of the gray-level distribution, and proximity of the grayscale to the GLCM diagonal.

The brightness and darkness variables of the image were then extracted from the TE and blastocoel + ICM images. The darkness of the image was measured as the sum of the top 10% darkest pixels, and the brightness as the sum of the remaining pixels. The mean brightness of each region was estimated using the pixels with values between the 10% darkest and the 10% brightest. Image metrics such as radius and area were extracted from the distribution of pixels obtained by the Matlab^® Image Processing Toolbox^TM [29].

2.10. Grayscale Intensification

Grayscale standardization alone cannot provide good texture segmentation results. In addition, it is important to intensify the light and dark tones of the image using an optimal threshold (ψ) specific for Geri^® and EmbryoScope^® incubators. Using Equation (3), it is possible to convert the values of a pixel matrix ($I_{ij}$) and increase the frequency of maxima and minima, yielding a new higher contrast matrix ($G{I}_{ij}$) as depicted in Figure 12.

$$G{I}_{ij}={\left(\frac{I_{ij}}{\psi }\right)}^2\ast \psi$$

(3)

• $I_{ij}$: grayscale pixel matrix.
• $\psi$: threshold.
• $G{I}_{ij}$: matrix containing the converted values of $I_{ij}$.

2.11. ICM Partial Isolation Based on the Grayscale

Segmentation based on grayscale is a step specific to images from the Geri^® incubator, which have greater contrast between light and dark tones. In this process, the weight of each pixel in the image is calculated as the absolute difference between the pixel value ($G_{ij}$) and a previously defined threshold ($\psi$) (Equation (4)). Using these calculated weights, it is possible to generate a new image ($I$) containing the partially segmented ICM. Resources from Matlab^® Image Processing Toolbox™ [29] were used for this process. The method is illustrated in Figure 13a, and the results of partial ICM segmentation are shown in Figure 13b.

$$I_{ij}=|G_{ij}-\psi |$$

(4)

2.12. Determination of Binary Distance

Another technique used only for Geri^® images is binary Euclidean distance ($BDij$) calculation (Equation (5)). The binary distance is defined as the distance from every pixel to the nearest nonzero pixel [41]. After binary transformation of the partially segmented ICM, the inverse of this image is submitted to the distance calculation process (Figure 14). This improves the quality of the final segmentation by eliminating noise present in the image.

$$BDij=\sqrt{{\left(x_1-x_2\right)}^2+{\left(y_1-y_2\right)}^2}$$

(5)

2.13. ICM Isolation by Gabor Filter

Segmentation by Gabor multichannel filters is an established method in image processing. This technique allows characterization and separation of textures by creating channels for different frequencies and orientations, similar to the human visual system. To segment the ICM and obtain the Gabor magnitude, we initially used the combinations of frequency and orientation presented by Jain and Farrokhnia [42]. After this process, Gaussian low-pass filters were applied (Figure 15e and Figure 16b), an important postprocessing step to normalize the different outputs and decrease the magnitude variations according to Kim and Kang [43]. Gaussian low-pass filters enhance smoothness and reduce noise (the difference between neighboring gray levels). In this way, it is possible to differentiate image regions in different planes (Figure 15f and Figure 16c) and fully segment the ICM (Figure 15g and Figure 16d). In the last steps, clustering is performed using the Matlab^® Statistics and Machine Learning Toolbox™ [44].

Additional variables were extracted by isolation of the ICM. The textures were then analyzed by the LBP algorithm and the GLCM. The area was calculated by the distribution of the pixels in this segmentation. The blastocoel area was calculated by subtraction of the area of the ICM (Figure 15g and Figure 16d) from the total area of the blastocoel + ICM image (Figure 15a and Figure 16a). Expansion of the embryo was estimated indirectly by calculating the area of embryo (without the TE) occupied by the ICM.

2.14. Evaluation of the Image Segmentation

The effectiveness of the segmentation (in correct or incorrect) was based on the visual screening by an experienced embryologist (mainly MFGN; 21, 25, 45, 46). This work did not use independent experts or human-annotated segmentation as comparison standards.

3. Results

The image processing steps summarized in Table 1 automatically generated a total of 33 variables, each representing a different feature of the EX or a subregion (ICM or TE). These variables can be divided into three main categories, EX, ICM, and TE, according to the ROI represented (Table 1).

These variables can also be divided into six categories according to the image property represented: texture, gray-level average; gray-level standard deviation, modal value (most frequent pixel value), relations, and light level.

Texture variables denote repeating random regular patterns that provide measures of structural arrangements on surfaces. They represent the different interactions among pixels, from the differences in gray level from pixel to pixel in local regions of the image to the spatial arrangement of gray levels throughout the image. Gray-level average, gray-level standard deviation, and modal value represent the overall brightness/darkness and variation in brightness of the embryo. Relations refer to associations between otherwise distinct features, such as radius and area of the blastocyst. Finally, the light level variables represent the brightness variation in different regions of the embryo (Supplementary File, Figures S3–S5).

Variable Definitions

Texture

1. 

Texture variation in ICM and blastocoel

Represents the texture gray-level differences in ICM and blastocoel. Reflects image sharpness and depth of grooves. Deeper grooves are related to greater variation and sharpness.

2. 

Texture similarities in ICM and blastocoel

Represents the linear dependency of gray levels.

3. 

Uniformity of the gray-level distribution in ICM and blastocoel

Represents the uniformity of the gray-level distribution through the sum of GLCM elements squared.

4. 

Proximity of the grayscale to the GLCM diagonal in ICM and blastocoel

Represents the distribution of elements in the GLCM according to the diagonal reference values.

5. 

Texture variation in ICM

6. 

Texture similarities in ICM

7. 

Uniformity of gray-level distribution in ICM

8. 

Proximity of the grayscale to the GLCM diagonal in ICM

These variables correspond to variables 1, 2, 3, and 4, respectively, but are calculated from the ICM segmented image.

9. 

Texture variation in TE

10. 

Texture similarities in TE

11. 

Uniformity of the gray-level distribution in TE

12. 

Proximity of the grayscale to the GLCM diagonal in TE

These variables correspond to variables 1, 2, 3, and 4, respectively, but are calculated from the TE segmented image.

13. 

Local texture descriptor in EX

Texture pattern descriptor used to describe the local texture patterns of an image. The algorithm extracts information about the local texture by establishing a threshold for “n” neighbors relative to the value of the central pixel. The result is a binary number between 0 and 255.

14. 

Local texture descriptor in ICM

15. 

Local texture descriptor in TE

These variables correspond to variable 13, but are calculated from ICM and TE segmented images, respectively.

Gray-level average

16. 

Gray-level average in ICM and blastocoel

Represents the average pixel intensity (gray-level value) within the segmented ICM plus blastocoel image.

17. 

Gray-level average in TE

Represents the average pixel intensity (gray-level value) within the segmented TE.

18. 

Blastocyst gray-level average

Represents the average pixel intensity (gray-level value) for the entire blastocyst with zona pellucida.

Gray-level standard deviation

19. 

Gray-level standard deviation in ICM and blastocoel

Represents the standard deviation of pixel intensity within the segmented ICM and blastocoel.

20. 

Gray-level standard deviation in TE

Represents the standard deviation of pixel intensity within the segmented TE.

Modal value

21. 

Modal value in ICM and blastocoel

Represents the most frequent light intensity value in the segmented ICM and blastocoel.

22. 

Modal value in TE

Represents the most frequent light intensity value in the segmented TE.

Relations

23. 

Blastocyst sum of binary image

The image threshold is given by the sum of the binary image divided by the area of the isolated blastocyst calculated from the radius.

24. 

Blastocyst radius

Radius of the blastocyst.

25. 

ICM area

Area of the segmented ICM.

26. 

Blastocoel area

Area of the segmented blastocoel.

27. 

Ratio between ICM and blastocoel

Represents the expansion of the blastocyst.

Light level

28. 

Mean luminosity in ICM and blastocoel

All pixels with a luminous intensity between 10 greater than and 10 less than the average intensity are counted, and this value is then divided by the total area of the blastocyst.

29. 

Mean luminosity in TE

Calculated as described for variable 28 but using values from the segregated TE.

30. 

Brightest region in ICM and blastocoel

Brightest area in the segmented ICM and blastocoel.

31. 

Brightest region in TE

Brightest area in the segmented TE.

32. 

Darkest region in ICM and blastocoel

Darkest area in the segmented ICM and blastocoel.

33. 

Darkest region in TE

Darkest area in the segmented TE.

While 33 variables were extracted, it is unclear if all can be modulated independently. Collinearity refers to a strong association between two independent variables, and multicollinearity refers to a strong association among more than two independent variables. Multicollinearity alters the grouping because the collinear variables end up implicitly receiving a greater weight. A collinearity analysis can be performed to identify if any of these 33 variables are related [45].

4. Discussion

In this work, we present a sequence of digital processing steps to segment human embryo images obtained from two distinct sources. The method produced 33 mathematical variables from the whole embryo and its specific regions (ICM, TE, and blastocoel). We speculate that these variables may help in identifying embryos most suitable for freezing and ultimate implantation for ART.

This sequence of steps was adapted from previous studies by our group on murine embryos [25] and bovine embryos [21,45,46]. The relatively high success rate of segmentation, especially for embryos photographed in the EmbryoScope^® incubator (90.3%), suggests that these previous findings were beneficial for analysis human of embryo images, a more challenging application due to both greater variability in embryo morphology compared to mice and cattle [47] and the more rigorous standards that must be met for ART.

In addition to work on nonhuman species, preliminary results have been published on digital processing of human blastocyst images from EmbryoScope^® for implantation selection [48,49]. In contrast, Geri^® images have been used much less extensively for AI development [50]. This paucity of research may reflect less ideal imaging conditions. Indeed, the rate of successful segmentation was much lower using Geri^® images (69.6%) and required extra steps.

Several of these proposed steps, such as application of the HT algorithm, have already been applied in previous image processing studies. Van den Heuvel et al. [51] used the HT to detect the center of the fetal skull and a dynamic programming algorithm to detect the outside of the fetal skull. They also fitted an ellipse to the dynamic programming result to measure head circumference, estimate gestational age, and monitor fetal growth. Huang et al. [52] also used the GLCM approach to identify textural features on ovarian adenocarcinoma cells indicative of chemoresistance. Specifically, the GLCM was used to calculate contrast, energy, entropy, and homogeneity, variables that collectively can reveal the disordered surface morphology characteristic of cancer cells. Alternatively, Nanni et al. [53] used the LBP algorithm to derive texture descriptors for classifying 2D HeLa images, cells in pap smear datasets, and pain levels from the facial images of newborns. Thus, our chosen algorithms constitute a set of mathematical tools already used for a wide range of image processing applications.

Some of the variables extracted may have little or no biological significance despite conveying information on texture, intensity of edges, and other features. However, these quantifiable features may still be advantageous over the subjective imaging criteria used to grade embryos [18,54]. As proof of concept for the potential utility of these variables in embryo selection, we are performing a pilot study (in progress) on prediction of pregnancy and aneuploidy. Briefly, the same 33 extracted variables were used as inputs for artificial neural network (ANN)-based programs to predict aneuploidy or pregnancy (fetal heartbeat) from blastocyst images. The ANNs were chosen based on a genetic algorithm search for the fittest model, similar to that described by Rocha et al. [21,45]. We recently obtained 96% predictive accuracy for the presence of a fetal heartbeat based on processing of 172 images from Embryoscope^® and receiver operating characteristic curve analysis as well 95% accuracy for prediction of aneuploidy based on 277 images using the same analytic processes. Fine-tuning the image processing protocol by creating high-resolution mosaic images of the whole embryo in segmentation masks to provide additional information on certain structures (e.g., ICM or TE) or features (texture, contrast, brightness, or edges) inscrutable by visual inspection may further enhance predictive efficacy.

In summary, we describe an image processing protocol that can successfully segment human blastocyst images from two distinct sources and extract 33 variables with potential utility in embryo selection for ART.