An assessment tool for computer-assisted semen analysis (CASA) algorithms.

Computer-Assisted Semen Analysis (CASA) enables reliable analysis of semen images, and is designed to process large number of images with high consistency, accuracy, and repeatability. Design and testing of CASA algorithms can be accelerated greatly if reliable simulations of semen images under a variety of conditions and sample quality modes are available. Using life-like simulation of semen images can quantify the performance of existing and proposed CASA algorithms, since the parameters of the simulated image are known and controllable. We present simulation models for sperm cell image and swimming modes observed in real 2D (top-down) images of sperm cells in laboratory specimen. The models simulate human sperm using four (4) types of swimming, namely linear mean, circular, hyperactive, and immotile (or dead). The simulation models are used in studying algorithms for segmentation, localization, and tracking of sperm cells. Several segmentation and localization algorithms were tested under varying levels of noise, and then compared using precision, recall, and the optimal subpattern assignment (OSPA) metric. Images of real human semen sample were used to validate the segmentation and localization observations obtained from simulations. An example is given of sperm cell tracking on simulated semen images of cells using the different tracking algorithms (nearest neighbor (NN), global nearest neighbor (GNN), probabilistic data association filter (PDAF), and joint probabilistic data association filter (JPDAF)). Tracking performance was evaluated through multi-object tracking precision (MOTP) and multi-object tracking accuracy (MOTA). Simulation models enable objective assessments of semen image processing algorithms. We demonstrate the use of a new simulation tool to assess and compare segmentation, localization, and tracking methods. The simulation software allows testing along a large spectrum of parameter values that control the appearance and behavior of simulated semen images. Users can generate scenarios of different characteristics and assess the effectiveness of different CASA algorithms in these environments. The simulation was used to assess and compare algorithms for segmentation and tracking of sperm cells in semen images.

Introduction

Computer-Assisted Semen Analysis (CASA) systems and their algorithms continue to be of great interest to clinicians and andrology researchers[1]. Modern CASA systems “have been designed to objectively and quantitatively measure several aspects of sperm structure and function, aiming to provide high levels of intra- and inter-laboratory consistency”[1,2]. To achieve this aim, methods of noise filtering, image segmentation, localization, multi-object tracking, and machine learning were employed[3-15].

A major challenge in developing and validating CASA systems is the accurate assessment and comparison of their semen analysis methods to the ground truth. In order to validate a CASA system, this process needs to take place across a representative sample of the expected semen images. For real-life samples, the ground truth is often unknown, motivating the use of high-quality image simulations with modifiable parameters for validation of CASA systems and algorithms. Simulations of this kind have the potential to help in developing automated semen analysis systems, and in comparing different candidate algorithms to each other.

To generate a simulation of semen sample video for assessment and validation, we need the following: (1) a model for the image of a sperm cell (to generate sperm cells in the simulated semen image) and (2) a model for the sperm cell movement that defines how the semen image changes over time. The development of sperm cell model in this paper mirrors similar work that was done for other cell types[16-21]. The development of sperm cell swim models is informed by several existing approaches. Some of these use a set of nonlinear equations of motion (Armon et al.[22] and Ch. 6 of Urbano (2015)[4]). Others use fluid dynamics to develop a more refined description of cell movements[23-31]. However, most existing simulation models of cell dynamics[25-31] did not integrate multiple cells into a comprehensive image such as the ones used in practice by CASA systems[32,33]. Integration of these single cell descriptions into full-scale multi-cell images is required before relevant image simulations can be developed and used.

In this paper, we provide models for the image of a sperm cell and for four (4) different swimming modes of sperm cells (circular, linear mean, hyperactive, and immotile). The computational requirements to simulate these models and integrate them into multiple-cell image are relatively low, allowing for generation of multiple images in a format similar to the images used by CASA systems. As an example, a snapshot of a real semen image and simulated image are shown side by side in Fig. 1. Figure 1a is an image of semen sample at magnification. Figure 1b is a simulated semen image generated by using parameters estimated from the semen image on the left.

Figure 1

(a) Image of real human semen sample. (b) Image of simulated semen sample.

The rest of this paper is organized as follows. In the “Methods” section, we present a model of a sperm cell and of the swimming modes. In the “Testing setup” section, we explain how simulated images were created to mimic real ones for visual and numerical comparisons. In the “Testing” section, we demonstrate use-cases for segmentation, localization, and tracking algorithms.

This paper is accompanied by simulation software which is publicly available at Choi et al.[34] as a stand-alone software and as MATLAB codes (https://github.com/JiwonChoi-NJIT/NJIT_sperm_simulator). The dataset supporting the conclusions of this article are available in the github repository (https://github.com/moshekam/NJIT-Semen-Images-Data-Fusion-Lab)[35].

Methods

The human sperm cell (see Fig. 2, retrieved from Villarreal (2006)[36]) has two main parts, the head and the flagellum. The flagellum is composed of midpiece, tail (also called principal piece), and end piece. The flagellum moves in a wave-like motion to propel the cell[37]. In the following sections we describe the 2-dimensional (2-D) image of the head and the flagellum of a sperm cell, and the models of sperm movements.

Figure 2

Diagram of a human sperm cell[36]. (Top) Main components of a sperm cell. (Bottom) Front and side view of sperm head and midpiece.

Simulating a sperm cell

We generate the image of the head and the flagellum of a sperm cell separately, and then combine them. In modeling the head and the flagellum of a sperm cell, we follow the morphology of sperm cells described in the WHO laboratory manual for the examination and processing of human semen[38]. The normal shape of the human sperm head is generally oval. The tail, or the principal piece, of the flagellum is a thin cylinder of uniform calibre. In our simulation, the flagellum is modeled by the tail only (we do not consider the midpiece or the end piece for the simulation model since they are often too small to be observed separately in the resolution used by CASA systems).

Figure 5

Flowchart of sperm head (center) image generation process. Image and point spread function are the inputs. The output is the simulated image of sperm head center . and is convolved and then scaled by to generate (Processes A-1 and A-2). The image is complemented to generate the image (Process A-3). Lastly, the background is added to image , generating an image of sperm head center (Process A-4).

Figure 7

Flowchart of the image generation process of the sperm membrane. Image and point spread function are the inputs. The output is the simulated image of sperm membrane . and is convolved and then scaled by to generate (Processes B-1 and B-2).

Figure 8

Diagram of sperm head image merging process. Images and are the inputs. The output is the simulated image of sperm head . and are added to generate (Process C-1).

Figure 9

Flowchart of flagellum generation process. Image and point spread function are the inputs. The output is the simulated image of flagellum . and are convolved and then scaled by to generate (Processes D-1 and D-2).

Figure 10

Diagram of sperm head and flagellum image merging process. Images and are the inputs. The output is the simulated image of sperm cell . and are added to generate (Process E-1).

The flowchart of the process of generating the sperm image is shown in Fig. 3. The inputs of the process are images (corresponding to the sperm head) and (corresponding to the curve of the flagellum), and point spread functions , , and . The image of sperm head center and membrane are generated using the operations labeled Process A and B, respectively. The inputs for Process A are image and point spread function , and the inputs for Process B are image and point spread function . The images of sperm head center and membrane are merged in Process C, producing a final image of sperm head . Processes A-C are therefore the sperm head generation process. The final image of the flagellum is generated using Process D, whose inputs are image and point spread function . Lastly, the image of sperm head and flagellum are merged in Process E, producing the complete image of the sperm cell, . A grayscale image of the simulated sperm cell is shown just below the full-color simulated cell image in Fig. 3.

Figure 3

Flowchart of sperm image generation process. Image and , and point spread functions , , and are the inputs. The output is the simulated sperm image. Details of each process are shown in Figs. 5, 7, 8, 9, and 10. Image shows the image of sperm head (Output of process C). Image shows the image of flagellum (Output of process D). Image shows the image of sperm head (Output of process E). Color bars next to - indicate the corresponding grayscale intensity values. Resulting grayscale image of simulated sperm cell below the colored plot of .

is an image of N by N pixels whereN is the size of the frame of the simulation image. Here, , , , are the locations of sperm heads of cells engaged in circular swimming, linear mean swimming, hyperactive swimming, or no swimming (immotile cells), respectively. is also an image of N by N pixels whereHere, , , , are the points along the curve of sperm flagellum of circular swimming, linear mean swimming, hyperactive, and immotile cells, respectively. Each flagellum consists of M points (typically ). The detailed placing of the sperm head (, , , ) and the points along the curve of sperm flagellum (, , , ) are provided below in the section titled “Swimming models”.

An image of a simulated sperm cell is shown in Fig. 4 with the points for the location of sperm head (where ) and the curve of the flagellum (where ) shown in green and blue, respectively. In the simulation, grayscale values are used in the range between 0 to 255 (256 levels). In following sections, the five processes (A-E) are explained in detail. The flowcharts of processes (A-E) are given in Figs. 5, 7, 8, 9 and 10. Images of - are accompanied by numbered colorbar that can be used to convert the image to grayscale.

Figure 4

Simulated sperm image with points for the location of the head (green) and the curve of the flagellum (blue).

Process A: sperm head (center) generation (Part 1)

The flowchart of Process A is shown in Fig. 5. To generate the oval-like structure of the sperm head, image is convolved with a 2-D normal distribution filter (Process A-1); the filter is used as a point spread function (Ch. 3–5 in Gonzalez[39]). The point spread function , a 2-D normal distribution, is defined asThe 3-D surface plot of is shown in Fig. 6a. Standard deviations and control the length and width of the cell head. In the examples shown in this section, px and px. The size of the filter is px (x and y ranges from − 12 to 12). The resulting output image is then scaled by a constant to set the peak intensity value to be 255 and is denoted (Process A-2 in Fig. 5). Constant can be changed to control the intensity value of cell head. The image is complemented to produce the image of dark-oval-like center on a white background (Process A-3). In Process A-4, the background for the simulation image is added. Let be the background that will be added to the simulated image. The size of is N by N. In the examples in this section, a uniform background is assumed, where every value in the pixel is equal to 204 (, 80% of 255). The value is subtracted from the image and any value below 0 is set to 0. The resulting image is denoted image ().

Figure 6

Point spread functions used in sperm cell image generation process; (a) , (b) , (c) .

Process B: sperm head (membrane) generation (Part 2)

The flowchart of Process B is shown in Fig. 7. In typical semen image samples, one can observe a halo-like membrane surrounding the sperm head, as shown on the left of Fig. 1 (referred to as “a horseshoe-shaped halo” in Urbano et al.[3]). To generate this halo-like membrane around the cell, a modified version of Laplacian of the 2-D normal distribution filter is used as the second point spread function on image . The Laplacian of the 2-D normal distribution filter is defined asThe point spread function is defined asThe resulting point spread function is shown in Fig. 6b. Standard deviations and control the length and width of the cell membrane. In the examples in this section, px and px. The size of the filter is px (x and y ranges from − 12 to 12). The point spread function is convolved with image , generating the image of a membrane (Process B-1). The resulting output image is then scaled by constant (Process B-2) to set the peak intensity to be 51 (peak intensity value of cell membrane, 20% of 255). This image of the membrane is denoted . Constant can be changed to control the intensity value of the cell membrane.

Process C: sperm head generation (Part 3)

The flowchart of Process C is shown in Fig. 8. In Process C, the images of the sperm head center and the image of the flagellum are added to complete the image of sperm head. The resulting sperm head image is shown in Fig. 8 and is denoted .

Process D: sperm flagellum generation

The flowchart of Process D is shown in Fig. 9. To generate the flagellum, image (Eq. 2) is convolved with point spread function (Eq. 6) to generate the image of uniform calibre look of the flagellum. The point spread function is defined as the Laplacian of Gaussian filter (equal to the Eq. (4) with ), namelyThe 3-D surface plot of is shown in Fig. 6c. Standard deviation controls the width of the flagellum. In the examples in this section, px. The size of the filter is px (x and y ranges from -12 to 12). The output image is then scaled by constant to set the peak intensity of the membrane as 13 (approximately 5% of 255). Constant can be changed to control the intensity of the flagellum. The resulting image of the flagellum is a long membrane-like curve. The scaled image of the flagellum is denoted .

Process E: sperm head and flagellum merge

Finally, in Process E, the image of cell head, , and the image of the flagellum, , are added to generate the image of a sperm cell ( in Figs. 3 and 10). Process E is shown in Fig. 10.

Swimming models

In this section, we describe the swim/movement models of the head and the flagellum of a sperm cell. The swimming models describe how the position of the head and the flagellum change over time. These positions define the points used to generate (cell head, Eq. (1)) and (cell flagellum, Eq. 2).

We categorize the sperm movement into the following four swimming modes: Figure 11 shows simulated tracks for the four swimming modes (4 seconds long). The tracks indicate the location of the head of a cell and the arrow indicates the direction of movement.

Figure 11

Simulated swimming path of the four swimming modes: circular swim, linear mean swim, hyperactive, and immotile. The color bar indicates the color of a track with respect to time. The duration of movement is 4 s.

Circular swim,

Linear mean swim,

Hyperactivated,

Immotile, or dead.

Circular swimming cells and linear mean swimming cells actively propel themselves forward, either on a large circular or linear path. Examples of (1) circular swim and (2) linear mean swim are shown on top of Fig. 11. The red arrow indicates the direction of movement for the cell. (3) Hyperactivated sperm cells do not travel significantly away from their initial location. The movements of hyperactive cells are described as “vigorous”, “whiplash type”, or “frantic”[40,41]. An example is shown in bottom left of Fig. 11. (4) Immotile movements define cells which show barely any movement at all (shown in bottom right of Fig. 11).

We synthesized the models to describe circular and linear mean movements based on the sperm tracks recorded in the literature[3,8,29,33,42]. Hyperactive swim, described to be “frantic” and “vigorous” is modeled through as Brownian motion to reflect the randomness in the cell movement[43]. In following sections, we describe each movement in detail.

Circular swim model

A sperm cell exhibiting a circular swim moves along a circular path with oscillations about this path caused by the beatings of the flagellum[33]. This circular swim can be represented as a sinusoidal modulated circular path (Fig. 12). The equations used to model the movement of the head and the flagellum of a circular swimming cell are given in Table 1. The 2-D position of the cell’s head is determined by Eqs. (7a), (7b). is the radius of the overall circular path, is the frequency of the sinusoid modulated on the circular path (Hz), a is the amplitude of the sinusoid modulated on the circular path, and is the frequency of the circular cycle (cycle/sec). and are horizontal and vertical offsets, respectively.

Figure 12

Trajectory of simulated circular swimming cell.

Table 1

Simulation of a circular swimming cell.

Circular swim model—head	x_H_C(t)=(r_c+asin(2πf_st))cos(2πf_ct)+C_x_c	(7a)
	y_H_C(t)=(r_c+asin(2πf_st))sin(2πf_ct)+C_y_c	(7a)
Circular swim model—flagellum	x_T_c(k,t)=b(k)sin2πkM-f_st	(8a)
	y_T_c(k)=-λkM	(8b)
	x_tail_Cy_tail_C=R(2πf_st)x_T_cy_T_c+x_H_Cy_H_C	(9)
	R(·)=cos(·)-sin(·)sin(·)cos(·)	(10)
	b(k)=aαλkM+β,α=0.02,β=0.8	(11)

head position of circular swimming cell.

radius of the circular path.

a :  amplitude of the sinusoid modulated on the circular path.

frequency of the sinusoid modulated on the circular path (Hz).

frequency of the circular cycle (cycle/s).

vertical and horizontal offset constant.

a set of k points along the center of the flagellum of circular swimming cell.

wavelength of flagellum (distance between the start and the end of the flagellum).

rotation matrix.

local variation in beating amplitude along the flagellum.

are Cartesian coordinates of the points along the center curve of the flagellum with , where in the default setting we used . The flagellum of circular swimming cells follows Dresdner’s model of sperm flagellum (Eq. 8a)[24]. Equation (8a) determines the oscillation of the flagellum (), and Eq. (8b) determines the end to end length (wavelength ) of the flagellum (distance between the start and the end of the flagellum). b(k) is the local variation in beating amplitude along the flagellum. The points are rotated using the rotation matrix to match the direction of movement, and shifted to the location of the sperm head () (Eq. 9); the points of the flagellum of a circular swimming cell are denoted . In the examples shown in this paper, the local variation in the beating amplitude, b(k), follows the affine function with for circular swimming cell. Here, a is the amplitude of a sinusoid modulated on the circular path, and is the wavelength of the flagellum.

An example of simulated flagellum of a circular swimming cell is shown in Fig. 13. On the left is a simulated image of a circular swimming cell with the path shown in blue line. On the right is a curve which represents the flagellum of that cell. The example in this image uses amplitude of pixels and wavelength of pixels.

Figure 13

(Left) Simulated image of circular swimming cell with past track shown in blue. (Right) Plot of the flagellum of circular swimming cell. In this example, wavelength () is 40 pixels and amplitude (a) is 4 pixels.

Linear mean swim model

In linear mean swim, a sperm cell moves along a straight-line path, rolling side to side as it propels itself forward. This feature causes a ribbon-like movement about the straight line. The equations used to model the movement of the head and the flagellum of linear mean swimming cell are given in Table 2.

Table 2

Simulation of a linear mean swimming cell.

Linear mean swim model—head	x_cy_c=Vcos(θ_r)t+C_x_LVsin(θ_r)t+C_y_L	(12)
	P=P_x(t)P_y(t)=r_v2sin(4πf_lt)r_hA_c2[sin(2πf_lt)+A_harsin(6πf_lt)]	(13)
	A_c=1maxθ(sin(θ)+A_harsin(3θ))	(14)
	x_H_Ly_H_L=R(θ_r)P+x_cy_c	(15)
Linear mean swim model—flagellum	x_o(k,t)y_o(k,t)=b_1(k)x_T(k,t)b_2(k)y_T(k,t)	(16)
	x_T(k,t)=r_v2sin4πkM+f_lt	(17a)
	y_T(k,t)=r_h2sin2πkM+f_lt	(17b)
	x_LM(k,t)y_LM(k,t)=x_o(k,t)y_o(k,t)-λkM0	(18)
	y_LM_2(k,t)=y_LM(k,t)-[P_y(t)-y_T(0,t)]b_3(ψ)	(19)
	x_tail_LMy_tail_LM=R(θ_r)x_LMy_LM_2+x_H_Ly_H_L	(20)
	b_1(k)=11+e^(αk/M+β)	(21a)
	b_2(k)=e^-γ_1k/M	(21b)
	b_3(k)=1-e^-γ_2k/M	(21c)

position of the sperm head of linear swimming cell .

horizontal and vertical offset constant (pixels).

V :  straight line path velocity (pixels/sec).

rate of change in ribbon angle (Hz).

width and height of ribbon (pixels).

user defined ratio between the first and the third harmonics.

correction constant for defined width of the ribbon.

direction of the forward movement (radian).

a set of k points along the center of the flagellum of linear mean swimming cell.

local horizontal and vertical variation in beating amplitude along the flagellum.

flagellum position correction function.

We generate the cell head movement (ribbon-like movement along a linear path) using Eqs. (12–15). The position along the linear path is denoted () and the position along the ribbon-like movement is denoted (). The resulting path is shown in Fig.  14a.

Figure 14

(a) Trajectory of simulated linear mean swimming cell (Eqs. 12–15). (b) Ribbon-like oscillatory movement along the straight-line path of linear mean swimming cells (Eq. 13).

The straight-line path is a linear function of time t with initial positions of . The straight-line path is a line segment from initial positions to the point () (as shown in Fig. 14). The horizontal and vertical positions along the linear path () change at the rate of horizontal and vertical velocity in the direction of forward movement (Eq. 12). The position along a “ribbon-like” path () is generated using Eqs. (13–14). The ribbon is shown in Fig. 14b, where is the height of the ribbon and is its width. and are periodic functions of time t. is a sinusoid with a frequency of and amplitude of . is a sum of two sinusoids of frequencies of and and amplitudes of and , respectively. is a correction constant to set the amplitude of function as (Eq. 14). The resulting ribbon-like path is shown in Fig. 14b. is the ratio between the two sinusoids of function . The examples in this paper use .

The position of the cell along the ribbon, (), is multiplied by rotation matrix resulting in a rotation by (Eq. (15), ). Lastly, The position on the ribbon-like path is added to the position along the straight-line path to define the position of the sperm head at all times (Eq. 15).

The flagellum is generated using Eqs. (16–21). Dresdner’s model is used to define the horizontal and vertical oscillations in the flagellum (). and are discrete set of points for realized at time t. The value of M in our simulation is 200. An example of simulated flagellum of linear mean swimming cell is shown in Fig. 15. The blue dots represent the successive location of the sperm head. We also show, for one of these locations, the associated flagellum in orange. We use the fundamental frequency of the ribbon in horizontal ( for ) and vertical directions ( for ) to determine the movement of the flagellum (Eqs. (13) and (17)). The local variation in beating amplitude for and are and , respectively (Eqs. (21a) and (21b)).

Figure 15

Example of sperm flagellum generation for linear mean swimming cell. (Orange) Simulated flagellum. (Blue) Track of linear mean swimming cell.

The function (Eq. (21a)) is a transformed sigmoid function, where the values of and determine the shift and the compression/expansion of the function; (Eq. 21b) is a decaying exponential function. Together, these functions provide the desirable realistic visual effect.

We define a line segment from the head to the end of the flagellum (). The oscillations and are added to this straight line (Eq. 18).

The differences in the vertical location of the head and the vertical position of the flagellum cause offset in placement of the flagellum because we have only considered the fundamental frequency component, , of to generate (). This mismatch is corrected by Eqs. (19) and (21c).

Lastly, the flagellum is rotated to match the direction of movement and shifted to the position of the cell head (Eq. 20). The constants in Eqs. (21a), (21b), and (21c) in the example of we provide here were set to: .

Hyperactive swim model

The equations used to model the movement of the head and the flagellum of hyperactive swimming cell are given in Table 3. For each dimension (x and y), the movement of a cell follows Eq. (22), where ( are drift coefficients of Brownian motion, () are diffusion coefficients of Brownian motion, and W(t) is the standard 1-dimensional Brownian motion (). Our simulation assumed no drift (), and the diffusion coefficients are assumed to be the same for the horizontal (x) and the vertical (y) directions (). The location of the sperm head of hyperactive cell is defined by Eq. (23), where T is the difference in time between each simulation frame ().

Table 3

Simulation of a hyperactive swimming cell.

Hyperactive swim model—head	x(t)=μ_xt+σ_x_bW(t)and	(22a)
	y(t)=μ_yt+σ_y_bW(t)	(22b)
	x_H_H(t)=x_H_H(t-T)+σ_bW(T)and	(23a)
	y_H_H(t)=y_H_H(t-T)+σ_bW(T),	(23b)
Hyperactive swim model - flagellum	(x_tail_H(k,t),y_tail_H(k,t))=       Linear interpolation on a set of points       [(x_H_L(0),y_H_L(0)),(x_H_L(1),y_H_L(1)),       (x_H_L(2),y_H_L(2)),…,(x_H_L(n),y_H_L(n))] wherex_H_L(n)=x_H_H(t-nT),       y_H_L(n)=y_H_H(t-nT)	(24)

drift coefficients of Brownian motion.

diffusion coefficients of Brownian motion.

W(t) :  standard 1-dimensional Brownian motion ().

position of the sperm head of hyperactive cell.

T :  the difference in time between each simulation frame.

a curve along the center of the flagellum of hyperactive cell.

The flagellum of the hyperactive swimming cell is defined to be the linear interpolation on a set of points (a set of location of the cell from time to t, Eq. (24)). In our simulations, the value n is set to or rounded up to the nearest integer value.

Immotile or dead cell model

Immotile or dead cells are simulated as non-moving objects. The equations used to model the movement of the head and the flagellum of immotile cell are given in Table 4.

Table 4

Simulation of an immotile cell.

Immotile cell model—head	x_H_I(t)=x_H_I(0)and	(25a)
	y_H_I(t)=y_H_I(0)	(25b)
Immotile cell model - flagellum	(x_tail_I(k,t),y_tail_I(k,t))=       Linear interpolation on a set of points       [(x_I_L(0),y_I_L(0)),(x_I_L(1),y_I_L(1)),       (x_I_L(2),y_I_L(2)),…,(x_I_L(n),y_I_L(n))] wherex_I_L(n)=x_H_H(-nT), y_I_L(n)=y_H_H(-nT)	(26)

position of the sperm head of immotile cell.

a curve along the center of the flagellum of immotile cell.

Throughout the simulation, the head and the flagellum of the cell will not move (Eqs. 25, 26). The flagellum is generated by taking a snapshot of the flagellum of a hyperactive cell. The flagellum of the immotile cell is generated using Eq. (26).

In Fig. 16, an example the simulated semen sample is shown, where the cells of four different swimming modes can be seen. The trajectories of each cell are shown in blue line. In the simulation software (available in Choi et al.[34]), parameters for cell concentration, appearance, and swimming modes can be changed to generate the user’s desired semen image for testing.

Figure 16

Example of simulated image with track of each cell shown in blue.

In Fig. 17, we show 3 additional scenarios which are simulation features. These simulation features are used to generate dynamic scenarios for testing of CASA algorithms and systems.

Figure 17

(a) Simulation image of two cells with additive random noise; tracks of the two cells shown in blue. The random noise in this example is Gaussian random variable. (b) Simulation image of sperm cells with variable intensity. (c) Simulation image of a sperm cell making transition to other swimming modes.

The first feature adds noise to the position of each sperm cell (example shown in Fig. 17a). The addition of noise is used to mimic the random movements caused by the cell and the surrounding fluid.

The second feature assigns a different intensity to each sperm cell (example shown in Fig. 17b). This feature is to simulate an environment where the sperm cells appear differently due to the depth of the chamber or slide used to observe the semen sample. The chamber that is used to observe semen is usually 10 to m deep[38]. Certain cells may be swimming much deeper within the chamber when other cells may be swimming close to the top of the chamber. This causes some cells to look faded. The second feature is used to mimic such scenario.

Sperm cells change their swimming mode as time progresses[32]. The third feature assigns transition probabilities to define how likely a cell will change from one swimming mode to the other. In the example shown in Fig. 17c, a sperm cell goes from linear mean swim to circular swim to hyperactive swim. As an example, one can assign probability that a cell undergoing linear mean swim has 85% chance of staying in linear mean swim, 10% chance of transitioning to circular swim, and 5% chance of transitioning to hyperactive swim after 1 second (0% chance to change into immotile). In Choi et al.[34], we provide sample images generated using these three different functionalities as examples.

Testing setup

For testing, we simulated images based on real semen image samples. The semen samples used in the study were collected by us for previous studies[3,4] and are available in the public domain[35]. The images of the samples were provided by the In-Vitro Fertilization laboratories at Penn Fertility Care. Each specimen was allowed to liquefy for 30–40 minutes at room temperature and was washed in media. The washed semen samples were pipetted on a 20 m deep Vitrolife MicroCell chamber for data collection. For detailed explanation on the preparation process of the samples, please refer to Urbano et al.[3] and Urbano[4].

The parameters we have extracted from the real images were the image background , the noise variance , the size of sperm head, the number of cells in the image , and the number of non-moving cells . The values of the parameters are given in Table 5. The sperm cells in the image were labeled manually. The process used to calculate these parameters are provided in the additional file [see Additional file 1].

Table 5

Simulation parameters for segmentation and localization testing.

	Noise variance σN2	Radius r_M (major axis)	Radius r_m (minor axis)	Number of cells N_C	Number of non-moving cells N_D
Sample 1	8.22×10^-6	2.86 px	1.86 px	10	3

Testing

Two use cases of the simulation were explored. First, we tested five sperm cell detection algorithms and assessed their performances. The cell detection algorithms consists of two major components, which are segmentation and localization[39]. In segmentation, regions in the image that contains sperm cells are separated from the image. In localization, the locations of sperm cells in the segmented regions are labeled. The performance of segmentation and localization was evaluated in terms of optimal subpattern assignment (OSPA) distance ()[44], precision and recall rates[45]. OSPA distance is a metric that quantifies errors in distance between the ground truth and detection and in cardinality (difference in number between ground truths and detections). Precision is the ratio between true positives (true matches) to the total number of detections. Recall is the ratio between true positives to total number of ground truth. Ideally, the value of OSPA distance is 0 and the values of precision and recall are 1. The calculation of OSPA distance, precision, and recall in real samples was made possible by manual labeling of sperm locations.

Second, we tested four tracking algorithms for sperm tracking (nearest neighbor (NN), global nearest neighbor (GNN), probabilistic data association filter (PDAF), and joint probabilistic data association filter (JPDAF))[46]. The purpose of a tracking algorithm is to track each of the sperm cells as time progresses. The calculated tracks are then used to provide parameters that describe the motility of sperm cells in the semen sample (e.g., percent of motile cells, velocity of sperm cells). The performance of the tracking algorithms was assessed in terms of Multiple Object Tracking Precision (MOTP) and Multiple Object Tracking Accuracy (MOTA) with cutoff distance of px[47]. MOTP is the sum of the distance between the matched ground truth and detected tracks. MOTA is a metric that quantifies the error in false positive (false alarm), false negative (missed detection), and track mismatch . The MOTA is defined to be . Ideally the value of MOTP is zero and MOTA is 1 (). Detailed explanation of the assessment metrics is given in the additional file [see Additional file 1].

Applying the simulation to assess segmentation, localization, and tracking

The real semen image and the simulated images were tested using five (5) different algorithms for detection (segmentation and localization). These algorithms were the following:

(Otsu) Binarization of image using Otsu’s thresholding followed by morphological enhancements (closing, dilation, and erosion)[48].

(Adaptive) Binarization of image using the adaptive thresholding method of Bradley[49] with the sensitivity defined to be 0.8, followed by morphological enhancements (closing, dilation, and erosion).

(Spot-enhancement) Binarization of spot-enhanced image using Otsu’s thresholding, followed by morphological enhancements (closing, dilation, and erosion) (method proposed in Urbano et al.[3]).

(Edge-detection) Edge detection of median-filtered image using Sobel operator (modified algorithm proposed in Abbiramy and Shanthi (2010)[5]) followed by morphological enhancements (dilation, closing, and erosion).

(GMM) Motion detection algorithm using Mixture of Gaussian Model (GMM) with the number of training frames set to 20, number of Gaussian modes in the mixture model set to 3, learning rate set to 0.005, and background ratio set to 0.7[11,50] followed by morphological enhancements (closing, dilation, and erosion).

A total of 20 different simulated images were generated for testing using the parameters obtained from a real human sample (sample 1). The size of the frame was px. In the simulated images, the non-moving cells were modeled as immotile cells. The moving cells () were modeled as either linear mean swimming () or circular swimming () cells. The algorithms were tested on real and simulated images. To each one of the images (real and simulated), we added zero-mean Gaussian noise with variance ranging from 0 to 1225 (standard deviation of 35 in grayscale). Figure 18 shows the OSPA distances, precision, and recall rates for varying levels of noise of real and simulated images.

Figure 18

OSPA distance, precision and recall rates of sample 1 for varying levels of additive Gaussian noise (a,c,e) real, (b,d,f) simulation.

The graphs show that the Adaptive algorithm exhibits the best performance in terms of OSPA (smallest distance away from ground truth) under low noise levels (noise variance approximately below 600), followed by the Spot-enhancement algorithm, the Edge-detection algorithm, the GMM algorithm, and the Otsu algorithm, respectively (Fig. 18a,b).

In terms of precision, sharp loss in performance is observed for the Otsu algorithm and the Spot-enhancement algorithm when the noise level increases beyond the noise variance of 100 and 300, respectively. Precision for the Adaptive algorithm and the Edge-detection algorithm also decreases as the noise level increases, where the Edge-detection algorithm becomes more robust than Adaptive algorithm in high levels of noise (approximately above noise variance of 600). The GMM algorithm shows the highest level of robustness, having almost no degradation in performance in precision (small number of false alarms).

In terms of recall, almost no change is observed between the noise variance of 0–1225 for the Adaptive algorithm, the Spot-enhancement algorithm, and the Edge-detection algorithm. The Otsu algorithm and the GMM algorithm showed degradation in performance when noise level increased, the Otsu algorithm showing sharp decline around noise variance of 100 and the GMM algorithm showing gradual decline between noise variance value of 0–1225.

The OSPA distance, precision and recall rates of the five different algorithms show similar trends for the real sample image and its simulated image. Overall, the Adaptive algorithm performed best in low level of noise (below noise variance of 600) and the Edge-detection algorithm has performed the best in high level of noise (above noise variance of 600). The Otsu algorithm showed the worst performance. Spot-enhancement algorithm performed well in low levels of noise, but failed when the noise level rose above noise variance of 300. GMM algorithm was the best algorithm in terms of precision, however did not show good performance in terms of OSPA distance and recall. Additional results for another real human sample (sample 2) are provided in the additional file [see Additional file 1].

Applying the simulation to assess tracking

Four (4) tracking algorithms (NN, GNN, PDAF, and JPDAF) were tested on simulated images. The codes for tracking algorithms were written by Leonardo Urbano and are available in Urbano et al.[51]. The output of each tracking algorithm is the track information of all the detected cells in the semen image.

We evaluated the performance of NN, GNN, PDAF, and JPDAF algorithms for varying numbers of cells. Each simulated image consisted of 20, 40, 100, or 200 cells. Each cell in each image was assigned to a swim type using equal probabilities: linear mean swim (), circular swim (), hyperactive swim (), and immotile (). The size of the frame was px. The framerate was set at 15 FPS and total of 10 seconds of each image were used for tracking. Twenty (20) different scenarios were generated for each number of cells for a total of 100 images. The background intensity level was 204 (80% of 255) and no noise was added to the video images. The simulation parameters used to generate the images are given in Table 6. Ground truth tracks for each image were provided by the simulation, and were compared to the estimated tracks calculated by NN, GNN, PDAF, and JPDAF algorithms.

Table 6

Simulation parameters used for tracking assessment.

Parameter	Value
r_c	80 pixel
f_c	50 deg/s
f_s	4 Hz
a	3 pixel
f_l	3 Hz
r_h	12 pixel
r_v	8 pixel
A_har	0.1
V	50 pixel/s
σ_x_b,σ_y_b	10 pixel

The mean MOTP, , , , and MOTA values for 20 different scenarios for each type of images (20, 40, 100, and 200 cells) are shown in Table 7. The values of multi-object tracking precision (MOTP) for the four tracking algorithms were approximately equal regardless of the number of cells in the image. As the number of cells increased, the value of multi-object tracking accuracy (MOTA) decreased and the false positive rate (), miss detection rate (), and mismatch rate () increased.

Table 7

Tracking performance of NN, GNN, PDAF, and JPDAF algorithms on varying number of cells.

# of cells	MOTP	FP¯	M¯	MME¯	MOTA
NN
20	1.4 px	0.2101	0.0941	0.0003	0.6955
40	1.4 px	0.2935	0.1143	0.0012	0.5910
100	1.4 px	0.3593	0.1609	0.0030	0.4768
200	1.4 px	0.4140	0.2163	0.0061	0.3635
GNN
20	1.4 px	0.0690	0.0812	0.0012	0.8487
40	1.4 px	0.0920	0.0886	0.0028	0.8165
100	1.4 px	0.1317	0.1090	0.0070	0.7524
200	1.5 px	0.1810	0.1343	0.0144	0.6703
PDAF
20	1.4 px	0.2076	0.0994	0.0003	0.6927
40	1.4 px	0.2934	0.1292	0.0009	0.5765
100	1.4 px	0.3636	0.1956	0.0027	0.4381
200	1.5 px	0.4448	0.2748	0.0049	0.2754
JPDAF
20	1.4 px	0.0714	0.0817	0.0012	0.8457
40	1.4 px	0.0969	0.0907	0.0026	0.8098
100	1.4 px	0.1291	0.1099	0.0058	0.7553
200	1.5 px	0.1717	0.1352	0.0118	0.6813

For large number of cells, GNN and JPDAF algorithms performed better than NN and PDAF tracking algorithms in terms of MOTA. The major factor for this difference was in the false positive rate (). Comparing NN and PDAF, MOTA for PDAF was lower than NN for all of the testing cases (20, 40, 100, and 200 cells). The GNN and JDPAF tracking algorithms showed similar performance in terms of MOTA. Overall, GNN and JPDAF showed best performance, which was followed by NN. The PDAF tracking algorithm had the worst performance.

Conclusion

We presented a model of 2-D (top-down) view of a sperm cell and of four (4) different swimming modes generated by observing the swimming paths of real human sperm cells. The simulation opens up opportunities for methodical study and comparison of different semen image processing algorithms, including algorithms for segmentation, localization, and tracking. In our examples, we tested five different segmentation and localization algorithms and used the simulation to rank the algorithms by their performance, obtaining ranking with simulation that appear similar to ranking using real images. In addition, we compared the performance of four different tracking algorithms (NN, GNN, PDAF, and JPDAF) on the simulated images using MOT metrics and have ranked them by their performances.

The simulation models and the software presented in this paper serve as a powerful new tool for developing and enhancing CASA systems and algorithms. Using this new tool, stronger and more robust CASA systems can be developed. The use of such systems is an attractive alternative to manual semen collection and assessment. Specifically, clinicians can generate and demonstrate to students and technicians a variety of images that represent different types of observed semen samples; this task can be accomplished without having to collect, process, and record human semen images. The simulated images can also be used to train technicians, by comparing the results of their analyses of the images to the simulations’ ground truth.

Supplementary Information.