Objective assessment of bull sperm motility parameters using computer vision algorithms

The quality control of frozen semen samples from cattle is established by parameters such as percentage of progressive motility (% MP) because it is related to the fertilization capacity of bulls. Nowadays, sperm quality test is performed by direct visual inspection of sperm through a microscope in andrology laboratories. However, there is a high subjectivity in the observation and assessment depending on the observer; thus, causing unreliable diagnoses and non-repetitiveness in results. The development of a low cost computer tool was proposed to identify the individual sperm motility in cattle through the application of artificial vision algorithms. The methodology consisted of: the acquisition and pre-processing of videos obtained from thawed cryopreserved samples, the segmentation, filtering and detection of spermatozoa using the Fisher Discriminant Analysis and Adaptive Gaussian Models, followed by the assignment and construction of sperm trajectories through the Munkres Algorithm and Kalman Filter. Finally, the characterization and assessment of sperm motility parameters were performed based on the criteria present in computer-aided semen analysis (CASA) commercial systems. The results obtained showed high correlation for the individual sperm motility with a determination coefficient of 0.8143 for 10 different samples of bull sperm with respect to manual analysis. Likewise, a concordance coefficient of 0.966 was found in the 95% confidence interval using the Bland-Altman test, indicating that the measures were highly similar. In this way, the methodology is a reliable technological support that contributes to the improvement of quality control of semen samples from cattle.


INTRODUCTION
In Colombia, agricultural activity and majority of the bovine livestock have a great impact on the internal economy of the national territory (DANE, 2016).The country has a large cattle market and, therefore, a demanding commercialization of semen straws of cattle.One of the most important aspects is the sperm quality because there is a high demand of bovine semen straws for artificial insemination processes (Livestock Context, *Corresponding author.E-mail: eeduardoroa@ucundinamarca.edu.co.Tel: 314 5482 315.Fax: (+57 1) 828 1483.Ext: 142.
Author(s) agree that this article remains permanently open access under the terms of the Creative Commons Attribution License 4.0 International License

2017).
Nowadays, semen analysis is usually done manually through physical observation of spermatozoa under the microscope by specialists in andrology for a long time period (World Health Organization, 2010).However, this conventional procedure can have errors due to factors such as: subjectivity in the determination of progressive and non-progressive moving spermatozoa, lack of repetitiveness in the results and ocular fatigue of the expert for long working hours, thus, affecting the result of the analysis and the objective determination of the bull fertility (Hidalgo et al., 2017).On the other hand, commercial systems of computer-aided semen analysis (CASA) (Lu et al., 2014) have been developed, which allow the identification of seminal parameters with automatically high precision.However, these systems are closed and expensive, hence their implementation is not common in the livestock sector.Although in Colombia, only the case of one commercial system is reported, researchers have already developed algorithms to determine sperm parameters and its relation with semen quality, as well as detection of velocities of cells by CASA (Nagy et al., 2015).
The main objective of this research was to develop a methodology that reduces errors in the results of sperm motility analysis, by implementing a computational tool that performs artificial vision techniques on recorded videos, to provide the specialist, a support technology based on results obtained objectively and makes it to be easily adaptable to many different laboratories.
The methodology implemented consists of four phases in general; phase 1 is the acquisition and pre-processing of the videos through a Leica DM500 led binocular microscope, phase 2 is the segmentation and extraction of motile and non-motile spermatozoa.Subsequently, in phase 3, the algorithm was implemented to assign detections and construct trajectories described by the spermatozoa in the video.Finally, the proposed method analyzes these trajectories of the spermatozoa and identifies the sperm motility parameters following the protocols established by the OMS manual.
In the present investigation, there is a report on the development results, implementation and validation of the computational tool for individual sperm motility assessment in bovines from videos of semen samples obtained from 10 specimens.

MATERIALS AND METHODS
Initially, preparation of the seminal samples was performed according to established protocols in andrology and veterinary medicine laboratories (Vincent et al., 2012).Considering that, this procedure can influence the results of the motility analysis (Contri et al., 2010).In the investigation, the authors only worked on samples of frozen semen for each breed of bull (white black ear), the frozen samples were withdrawn from the liquid nitrogen vessel, each straw was thawed by a water bath for 60 s and 5 μL of semen were placed on a preheated slide at 37°C.Subsequently, the videos were captured using a Leica DM500 led binocular microscope, with digital The proposed methodology for the identification of individual sperm motility consists of 4 phases, as shown in Figure 1.In the first phase, video acquisition and pre-processing of the video were performed.In the second phase, the segmentation and detection of spermatozoa in motion were performed through adaptive Gaussian models.Subsequently, in the third stage, the allocation and construction of trajectories were performed.Finally, validation of the computational tool was performed through sperm motility indexes.Each of the four phases implemented are described below.

Acquisition and pre-processing
When the sample with 5 μL of semen was prepared on a slide and covered with a coverslip, it was focused with a 40x eyepiece using the Leica DM500 led binocular microscope made in Germany and a laptop with Intel core i5 processor and 4 GB memory ram, respectively.Afterwards, the camera software obtained the videos that were later stored in a computer, with the following attributes.Format: WMV; resolution: 640 x 480 pixels; frame rate: 30 Hz; data rate: 1572 kbp; duration: 5 s.
After obtaining the videos, a pre-processing stage was applied which allowed adjusting the contrast of the frame in order to increase the difference between light and dark pixels, thus improving the performance of the segmentation process.

Segmentation and detection of spermatozoa
For extraction of the spermatozoa, a similar technique of frame differentiation was used as developed by Stauffer and Grimson (1999), called adaptive background Gaussian models, with comparison of the value of each pixel with respect to Gaussian functions that are modeled adaptively.It is important to define the probability of each pixel belonging to a background image region.The variance of the pixel values over time allows the modeling of Gaussian functions present in an image.Therefore, these models make prediction possible when a particular pixel is part of a moving object or the background of the image.The probability of observing the pixel value is obtained from Equation 1: Where, K is the number of Gaussian models, ω is an estimated weight coefficient for a model i per unit of time t.The variable η represents the conventional Gaussian probability density, where takes the value of the pixels in the image.μ is the mean value and Σ is the covariance existing in the combination of models.Likewise, by Equation 2, Gaussian probability density was determined adaptively, allowing differentiation of objects that have motion on a video.
( ) Where, n is the domain size of the normal distribution.
Since the resulting binary image exhibits in most cases, unwanted particles and objects due to the movement of the seminal sample, it is necessary to carry out a filtering step in order to optimize the image and prevent further erroneous detections in the following process stages (Maintz, 2005).Later, an analysis of connected components was implemented in order to recognize the regions of the binary image and to develop an analysis of their geometric parameters (Gonzalez and Woods, 2008), in which fundamental data such as area and center of mass (centroid) are abstracted.The centroid denotes the coordinate on the image which is located in the region of pixels of interest, which in this case, would correspond to the head of a spermatozoon.
Considering that the adaptive Gaussian modeling method can only segment spermatozoa in motion, the linear discriminant analysis method described by Bishop (2016) aims to segment all the spermatozoa present in the frames, regardless of whether they are in motion or not, in order to determine the total motility of the sample.This analysis allows the classification of data between two or more classes.Its function is to maximize inter-class variance and minimize intra-class variance.In this way, a discriminant vector is created to maximize the variance between the pixels belonging to the bottom of the head of a spermatozoon by Equations 3 and 4.
Where, μ is the mean of all pixels, μc is the mean of each class, NC is the number of patterns in the class and xi is the value of each pixel.Subsequently, the discriminant vector is obtained with Equation 5 that separates the two classes and allows segmentation.

Assignment and construction of trajectories
When applying the segmentation algorithm to the video, multiple sperm detections are evident in each frame.These must be interpreted and organized under a set of real trajectories, in order to ensure that different detections correspond to a common trajectory along the frames and, in this way, registering a record of the coordinates of each moving target.
The assignment and construction of the tracks followed by the spermatozoa are developed at first through the application of Munkres assignment algorithm, followed by a Kalman filter.The first technique corresponds to variation of the assignment algorithm of Munkres (1957).This mathematical model aims to make a direct assignment of m elements to n destinations by means of an m x n matrix called the cost matrix.These costs are values related to the probability of assigning an element of row m to one of column n whose total cost is minimal.
Where m represents the detections in the binary image, n represents the sperm trajectories and each element of the matrix relates to the difference between estimated position and the actual position of the sperm.On the other hand, the Kalman filter is a set of mathematical equations that provide efficient computational means to estimate the state of a process, in such a way that it minimizes the mean square error.This filter allows estimations of states in past, present and even future.This algorithm was designed by Kalman (1960) which describes its operation as a predictive and corrective cycle that communicates with each other through a feedback.In Equations 6 and 7, the Kalman filter prediction stage is described in discrete time: Where, ̂ is the current state prediction, ̂ represents the previous state prediction, is the position measurement, A and B are state transfer matrices, is the error covariance a priori and Q is the process noise covariance.Subsequently, the position is updated through the correction stage described by Equations 8, 9 and 10, which consist of obtaining new measurements of the state variable and its conjugation with the a priori estimation of the prediction.

(
) Where, is the Kalman gain, H is the measurement state matrix, R is the covariance of the measurement error and is the current system measurement.
In this way, the positions that define the trajectories of the spermatozoa are determined, considering both the observed measurements and the coordinates estimation.All previous measurements constitute a criterion for predicting future measurements.The correction, in turn, updates the data and incorporates it into the information history through a feedback loop.In the cost matrix, the difference between the current position measurement and the Kalman prediction is stored, and eventually, a given detection is assigned with a greater probability of belonging to a real trajectory, as shown in the Figure 2.
The first detections presented in the initial frame of a video are classified in initial positions of trajectory and the unassigned detections become a new trajectory.In the case of the assigned detections, the update of the measurement in each frame is given by assigning a new element to the track.On the other hand, those tracks that do not comply with a minimum number of points are automatically deleted by the algorithm, because they are insignificant or indicate a detection error.

Characterization of sperm motility
Finally, the trajectory of the spermatozoon is defined by means of a matrix that represents the coordinates X and Y of the positions through time within the image, which allows determination of the most common parameters of sperm motility as realized in the systems CASA (Agarwal and Sharma, 2017) which have: curvilinear velocity (VCL), average velocity (VAP), straight line velocity (VSL), linearity index (LIN), straightness index (STR), wobbling index (WOB), percentage of total motility (% MT) and progressive motility (% MP).
The curvilinear velocity is a fundamental parameter in individual sperm motility.The other parameters are calculated from this value.The total distance was obtained by summing the distances between each coordinate k in the space of the image and its multiplication by a factor (A) that converts Cartesian coordinates ( ) into micrometers dominion (μm).The velocity was determined in Equation 11 by multiplying the total distance by the capture frequency of the video (s -1 ) as follows: Because the average velocity corresponds to the smoothed real trajectory, a mean filter was applied for each column of the coordinate matrix and subsequently computed.Mathematically, the sequential mean was calculated for each element of the column vector ( ), obtaining another average column vector ( ) as shown in Equations 12 and 13.Then, the same equation previously used was applied, but adapted to the ne values (Equation 14): Also, the rectilinear velocity was computed from the initial point ( ) to the final point ( ) of the trajectory described by the spermatozoon, and likewise using Equation 15, the velocity was converted to units of µm/s: Subsequently, the linearity indices were calculated in Equation 16, straightness in Equation 17and wobbling in Equation 18: Moreover, a characterization stage was designed in order to classify each spermatozoid into a motility category depending on its speed and trajectory.Finally, the percentages of total motility and progressive motility are determined by Equations 19 and 20, being the progressive motility percentage and one of the most important descriptors of the quality of a seminal sample as described in Vincent et al. (2012).
The criterion implemented for the classification of sperm motility in the categories of fast progressive (grade A), slow progressive (grade B), non-progressive (grade C) and non-motile (grade D) was based on criteria implemented in CASA systems for the analysis of frozen samples of bovine semen.However, the thresholds used for the classification were previously characterized, considering the specific conditions of our design, since for CASA systems, there is no official standard (Simonik et al., 2015).19 20

Acquisition and pre-processing
The first phase tested was the acquisition and preprocessing of the images.Figure 3C shows the results of the process in which contrast adjustment was applied in order to increase the background contrast with respect to the spermatozoa.

Segmentation and detection of spermatozoids
Segmentation by a combination of adaptive Gaussian models had effect shown in Figures 3B1 to B3.It is evident that these spermatozoids with displacement highlighted by red circles are extracted from the background together with some particles generated by the vibration of the sperm liquid, caused by the turbulence proper to the flagellar movement of spermatozoa.Some sperm tails are also extracted, either completely or partially; also, these spermatozoa have very poor or near zero displacement.Figures 3A1 to A3 highlights spermatozoids in yellow, do not exhibit movement at all and are not segmented by this technique, as expected, because the pixels that constitute them do not show variability or displacement in the video and are practically taken as part of the background.
Later, the application of morphological operators to the binary image resulting from the segmentation method allowed eliminating the small particles generated by the movement of the background.Also, the tails of the spermatozoa are reduced until they disappear, thus leaving only the heads of the spermatozoids in motion.However, the filtering process is not enough and sometimes it allows the presence of other objects of larger size in the images that do not belong to spermatic heads.For this reason, an aperture operator was applied followed by a closure operator with a circular structuring object of 2-and 4-pixel radius, respectively, in order to obtain an image with objects that belong to the heads of the spermatozoids without disturbances.
Since the adaptive Gaussian models did not segment the sperm without motility, another method was required.For this, Fisher's linear discriminant analysis was used, which allowed the sperm cell heads to be segmented independently of their mobility or vitality, in order to analyze all spermatozoa in the video.
Figure 3C shows the final frame of a test video in which this analysis was done.The corresponding result is shown in Figure 3D, where it is shown that the sperm heads are extracted from the background effectively.Fisher's method segments all sperm heads independently of their movement, since it is based only on their color composition, allowing the detection of these spermatozoa without any movement.

Assignment and construction of trajectories
In Figure 4, the effect of application of the Kalman filter and Munkres algorithm is illustrated.In Figure 4A1 to A3 the transition of three frames at different times is observed, showing the positions of all the detected centroids, labeled with red points.Whereas, in Figure 4B1 to B3, the assignment of these positions to the actual trajectories of the spermatozoids is shown with blue lines.
Likewise, the superposition of two sperm trajectories is shown in Figure 4A3 and B3.However, the Munkres assignment allows proper differentiation of these two trajectories.It was also observed that when the current position of the spermatozoon was not detected, either due to segmentation errors, the Kalman filter enabled the prediction of the position in each frame; therefore, when the target is detected again, the correction of the position should be done to compute the prediction model with the current measurement.

Characterization of sperm motility
With the Kalman filter, it was possible to define the trajectory of each moving sperm present in the video.A trajectory is described as a series of Cartesian coordinates in the image; this allows characterization and classification of each spermatozoon within a degree of motility just like the commercial CASA systems.Figure 4C shows the final frame of the trajectory as described in Figure 4A3, where it is observed that for each trajectory, the VCL, VAP and VSL are calculated to proceed to the classification stage.
Finally, in Figure 4D, the expansion of a trajectory was performed to show the classification process.In this case, a very undulating movement of a spermatozoid (WOB = 47.3%)combined with a very low linearity (LIN = 39.1%) was observed, which makes it a spermatozoon without progressivity belonging to category C.

Validation of the computational tool
The evaluation performed by the andrology expert was based on the same classification criteria described in the state of the art except for categories A and B. Categories A and B are linked in a single class (A + B) for practical effects, because, according to the expert, it is often unnecessary to distinguish between fast and low progressive spermatozoa.Taking this into account, the classification is as follows: 1. Class 1: spermatozoa with progressive motility (Category A + B); 2. Class 2: spermatozoa with non-progressive motility (Category C); 3. Class 3: sperm without motility (Category D).
From these classes, the calculation of the total sperm detected (total), percentage of progressive motility (% MP) and percentage of total motility (% MT) were derived.Thus, the observation and assessment obtained from the 10 videos by the andrology expert are shown in Table 1.It shows the number of sperm classified in each category and percentages of progressive and total motility.
On the other hand, the analysis of the 10 videos of seminal samples from different bulls was arranged and the algorithm was executed to show results in terms of the parameters of motility observed by the expert.Table 2 shows the results obtained by the computational tool.
Likewise, the determination coefficients were calculated in order to find the level of success of the computational tool developed with respect to the analysis performed by the expert.Figure 5A shows the correlation of spermatozoa in the (A + B) category for the expert and the algorithm.It is shown that the computational tool was able to classify most progressive spermatozoa with a high determination coefficient (R2 = 0.8542).
However, there were differences in the detection and classification of spermatozoa belonging to category C. Figure 5B shows that the correlation was low with a coefficient of determination of R2 = 0.3952.This was mainly due to the difficulty in adjusting a minimum sensitivity threshold in order to sort a moving and nonprogressive sperm, because the visual inspection of the expert is very subjective and there is no strict regimen, for example, when considering a spermatozoon with almost null displacement or with extremely short displacement.On the other hand, for the detection of non-motile spermatozoa, a good correlation (R2=0.8056)was observed, therefore, Fisher's discriminant analysis had a good performance, as shown in Figure 5C.
In general terms, it can be stated that the motile classification stage of the algorithm presented good agreement with the expert's manual analysis, since by combining the data of the three classes, a percentage of progressive motility of R2 = 0.8143 and high correlation with a coefficient of determination of R2 = 0.8754 can be  5D and E, respectively.Finally, using the Bland-Altman test shown in Figure 5F, which represents the difference between the methodologies implemented, no significant variability was found between the sperm motility measurement by the expert and the implemented computational tool.The outcomes show the differences of the manual measure with respect to the algorithm with a relation of 0.966 in a 95% confidence interval of (3.058 to 3.19).They indicated that the measurements were very similar, if the two methodologies were the same, the expected ratio would be 1.

DISCUSSION
In the present work, the versatility of the artificial vision is demonstrated in processes that require observation of small objects on microscopic images for long periods of time.This application provides a technological support to the manual analysis by the expert, and in this way, the intrinsic drawbacks of his intervention are greatly reduced.However, it is necessary to use an expert to properly focus the sample of the ejaculate and obtain videos of high quality.
It is important to note that the present methodology aimed only at post-thawed samples and fresh ejaculates were not supported.Though fresh spermatozoa have similar performance than the frozen counterpart (Szeptycki and Bentov, 2016), their captures under microscope and kinetics are different.
The preparation of the semen samples was carried out according to the protocols established in the World Health Organization (WHO, 2010) laboratory manual, in order to remove agglomerations of cells that could impair the procedure.Moreover, the image capture was obtained using a suitable illumination and focus configuration of a binocular microscope Leica model DM500 led, with digital camera, Wi-Fi model ICC50W.Each capture was taken immediately after thawing because sperm vitality decreases over time (Ibȃnescu et al., 2016).
Disadvantages were evidenced in some samples due to the cases in which some spermatozoa were defocused by the change of their height within the sample, which caused erroneous detections in the segmentation.The contrast adjustment allowed this situation to be improved by facilitating the method to recognize variation in pixel values along the frames.
Due to the vibration of the seminal fluid caused by the oscillatory movement of the sperm tails, there was an exaggerated detection of particles not corresponding to sperm heads; therefore, morphological operators were implemented to eliminate them in the field of the binary image.This could be avoided if a non-linear diffusion filter is used as shown by Imani et al. (2014), but morphological R 2 = 0.8154 operators were enough.In this way, only the regions belonging to the spermatic heads were obtained, which through an analysis of connected components were extracted from the values of the position coordinates in each frame.
The adaptive mixture of Gaussian models algorithm is one of the useful techniques for background subtraction as summarized by Piccardi (2004).It allowed a proper detection of movement in the cells that presented displacements, but those immobile ones were not segmented.Therefore, a Fisher´s linear discriminant analysis was implemented to extract all the sperm cells using an experimentally extracted database to perform a complete motility analysis.Fisher method offers great precision of segmentation, but consumes a high amount of memory and processor (Zhiwei et al., 2012); thus, for higher resolutions, it could take significant periods of time.
Furthermore, the criterion of eccentricity made it possible to discard non-spermatozoa objects in the segmentation process, making the method more reliable.However, in certain specific cases, segmentation errors were observed, where two sperm heads that are in contact were taken as a single object and were therefore ignored by the area criterion.Watershed algorithms and particle filters were able to break this contact in the segmentation process (Ravanfar et al., 2014).
The reconstruction of each sperm trajectory in the video was possible due to the implementation of the Munkres assignment algorithm; it could assign the new centroids to the existing trajectories using the Kalman filter prediction and correction criteria.This filter, whose initial variables were adjusted experimentally, worked correctly in the individual tracking process.A common issue found in multi-tracking algorithms is the problem derived from the superposition of targets.Since each sperm trajectory develops its own Kalman filter, this issue was solved (Urbano et al., 2017).
Despite not having an official protocol for the CASA systems, one was adapted for the present design, by characterizing the coordinates that describe the trajectory of each motile spermatozoon in the video.This allowed calculation of the parameters of individual motility according to the commercial systems CASA.Starting from the acquisition of these parameters and the classification of motility, for each video, a percentage of progressive motility was obtained, which was submitted to validation having as reference, the assessment by an expert in andrology.
The validation regarding observation and manual analysis by a technician showed high concordance in the classification of progressive motility.However, the classification of non-progressive motility was affected, mainly by the complexity of the objective determination of the speed limits or cut-off values of non-progressive and almost immotile spermatozoa (Wilson-Leedy and Ingermann, 2007), because the sensitivity of the algorithm is quite different as compared to the human eye and the conventional manual procedure presents downsides (Baracaldo et al., 2007).This could be overcome but it would require higher frames per second.
This affected the total motility analysis, but the results showed that they followed the trend line of the data provided by the expert.On the other hand, the parameter of progressive motility, one of the significant descriptors for fertility together with VAP and VSL (Ahmed et al., 2016) (Rezagholizadeh et al., 2015), had a high coefficient of determination (R2=0.8143)and, therefore, a very good correlation.
As also demonstrated by Giaretta et al. (2017), open alternatives are as competitive as commercial systems like Hamilton-Thorne IVOS and present similar performance for sperm motility analysis, even though there is no official standard for computer-assisted tools.
Thus, a semi-automatic, reconfigurable and modular computational tool was developed for the assisted analysis of individual sperm motility of bovine semen samples.It has the capacity to reduce human errors in the interpretation of motility parameter and thus provide reliable results.In Colombia, this tool would provide an economical and accessible alternative for quality analysis, academy and research in the field of andrology and veterinary medicine.
Future developments involve, firstly, the minimization of errors presented in this study, then, the design of a complete CASA system with the capacity to analyze motility and viability, concentration and sperm morphology.In addition, adjusting the design to provide the user with the option of evaluating both bovine semen samples and other animal species, is proposed.

Figure 2 .
Figure 2. Representation of the assignment of sperm detections to their trajectories.

Figure 3 .
Figure 3. (A1-A3) Sequence of original frames.(B1-B3) Sequence of frames showing segmentation of three spermatozoids in motion by the adaptive method of Gaussian models.(C) Original frame.(D) Fisher discriminant analysis segmentation.

Figure 4 .
Figure 4. Frames at different times: 1, 2, 3. (A) in red dots, unassigned centroid detections are marked; (B) in blue lines, the constructed trajectories are illustrated; (C) trajectories are characterized; VLC in green; VAP in blue; VSL in red; (D) characterized trajectory of a random spermatozoon.

Figure 5 .
Figure 5. (A) Linear regression and coefficient of determination of spermatozoa classified in A + B, (B) Linear regression and coefficient of determination of spermatozoa classified in C, (C) Linear regression and coefficient of determination of spermatozoa classified in D, (D) Linear regression and coefficient of determination of the percentage of progressive motility calculation, (E) Correlation between the computational tool vs. the expert technician, (F) Bland-Altman test for classes A+B, C and D.

Table 1 .
Results of manual analysis by the expert technician, elaborated by the author.

Table 2 .
Results obtained with the computational tool, elaborated by the author.