of performance for this class of "pipeline" DCs counters. Section 2.4 summarizes the level
of effort required to build an ADS that relies on conventional computer vision techniques.

DL-based object detection models like Faster R-CNN [9] might require less design work and perform better, but large labelled datasets are required. Therefore, the next section investigates chromosome generation to synthesize labelled datasets.

**3.3 Simulation of chromosome patches**

Investigating chromosome generation is appealing for two main reasons. First, for VAEs, investigating links between latent variables and morphological variations might help one to build a better classifier. Instead of relying on a logistic regression, simple functions of some latent variables might be enough for classification. Second, a good chromosome simulator could be used to generate fully annotated metaphase images, which could be used to pre-train aberration detection models, or evaluate architectures on artificial data.

**3.3.1 Simulating chromosomes with VAEs**

Simulating chromosomes with a VAE can be done by sampling the latent space prior (which is a Gaussian with zero mean and unit variance), and decoding those samples into image space. While this is what the authors of [36] suggest, there are several caveats with this approach. First, as pointed in the introduction of [44], the true posterior over the latent variables rarely matches the prior. Therefore, sampling from this unit variance Gaussian prior might lead to low quality samples. To circumvent this issue, we investigate two alternative methods to sample from this posterior.

First, we projected latent space samples in an orthogonal basis, using PCA. We estimated the probability distribution of PCA loadings for every coordinate of the latent space using a Kernel Density Estimation (KDE). We estimate the bandwidth using the cross-validation algorithm implemented in [100]. A grid of chromosome samples is displayed in Figure 3.7 a). Those samples show low diversity, and are not very realistic. As the dimension of the latent space is quite high (256), density estimation methods struggle because of the high distance between samples. Therefore, we fitted a Gaussian Mixture Model (GMM) on the latent space of the VAE. Again, we sampled from this GMM, and used those samples as an input for the decoder of a pre-trained VAE. A grid of samples is shown in Figure 3.7. This improves samples quality and diversity, but introduces arti-
<figure><img src="image_53.png" /><figcaption>Figure 3.7 – a): Samples of the latent distribution estimated with KDE. b): Samples of the latent distribution estimated with a Gaussian Mixture Model. In both cases, those samples were decoded with the decoder of the VAE trained earlier.</figcaption></figure>
<figure><img src="image_54.png" /><figcaption>Figure 3.8 – 8 DC samples from pix2pix. Top 4 rows are "skeleton" inputs, bottom 4 rows are chromosome outputs.</figcaption></figure>
facts. Overall, we found generating samples in a completely unsupervised fashion to be
very difficult.

3.3.2 Simulating chromosomes with pix2pix

We use the pix2pix architecture [43] to synthesize chromosomes as illustrated by Fig-
ure 3.8 which shows a grid of samples from pix2pix. While chromosome samples have no
artifacts, distinguishing MC and DC chromosomes remains difficult. Thin chromosomes
often have completely indistinguishable chromatides which makes telling wether the syn-
thetic chromosome is a MC or a DC basically impossible. This prevents one from using
this technique to synthesize artificial metaphases to pre-train object detection models.

3.4 Training Faster R-CNN on a synthetic dataset

We trained Faster R-CNN (see Section 1.4.1 for a description of this model) on the synthetic datasets introduced in Section 3.1.2. To evaluate the relationship between model performance and the number of training image, we trained several models on 100, 200, 300, 400, 500, 1000 and 2000 images. Furthermore, we evaluated models trained on a specific dose on the test sets of other doses. For example, the 5Gy model was evaluated on the 1 Gy and 2 Gy test sets. This lets us see how much a change in dose affects model performance. This performance evaluation is summarized in Figure 3.9.

Because this is a simulation study intended to study the scaling behaviour of Faster R-CNN in rare object detection, we did not compute the counting distributions produced by our models. Instead, we computed a conventional object detection metric called Average Precision (AP), see Section 4.3.2 and[101] for detailed explanations on AP computation.

Overall, the conclusions are fairly intuitive: performance improves with the number of images in the training dataset. Cross-evaluations reveal that a model trained on a specific dose does not perform as well on other dose levels. Models always perform better in terms of AP at higher dose levels, where the proportion of DCs is higher.

We also note that training models at lower doses is difficult, because DCs become even
more under-represented as the dose get lower. At 1 Gy, there might only be one DC every
20 metaphases, which means that the model does not get any gradient feedback for most
images. We found that training models at doses lower than 1 Gy was unfeasible, as the
trained model would not predict any bounding boxes. Another important detail is that
<figure><img src="image_55.png" /><figcaption>Figure 3.9 – Performance curves on simulated data, for models trained at 1, 2 and 5 Gy, with training dataset sizes ranging between 100 and 2000 images. x-axis is number of images, y-axis is Average Precision (AP), see Section 4.3.2 for an explanation of this object detection metric.</figcaption></figure>

at low doses (less than 0.5 Gy), small datasets might not contain any aberration, making
training completely impossible.

**3.5 Conclusion**

Training classifiers on datasets of chromosome patches provided limited result. Avoiding label leakage and building a truly representative datasets are very difficult tasks. This is further confirmed by the review presented in Chapter 2, which suggests that building
a two-stage ADS system requires a significant amount of engineering to deal with edge
cases.

Training a DL-based object detection model to detect chromosomal aberrations seems much easier at first glance, as those can be trained end-to-end, and have been shown to reach a high level of performance in a wide variety of domains. However, a large labelled dataset is required. We used our classification dataset as a training dataset for image generation methods but the performance was unsatisfying. Both VAEs and GANs showcased low diversity, and low sample quality.

In the end, we used a "skeletonized" variant of this classification dataset to produce a supervised metaphase-level dataset. This was used to run a simulation study on the scaling performance of Faster R-CNN for DC detection. Faster R-CNN reached a non-trivial level of performance using as little as 1000 annotated images, without any architecture search or hyperparameter optimization. However, this result depends strongly on the dose of the dataset the model is trained on. Training on representative lower dose datasets performs worse, and requires a larger number of images. Furthermore, models trained on a specific dose perform differently on other doses, because of the difference in dicentric yields.

While our simulated data is much simpler than real world metaphases, this simulation
study alleviated the largest concern, which was that DL-based object detectors would not
be trainable with an attainable number of labelled images.
# END-TO-END CHROMOSOMAL ABERRATION DETECTION IN METAPHASE IMAGES

As we saw in Chapters 2 and 3, designing a "pipeline" ADS system is difficult. The performance of every step is linked, sometimes in unexpected ways. ADS systems relying on conventional computer vision techniques may experience sudden performance degradation because of minute variations in image quality, as shown in the papers discussing DCScore performance in 2.6. To prevent this, ADCI relies on sophisticated image selection models, as explained in Section 2.4. On the other hand, DL-based models can be trained to be robust to changes in image quality. Our goal in this chapter is to use this ability to build simpler DC detection models.

We use a labelled object detection dataset to build a simpler, more effective model taking inspiration from keypoint regression. Our model predicts Gaussian "heatmaps" composed of circular spots localized in the center of the chromosomal aberration bounding box. Predicting one heatmap per class makes dealing with several aberration types easy. This model performs localization and classification in a single step and reaches a high level of performance, both in terms of object detection metrics and calibration curve fitting. We improve performance and provide a simple way to model prediction uncertainty by considering an *ensemble* of Unets.

The rest of this Chapter is organized as follows. First, we present some of our inspirations for our ensemble-based architecture. Then, we describe our training procedure, the evaluation of our model, and its object detection performance. We discuss the high requirements of biological dosimetry, and why simple object detection metrics are not necessarily informative in this context. Using domain knowledge, we improve the false positives performance of our ensemble, and present state-of-the-art results for fully-automated dose estimation. Finally, we provide visual explanations for some key abilities
of our model, namely its ability to reject non-chromosome objects, and the uncertainty
modelling brought by ensembling.

**4.1 Related works**

**4.1.1 Key-point regression in deep learning**

In keypoint regression, Gaussian spots are predicted over specific landmarks of the image, like the eyes on a human face or the joints of a skeleton. It is a subtask of a large number of computer vision tasks like facial recognition or pose estimation [102], [103]. It can be solved with common image segmentation models by minimizing the $L_2$ distance between a predicted heatmap and the corresponding ground truth over a training dataset. While heatmap regression and object detection are different tasks at first glance, and heatmap regression does not predict bounding box extent, several authors have proposed a unified method that consists in predicting a center point and bounding box dimensions [67], or the corners of a bounding box as key-points [69] (see Section 1.4.2 for a detailed discussion). To our knowledge, key-point regression has never been used in biological dosimetry.

**4.1.2 Object detection and counting**

Object counting is a common computer vision task, and a wide variety of solutions
have been suggested in the literature. Density-based methods aim to predict counts by
integrating a density map [104], [105]. This class of methods is usually simple to imple-
ment and reaches a high level of performance. However, it does not detect objects, and
summary statistics like Precision and Recall cannot be computed. Alternatively, counting
can be solved as a subtask of object detection, by enumerating the bounding boxes be-
longing to a certain class. However, accurately predicting bounding boxes is more difficult
than predicting a density map, and detection-based methods usually perform worse than
density-based ones [106]. In [107], the authors propose a detection-based method that
relies on predicting a Gaussian spot centered on the object, but does not predict bound-
ing box extent. In [69], the author propose a model that uses the same heatmap-based
technique, but also predicts a (height, width) tuple for a bounding box.
4.1.3 Model aggregation

A large number of papers have studied ensemble methods for neural networks, either to improve model calibration, for uncertainty modelling or to improve performance. The authors focus either on sampling a diverse set of models to build an ensemble, or on the properties of the aggregated prediction.

For classification models, Lakshminarayanan *et al.* [57] showed that ensembles of neural networks improved on the performance and calibration of single models. As an ensemble of *M* models requires *M* training runs, [58] showed that a carefully chosen learning rate schedule could encourage loss landscape exploration to get a collection of models with high diversity in a single training run by retrieving checkpoints. In [108], samples and checkpoints are re-weighted according to their performance, like AdaBoost [109].

In semantic segmentation, ensemble methods have received attention because they
improve performance and provide a localized measure of uncertainty, usually by computing
the entropy of the ensemble average for every pixel in the image. In [110], an ensemble of
fully convolutional neural networks is used to segment aerial images. In [111], a diverse
ensemble is built by training several Unets [10] with different encoders, and predictions
are aggregated with a weighted average. In [112], the authors provide a review of ensemble
methods for polyp segmentation.

For bounding-box based detection models, aggregation is required, as several boxes
localizing the same object may overlap. Non-Maximum Suppression (NMS) [113] is then
used and consists in sorting the boxes with respect to their confidence levels. Lower
confidence boxes overlapping a high confidence box beyond a specific IoU (Jaccard index)
threshold are discarded. In [114], box merging algorithms are explicitly considered in a
particular model aggregation framework. The authors suggest computing an "average"
box by weighting coordinates based on the box confidence.

Approximate Bayesian deep learning

Common techniques used to sample the posterior distribution are intractable for modern neural networks given their large parameter counts. In [7], Mandt *et al.* demonstrate that SGD can be seen as an Orstein-Uhlenbeck process with some limit Gaussian distribution. SGD can be seen as Langevin sampling of the posterior weight distribution, and SGD samples can be used to estimate the mean and covariance of the Gaussian posterior distribution. This provides a cheap and simple way to retrieve samples of this distribution
during the training procedure.

in [59], Garipov et al. al show that local minima are connected by low-loss paths, and that one could average weight vectors along deterministic trajectories between those local minima to improve performance. In [13], Izmailov et al. confirm the theoretical analysis of [7] by showing that averaging SGD iterates leads to wider minima and improved generalization in practice.

Using the theoretical insights explained in [7], Maddox et al. [6] approximate the limit posterior weight distribution with a Gaussian distribution, where the covariance matrix is defined as the covariance of the last gradient descent iterates. New sets of weights can be sampled from this posterior distribution for ensemble and uncertainty estimation. In [115], the authors go further by sampling an ensemble from several modes of the posterior weight distribution to increase ensemble diversity and therefore performance.

Finally, other approximations of the posterior weight distribution have been proposed. In [116], the authors use Kalman filtering to derive a sequential estimate of the posterior distribution over the weights.

## 4.2 Methods

In this section, we first give a precise definition of the data and model, including the loss for training. Second, we explain our aggregation procedure, going from a set of continuous heatmap predictions for a single image to a set of binary decisions maps and finally, to an aggregated ensemble-level prediction. Finally, we explain the PCA-based visualization techniques used to justify the performance gains of the ensemble and its robustness to the presence of debris in metaphase images.

### 4.2.1 Keypoint regression with heatmap regression models

In Heatmap Regression Models (HRMs), objects of interest are represented as Gaussian spots. The model is trained to predict spot positions in the image domain, with a labelled dataset $\mathcal{D} = \{(x_1, y_1), \cdots, (x_n, y_n)\}$ comprised of $n$ realizations of a pair of random variables $(X, Y)$. For each image $x_i$, we have $x_i(u, v) \in [0, 1]$ at each location $(u, v) \in \Omega$, where $\Omega$ denotes the image grid of size $|\Omega| = H \times W$.
<figure><img src="image_56.png" /><figcaption>Figure 4.1 – Visual comparison between regularized and unregularized model. First image from the right: input image, second: bottom right crop, third: gradient norm of the prediction for the unregularized model, fourth: gradient norm of the prediction for the regularized model.</figcaption></figure>

## Model design and sparsity promoting loss function

Our heatmap regression model is a convolutional neural network $\phi_{\theta}(x) : x \in [0, 1]^{H \times W} \rightarrow y \in [0, 1]^{H \times W}$. The final output is constrained between 0 and 1 with a sigmoid layer. We use the Unet architecture [10] to predict a low resolution heatmap. The image $y_L$ is of size $(H/L, W/L)$ for some arbitrary downsampling factor $L$, as the location accuracy provided by the highest resolution output is not useful. For Unet-based architectures, images are usually downsampled (or upsampled, in the decoder) by a factor of 2 at each layer: at layer $l$ of the encoder, features have a spatial size of $H/2^l \times W/2^l$. For the sake of simplicity, notations are given in the single-channel case. Additional classes of aberrations (like fragments) are modeled with additional channels, so that a third index is added. In the remainder of this Chapter, we consider two aberration classes, DCs and fragments. The parameters $\theta$ are learned by solving the following optimization problem:

$$ \theta^* = \arg \min_{\theta} \sum_{i=1}^{n} \mathcal{L}(\phi_{\theta}(x_i), y_i) + \lambda \mathcal{R}(\phi_{\theta}(x_i)), \quad (4.1) $$

where $\lambda$ is an hyperparameter that balances the data fidelity term and the regularization term. In our modelling approach, the data fidelity term $\mathcal{L}$ has the following form:

$$ \mathcal{L}(\phi_{\theta}(x_i), y_i) = \|\phi_{\theta}(x_i) - y_i\|_2^2 \quad (4.2) $$

Because the number of aberrations is very low compared to the number of pixels in the image, the background is expected to be 0, except in a small number of "hot" spots corresponding to locations containing aberrations. The Sparse Variation (SV) regularizer
<figure><img src="image_57.png" /><figcaption>Figure 4.2 – For each image, we compute the total variation of its prediction for the regularized and unregularized model, as shown in Figure 4.1. This figure represents the Cumulative Distribution Function (cdf) of the total variation over all images in the test set.</figcaption></figure>

(4.3) has been specifically considered here to encourage the emergence of a very small number of "hot" spots as aberrations are rare events in Giemsa-stained images. This regularizer is defined as [11]:

$$ \mathcal{R}(\phi_{\theta}(x_i)) = \sum_{(u,v) \in \Omega} \sqrt{\rho^2 \|\nabla_{u,v}\phi_{\theta}(x_i)\|_2^2 + (1-\rho)^2 \phi_{\theta}(x_i)^2(u, v)}, \quad (4.3) $$

where $\rho$ is a parameter that balances the sparsity and the smoothness terms in the predicted heatmap. The components of the gradient vector are computed with respect to the image coordinate axes as follows:

$$ \nabla_{u,v}\phi_{\theta}(x_i) = \begin{bmatrix} \phi_{\theta}(x_i)(u, v) - \phi_{\theta}(x_i)(u + 1, v) \\ \phi_{\theta}(x_i)(u, v) - \phi_{\theta}(x_i)(u, v + 1) \end{bmatrix}. \quad (4.4) $$

As this is a simple linear transformation of the image, computational overhead is minimal.

The criterion (4.2) is highly non-convex because of non-linearities in $\phi_{\theta}$. Therefore, finding a global minimum is hopeless. Nevertheless, a good local minima may be found
using iterative first order methods, usually some variant of SGD. The exact gradient of the training criterion with respect to $\theta$ is estimated on a random subset $\mathcal{J}$ of the complete dataset $\mathcal{D}$, because of memory constraints:

$$ \nabla_{\theta} \sum_{i=1}^{n} \mathcal{L}(\phi_{\theta}(x_i), y_i) + \mathcal{R}(\phi_{\theta}(x_i)) \approx \nabla_{\theta} \sum_{j=1}^{|\mathcal{J}|} \mathcal{L}(\phi_{\theta}(x_j), y_j) + \mathcal{R}(\phi_{\theta}(x_j)). \quad (4.5) $$

In Figure 4.2, we reported two curves corresponding the cumulative distributions functions of the total variation images computed over the prediction maps regularized with the sparse variation regularizer (blue curve) and without (red curve). See 4.1 for an example of images on which this total variation is computed.

### 4.2.2 Implicit ensembling of neural networks

<figure><img src="image_58.png" /><figcaption>Figure 4.3 – Vote-based aggregation of checkpoints. DCs predictions are plotted in red, and fragment predictions are plotted in green. For every image $x_i$, the heatmap prediction $\phi_{\theta}(x_i)$ is binarized (giving $\Lambda_{\theta}(x_i)$) with a confidence threshold $T_C$. Those maps are summed (giving $S_i$), and regions of the image receiving more than $T_A$ votes are considered as detections. Darker shades of red and green indicates region of the images receiving more votes.</figcaption></figure>

To reduce the number of false positives (i.e, improve Precision) at a fixed recall level, we investigated the ensemble method introduced in [57]. Because training a deep neural network is a stochastic process, successive training runs of the same model tend to explore different regions of the parameter space. Those local minimas are usually very close in terms of validation loss, but their predictions are not identical. This behavior has received significant attention in the literature [58], [59]. As shown in [58], it may not even be necessary to run several successive training runs. A carefully chosen learning rate schedule
may be enough to achieve enough parameter space exploration to build a diverse ensemble from checkpoints of a single training run. In our case, we even find that the gradient noise introduced by stochastic batch sampling leads to sufficient checkpoint diversity for aggregation to be worthwile without any specific learning rate schedule. Therefore, we do not have to deal with the training instability mentioned in [6].

More formally, SGD can be seen as a Langevin sampling of the posterior weight distribution $p(\theta | \mathcal{D})$ [7]. To avoid excessive autocorrelation between samples, we store the vectors $\theta$ at the end of each epoch, instead of every gradient step. This can be interpreted as a variant of thinning, also used in Monte Carlo Markov Chain inference. For each image $x_i$, $i \in \{1, \dots, n\}$ in the test set, we consider a set $\{\phi_{\theta_1}(x_i), \dots, \phi_{\theta_M}(x_i)\}$ of predictions, where $M$ is the number of predictions (or models). Using a confidence threshold $T_C$, we build a set $\{\Lambda_{\theta_1}(x_i), \dots, \Lambda_{\theta_M}(x_i)\}$ of $M$ different binary predictions for each test image $x_i$, that we sum over all members of the ensemble at each location $(u, v) \in \Omega$ as follows:

$$S_i(u, v) = \sum_{m=1}^{M} \Lambda_{\theta_m}(x_i)(u, v), \quad (u, v) \in \Omega. \qquad (4.6)$$

Finally, we set an agreement threshold $T_A$ to compute the agreement between the artificial "experts". The setting of thresholds $T_A$ and $T_C$ impact the final decision. If the confidence threshold $T_C$ is high and the voting threshold $T_A$ is low, the decision will be made from a small set of experts. Otherwise, a small value of $T_C$ but a high voting threshold $T_A$ means that low confidence predictions are considered, but a higher agreement between them is needed to confirm a detection. In the end, we get a precision surface depending on $T_C$ and $T_A$. Our aggregated decision for any image $x_i$ is a binary image $D_i$ such that value at location $(u, v)$ is 0 if no aberration is predicted, and 1 otherwise (See Figure 4.3 for illustration):

$$D_i(u, v) = 1[\mathbf{S}_i(u, v) > T_A], \quad (u, v) \in \Omega. \qquad (4.7)$$

The agreement threshold $T_A$ can be adjusted by the end user to optimize either Precision or Recall scores, like the confidence threshold $T_C$. We discuss the effects of choosing a specific threshold in Section 4.3.2.

### 4.2.3 Setting of model parameters

We trained Unet for $N_e = 100$ epochs with Adam [12], with a constant learning rate of $3 \times 10^{-4}$, a weight decay parameter of 0.1 and a batch size of 12 on a single
Tesla V100. The learning rate was unchanged during training to ensure parameter space exploration, using an analogous reasoning to the one provided in [13].

We did not use data augmentation for two reasons. First, we found that the wide variety of chromosome morphology and orientations in our dataset was enough for our model to learn this invariance. We did not observe detection failures based on object orientation. Second, more agressive data augmentation like noise or blurring quickly made monocentrics and dicentrics indistinguishable. Training stability was very sensitive to the variance of the blurring kernel or the Gaussian noise, because accurate chromosome classification relies on very small details.

We predict a lower resolution heatmap of size $H' = 224, W' = 252$, where height and width are downsampled by a factor of 4. While batch sampling (and therefore parameter space exploration) is randomized, parameter initialization is fixed between training runs. We ran a grid search with $log_{10}$ spacing for $\lambda$ regularization parameter with 10 and $10^{-4}$ as upper and lower bound of the search interval. Training was implemented in PyTorch [20], and uses `segmentation_modelsPytorch` implementation of Unet.

### 4.2.4 Visualization of the training dynamics of single model

As an additional visual explanation for aggregation performance gain, it may be informative to display training trajectories in feature space. Classification networks output a single classification vector per image, so that plotting training dynamics over time using dimensionality reduction is relatively easy (see [117]). A scatterplot of classification vectors embedded in a lower dimension at each epoch provides a good view of how classes are progressively separated during training. Usually, UMAP [118] is used, which requires a pairwise distance matrix between classification vectors. However, this visualization does not work for models outputting a probability distribution for all locations $u, v$ of the input image, like Unet.

To address this issue, we adapt the approach [117] to our context. We consider feature maps as bags of independent feature vectors. Figure 4.4 provides a visual summary of our approach. We do not retain feature vectors for all locations $u, v$ in the feature map. Instead, we only select feature vectors corresponding to the locations of aberrations, and retrieve some feature vectors at random 'background' (i.e., where there are no aberrations) locations. This bag of features is projected on the 2D plane using its PCA decomposition. By retrieving the same locations across several training steps, we can visualize how the both aberration classes and the background are separated during training. An SVM clas-
<figure><img src="image_59.png" /><figcaption>Figure 4.4 – Procedure used to display feature separation in the latent space of the last decoder block for a single epoch (i.e a single weight vector $\theta$). For all images $x_1, \dots, x_n$, the feature maps produced by the last layer of the decoder are retrieved and treated as a set of independent, $C_l$-dimensional feature vectors. Using PCA dimension reduction, we produce a 2D scatterplot that shows how the model separates the different classes (background, dicentrics, fragments) across training epochs. Note that the eigenvectors used for this dimensionality reduction are computed over all epochs of training.</figcaption></figure>

sifier fitted on the embeddings retrieved for a single epoch is used to map regions of the latent space to a specific class, which helps visualize the dynamics of training. The rest of this section gives a formal overview of our visualization technique.

Formally, for an input image $x_i$ of size $H \times W$, the $\ell$-th layer of our Unet produces a feature volume $f_i^\ell$ of size $H' \times W' \times N_\ell$ (see Section 4.4.3 and Figure 4.4 (a)). Here, we choose the second-to-last layer of the decoder and we drop the superscript $\ell$ to improve readability, so that $f_i^\ell = f_i$. For the last layer, $N_\ell = 128$ (see illustration in Figure 4.4). We retrieve the set of feature vectors that correspond to the spatial locations of aberrations (dicentrics and fragments) in image $x_i$. For a set of $N_a^i$ aberrations located at $\mathcal{X}_a^i = \{(u_1, v_1), \dots, (u_{N_a^i}, v_{N_a^i})\}$ in image $x_i$, we build a set of feature vectors
$f_a^i = \{f_i(u_1, v_1), \cdots, f_i(u_{N_a^i}, v_{N_a^i})\} \in \mathbb{R}^{N_a^i \times N_e}$. We retrieve the feature vectors corresponding to all aberrations in every image the test set. We also sample an additional set of $N_b^i$ background pixels, denoted as $f_b^i$ at locations $\mathcal{X}_b^i = \{(u'_1, v'_1), \cdots, (u'_{N_b^i}, v'_{N_b^i})\}$. Those locations are randomly sampled, provided the locations do not correspond to aberration pixels. Therefore, they correspond to background, monocentric or debris. Finally, we define $f_{a,b}^i = f_a^i \cup f_b^i$ so that the total number of feature vectors in $f_i$ is $N_a^i + N_b^i$, as shown in Figure 4.4b). It is worth noting that $f_{a,b}^i$ is a subset of the complete feature map $f_i$. This makes the visualizations described less cluttered, and reduces computation time. Finally, this feature set $f_{a,b}^i$ is retrieved for each image $x_i$ in the test set, to build a large feature set $\mathcal{F}_e = \{f_{a,b}^0 \cdots f_{a,b}^n\}$, where $e$ corresponds to the set of model parameters retrieved at epoch $e$.

The set of feature vector $\mathcal{F}_e$ is retrieved for each epoch $e \in \{1, \cdots, N_e\}$. These sets are concatenated in global set $\mathcal{F} = \{\mathcal{F}_1, \cdots, \mathcal{F}_{N_e}\}$. Although all the subsets of $\mathcal{F}$ correspond to the same locations, they are different at each epoch $e$ because of the stochasticity of gradient descent. PCA can be used to visualize those feature sets as 2D scatterplots, and in particular to check how well aberrations are separated from the background at every epoch. The set $\mathcal{F}_1$ is chosen as a reference feature set and $\mathcal{F}_e, e \in \{2, \cdots, N_e\}$ is registered with respect to $\mathcal{F}_1$ using the Procrustes method [119]. This guarantees that latent space scale shifts or rotations are removed for visualization. A PCA decomposition is computed on $\mathcal{F}$, and for each epoch $e$, each feature set in $\mathcal{F}_e$ is projected on the first two principal components of this decomposition. This provides a 2D visualization of the trajectory of each feature vectors during training.

Furthermore, we train a kernel SVM classifier $p_e$ on the 2D embeddings of $\mathcal{F}_e$ to predict which aberration class corresponds to a location in 2D embedding space. This classifier takes a 2D embedding (a vector of $\mathcal{F}_e$) as an input, and outputs a probability distribution over three classes: background, DCs and fragments. For all epochs $1 \le e \le N_e$, we train a different classifier, and predict a probability distribution over a grid that samples the 2D aberration space uniformly. For all positions $(u, v)$ of this grid, a probability distribution over the total number $N_r$ aberration classes $p_e(u, v, r) \in [0, 1], r \in \{0, \cdots, N_r\}, \sum_{r=1}^{N_r} p_e(u, v, r) = 1$ is predicted at epoch $e$. As all the point clouds are aligned and a single set of principal components is computed for all time steps, the changes in the decision boundary from one epoch to the next can be solely attributed to the dynamics of training.

Finally, to visualize the displacement of class boundaries across training, we define the
averaged classifier:

$$ \bar{p}(u, v) = \frac{1}{N_e} \sum_{e=1}^{N_e} p_e(u, v). \qquad (4.8) $$

Visualizing the spread of the distribution of $\bar{p}$ can be done by computing the entropy of the distribution predicted by the averaged classifier:

$$ H(\bar{p})(u, v) = - \sum_{r=1}^{N_r} \bar{p}(u, v, r) \log \bar{p}(u, v, r). \qquad (4.9) $$

## 4.3 Materials

### 4.3.1 Data description

<figure><img src="image_60.png" /><figcaption>Figure 4.5 – Repartition of images into aberration counts bins.</figcaption></figure>

Our training dataset is composed of 5430 labelled images of size $H = 888, W = 1008$, padded to $H = 896, W = 1008$ to ensure that downscaling has an integer height and width. Labels are binary images with size $H' = 202, W' = 252$, taking value 0 everywhere except at the center of chromosomal aberrations (roughly between the two centromeres for a DC), where it takes value 1. There is one binary image per aberration classes for each image, so that aberration classification is possible. Chromosomal aberrations are
<figure><img src="image_61.png" /><figcaption>Figure 4.6 – Sketch of model evaluation. The intersection $I_i$ between the binarized ground truth $\Lambda_{GT}(y_i)$ and the binarized prediction map $\Lambda_{\theta}(x_i)$ is computed. Objects appearing in both are true positives, objects appearing only in $\Lambda_{\theta}(x_i)$ are false positives, objects appearing only in $\Lambda_{GT}(y_i)$ are false negatives. In this case, we have two true positives, 1 false negative and 1 false positive, so that Precision is $TP/(TP + TP) = 2/3$ and Recall is $TP/(TP + FN) = 2/3$</figcaption></figure>

the only labelled objects, neither debris nor monocentric chromosomes are labelled. We chose this labelling scheme instead of semantic segmentation or bounding boxes as it lead to the lowest labelling time, which in turn meant a greater number of images could be labelled for the same labelling budget. For the same reason, we did not label debris or monocentric chromosomes. This also prevented the discovery of trivial models where chromosomes would be detected but always labelled as monocentric, as they outnumber dicentric ones by an extremely large margin. As explained in Section 4.2.1, this binary image is blurred with a Gaussian kernel, to reduce the underrepresentation of the labels against the background.

Images have been selected so that each image contains only the chromosomes corresponding to a single cell, i.e., no image contains more (or less) than 46 chromosomes. Our dataset contains 5021 dicentrics and 7540 fragments, Figure 4.5 shows the repartition of images of images into aberration counts bins. Images with a high aberration count are
much rarer than images with a low aberration count. This training dataset is a subset of a much larger archive of patient data (containing around 80k images). It contains all images with at least one aberration in the archive. This reduces training time considerably, and prevents the discovery of trivial models where no object is ever predicted. However, this also means that it is not an accurate representation of real-world metaphase. As a normal metaphase contains 23 chromosome pairs, this means that even in this case, the overwhelming majority of chromosomes are healthy (i.e monocentric) ones. In our evaluation setting, the training set consists of 80% of those image. 10% of the images are retained for a validation set, used for hyperparameter selection. Finally, 10% of the data is held as a test set for a fair performance evaluation.

As this training dataset is not representative of a real-world exposition, we use another dataset to estimate the calibration curve of our model. This dataset was built by collecting metaphases from samples irradiated at specific, known doses. The aberrations in this dataset are not labelled. It contains 21215 metaphases taken from samples irradiated at 0 Gy, 0.1 Gy, 0.2 Gy, 0.3 Gy, 0.5 Gy, 0.7 Gy, 0.9 Gy, 1 Gy, 1.5 Gy, 2 Gy, 3 Gy and 4 Gy.

### 4.3.2 Evaluation metrics

#### Performance of a single model

While our predicted heatmap $\phi_{\theta}(x_i)$ can take any value between 0 and 1 at each spatial position $(u, v)$ in the image, ultimately a binary decision needs to be taken with regard to the presence or absence of aberration at location $(u, v) \in \Omega$. In the next step $\phi_{\theta}(x_i)$ is used to build a binary map $\Lambda_{\theta}(x_i)$ given an arbitrary threshold $T_C$:

$$ \Lambda_{\theta}(x_i)(u, v) = 1[\phi_{\theta}(x_i)(u, v) > T_C], \quad (u, v) \in \Omega, x_i \in \{x_1, \dots, x_n\}, \qquad (4.10) $$

where $1[\cdot]$ is the indicator function. Furthermore, a ground truth binary map is defined given a confidence threshold $T_{GT}$, which is set to a small fixed value (e.g., 0.01).

$$ \Lambda_{GT}(y_i)(u, v) = 1[y_i(u, v) > T_{GT}], \quad (u, v) \in \Omega, y_i \in \{y_1, \dots, y_n\}. \qquad (4.11) $$

Once a binary decision for the presence of an object is made at each location in the image, we define true positives, false positives and false negatives. In object detection, true positives are usually defined up to a small location error, as matching the ground truth perfectly would be too stringent. In our case, the spots are small compared to the
object size so that the position error remains very small even in the cases where the intersection between the predicted and ground truth spot is the smallest possible one (one pixel) (see Figure 4.6). Therefore, we consider any overlap between a prediction and a ground truth spot to be a true positive, as long as this ground truth spot has not been predicted before. If this is the case, the prediction is considered to be a false positive. Predicted objects that do not overlap ground truth spots are also considered to be false positives. Finally, objects in the ground truth heatmap that are not predicted by the model are considered to be false negatives. True negatives are ill-defined in object detection, and are not considered. With those three values, we can compute Precision and Recall.

More formally, the intersection image $I_i$ between $\Lambda_\theta(x_i)$ and $\Lambda_{GT}(y_i)$ is computed as the pixel-wise product of ground truth and prediction

$$I_i(u, v) = \Lambda_\theta(x_i)(u, v) \times \Lambda_{GT}(y_i)(u, v), \quad (u, v) \in \Omega, \ i \in \{1, \dots, n\}. \qquad (4.12)$$

Finally, we compute the number of connected components $N_{cc}(I_i)$, $N_{cc}(\Lambda_\theta(x_i))$ and $N_{cc}(\Lambda_{GT}(y_i))$ in $I_i$, $\Lambda_\theta(x_i)$ and $\Lambda_{GT}(y_i)$, respectively. The number of true positives, false positives and false negatives are defined as follows:

$$\left\{
\begin{aligned}
TP &= N_{cc}(I_i), \\
FP &= \max(0, N_{cc}(\Lambda_\theta(x_i)) - TP), \\
FN &= \max(0, N_{cc}(\Lambda_{GT}(y_i)) - TP).
\end{aligned}
\right.
\qquad (4.13)$$

Hence, we compute Precision = $\frac{TP}{TP+FP}$ and Recall = $\frac{TP}{TP+FN}$ for a single $(\Lambda_\theta(x_i), y_i)$ pair. Note that Precision and Recall are functions of the chosen confidence level $T_C$; a higher confidence threshold increases Precision but decreases Recall. In what follows, we report results for a set of confidence thresholds $\{T_{C_1}, \dots, T_{C_{10}}\}$ to analyze this tradeoff. Moreover, we decided to report the Precision and Recall scores separately instead of providing an aggregated metric like Average Precision (AP), as the tradeoff between those metrics in the context of biological dosimetry is especially important.

Because training is a stochastic process, the maximum predicted probability over the test set is not exactly 1; different models may get a different maximum confidence value. In other words, two sets of weights $\theta$ and $\theta'$ will give two different maximum probabilities $p_\theta$ and $p_{\theta'}$, so that for each model, the performance metrics are computed over different confidence thresholds. Therefore, each performance metric curve is linearly interpolated
over a common confidence grid {0.1, 0.2, ..., 0.9}. For the set of confidence thresholds that exceed the maximum probability over the complete test set, Precision is not defined. In this case, it is arbitrarily chosen to be 1 (and Recall is 0). To provide a metric showcasing performance variation across training, we reported performance quantiles (5% and 95%) for each threshold. This means that, for each threshold, 50 Precision and Recall values are computed (one for every considered epoch) and the aforementioned quantiles of those values are reported.

## Performance of model ensemble

For the ensemble, we use the same evaluation procedure as in the single model case (described in Section 4.3.2), but with the aggregated binary decision map $D_i$ described in Section 4.2.2. However, as explained in Section 4.2.2, the performance of the ensemble depends on an agreement threshold $T_A$ and a confidence threshold $T_C$. Therefore, instead of a Precision curve, we get a Precision surface, which makes comparison with the single model case more difficult. Instead, we set a specific vote threshold, and report the same Precision curve as in the single model case. Finally, to evaluate the sensitivity of the performance to the sampling of the ensemble, we sampled 100 random ensembles, and computed $q_{05}$ and $q_{95}$ for every confidence threshold in the grid mentioned in the Section 4.3.2.

# 4.4 Experimental results

In this section, we discuss model performance, both in term of object detection and calibration curve estimation. First, we discuss the results of the single-model training and the impact of regularization term. We also show a visualization which suggests that our model is robust to the presence of debris, without the need for specific labelling or architectural choices. Finally, we discuss the performance improvements obtained with the ensemble approach, and we provide PCA-based approach to visualize model feature trajectories during training.

## 4.4.1 Performance of single model

In Table 4.2, the single non-regularized model already achieves significant gain in terms of Precision and Recall when compared to Metafer in the case of MC versus DC classifica-
<figure><img src="image_62.png" /><figcaption>Figure 4.7 – Precision, Recall and False Discovery Rate (FDR) as functions of confidence for DCs (left) and fragments (right). Top: performance summary for the unregularized model (i.e $\lambda = 0$ for the sparse variation term). Middle: performance summary for $\lambda = 0.2, \rho = 0.1$. Bottom: performance summary for the ensemble of checkpoints from the training of the regularized model for a threshold of 2 votes. Shaded area indicates the [$q_{0.05}, q_{0.95}$] inter-quantile interval, computed respectively over the last 50 checkpoints for single models, and over a 100 random samples of 10 checkpoints for the bottom plot (ensemble).</figcaption></figure>
tion (Metafer is not designed to detect fragments), although some performance variation is noticeable during training, as confirmed by the inter-quantile range of performance. This variation suggests that there is sufficient parameter space exploration to get enough prediction diversity for aggregation to be worthwhile.

Metafer relies on conventional computer vision techniques. The chromosome objects in the metaphase are segmented and then classified. Segmentation errors can lead to classification errors, for example when one chromosome is split during segmentation. Overall, the performance of both tasks (segmentation and classification) is not very robust to variations in image quality induced by variations in acquisition circumstances like illumination or staining quality. For example, for the segmentation task, the Recall of Metafer lies between 35% and 75%. For the classification task, The Recall is around 35%, and the Precision around 40%. Note that the performance figures given for Metafer in Table 4.2 are not computed on the dataset mentioned in Section 4.3.1, but taken from [93] instead. Although the performance of Metafer is uncertain, and we did not compute it on our dataset, we feel confident that our model brings a very significant improvement in chromosomal aberration, as the uncertainty around the ensemble performance is low enough that even in the worst case scenario, it achieves very significant performance improvements over Metafer.

Finally, training with a loss that promotes sparsity improves Precision and Recall for fragments and DCs. this improvement in precision is especially relevant for automated use, as keeping the number of false positives low is required to avoid overloading care facilities. Because all chromosomes may not be correctly retrieved (Recall is less than 1), this model tends to underestimate doses.

<table><caption>Table 4.1 – Comparison between model 1 (denoted as M1, $T_C = 0.6, \lambda = 0$), model 2 (denoted as M2, $T_C = 0.6, \lambda = 0.2, \rho = 0.1$) and ensemble (4 votes, $T_C = 0.5$) for the fragment class. [$q_{5\%}$, $q_{95\%}$] interval is reported in parenthesis. For model 1 and model 2, this interval is computed over the last 50 epochs of training. For the ensemble, it is computed over a 100 randomly sampled ensembles (ensembles are sampled from checkpoints during training). All performance metrics are computed on the test set.</caption><thead><tr><th></th><th>M1 (unregularized)</th><th>M2 (regularized)</th><th>ensemble (3 votes, T<sub>C</sub> = 0.5)</th></tr></thead><tbody><tr><th>Pr (frags)</th><td>70.8% (62.0,77.4)</td><td>79.3% (71.8, 84.2)</td><td>81.0% (78.5, 83.3)</td></tr><tr><th>Rc (frags)</th><td>49.4% (37.8, 59.7)</td><td>55.0% (46.0, 66.3)</td><td>65.5% (62.4, 68.5)</td></tr></tbody></table>
<table><thead><tr><th></th><th>M1 (unregularized)</th><th>M2 (regularized)</th><th>ensemble (3 votes, T<sub>C</sub> = 0.5)</th><th>Metafer</th></tr></thead><tbody><tr><td>Pr (dics)</td><td>76.6% (68.8, 84.7)</td><td>83.6% (75.4, 92.2)</td><td><b>85.8%</b> (83.7, 88.3)</td><td>~ 40%</td></tr><tr><td>Rc (dics)</td><td>45.2% (31.6, 56.4)</td><td>51.2% (37.0, 62.7)</td><td><b>61.7%</b> (57.1, 65.8)</td><td>~ 35%</td></tr></tbody></table>

Table 4.2 – Comparison between model 1 (denoted as M1, $T_C = 0.6, \lambda = 0$), model 2 (denoted as M2, $T_c = 0.6, \lambda = 0.2, \rho = 0.1$), ensemble (4 votes, $T_C = 0.5$), and Metafer for the DC class (where Metafer performance is available). All performance metrics are computed on a separate test set. Metafer performance is retrieved from previous work, and was not evaluated on the dataset used in this Chapter. Metafer performance should only be taken as a rough point of reference, see Section 4.4.1 for additional discussion of this point.

### 4.4.2 Performance of model ensemble

<figure><img src="image_63.png" /></figure>

Figure 4.8 – Prediction diversity for a set of four different models. The first model does not predict the fragment in the bottom left corner (a very low confidence threshold would be required to consider it as a detection), but the three next model do. All four models predict the DC near the center of the image correctly.

In this section, we report the results obtained with the ensemble procedure. The parameter space of the model is explored through SGD. Although the selected checkpoints reach a similar validation performance, the predictions they yield vary, as shown by Figure 4.8. As mentioned earlier, this can be used to filter spurious predictions, as those are less likely to be present in all models of the ensemble.

Precision is a monotonously increasing function of voting and confidence thresholds, while Recall decreases with higher confidence and higher voting thresholds. End-users may fine-tune the balance between both of those metrics depending on their goal by choosing a specific ($T_C, T_A$) pair. To estimate the sensitivity of Precision and Recall to the sampling of the checkpoints, we evaluate those metrics for 100 samples ensembles and report the $q_{05}, q_{95}$ interval for Precision and Recall in Figure 4.7).
To ensure that ensemble results are easily readable, we do not report surfaces for Precision and Recall. Instead, we select a single voting threshold and report the results over all agreement thresholds, like with the single model results. The ensemble provides a significant performance improvement over the single-model baseline; aggregation does help to filter out spurious detections and improves Precision and Recall, as reported in Tables 4.1 and 4.2. We chose the ($T_C, T_A$) parameters to keep Precision broadly similar across all models, for DCs and fragments. This makes Recall improvements more salient, but one could choose other values for confidence and vote thresholds. Overall, there is a wide set of threshold combinations that yield large performance improvements over the Metafer baseline. Furthermore, ensembling also reduces performance variation: the performance is closer between different ensembles than between single checkpoints. This suggests that our results are not dependent on a specific sampling or selection of the ensemble.

### 4.4.3 Robustness to non-chromosome objects in metaphase images

<figure><img src="image_64.png" /><figcaption>Figure 4.9 – Rejection of nuclei depending on model layer. First image from the left shows the input image and ground truth Gaussian heatmap. Second image shows PCA embedding of features at the output of the first encoder layer. Third image shows PCA embeddings of features at the last encoder layer. Rightmost image shows embeddings of features at the first decoder layer. The embedded feature maps are resized from $H', W'$ to $H, W$ so that every image has the same size.</figcaption></figure>

An automated dosimetry system should always distinguish between chromosome and non-chromosome (nuclei, debris) objects. Current automated methods reject non-chromosome objects using explicit shape analysis. For example, a nuclei can be rejected using the fact that it is broadly circular and has a uniform texture. However, the shape of debris is usually more complex, which makes metaphase selection more difficult in most dosimetry
systems.

In our approach, we do not detect debris and nuclei explicitly, as they are not annotated in our dataset. Instead, Unet learns to reject debris from the training data, without the need for specific annotations or handcrafted object detection algorithms. In the rest of this section, we propose a simple visualization of this fact, using PCA.

For an image $x_i$ and a set of parameters $\theta$, the activation volume provided by the $\ell$-th layer of the neural network is defined as:

$$f_i^\ell = \sigma(\theta_\ell f_i^{\ell-1} + b_\ell), \quad f_i^\ell \in \mathbb{R}^{H' \times W' \times N_\ell}, \qquad (4.14)$$

where $N_\ell$ the number of convolution filters in layer $\ell$, i.e the number of channels of the output of this convolution layer. For a set of $n$ images $\{x_1, \cdots, x_n\}$, the corresponding activation volume is of size $\mathbb{R}^{n \times H' \times W' \times N_\ell}$. This volume can be flattened in a matrix $E_\ell \in \mathbb{R}^{(n \times H' \times W') \times N_\ell}$, containing $n \times H' \times W'$ samples (rows) of a random vector of size $N_\ell$. Once $E_\ell$ is centered and standardized, the eigenvectors of $E_\ell^T E_\ell$ form the usual PCA orthogonal basis $Q_\ell \in \mathbb{R}^{N_\ell \times N_\ell}$, so that $E_\ell Q_\ell^K$ (with $Q_\ell^K \in \mathbb{R}^{N_\ell \times K}$) projects each pixel of the activation volume ($N_\ell$-sized vectors) onto the first $K$ components (columns) of this basis.

By choosing $K = 3$ and normalizing the projected vectors to sum to one, we can plot a visualization of the activations in RGB space (Figure 4.9). Note that a different $Q_\ell$ is computed for each layer considered in Figure 4.9. Therefore, the colors do not have any specific meaning. The second image of Figure 4.9 from the left shows that the projection of feature vectors belonging to chromosomes and nuclei on the PC basis are highly similar, as the color of those regions is identical. The main reason is that the convolution filters in early layers have a small receptive field (see [120]) and mostly capture texture and edge information, which is close for subsets of chromosomes and nuclei. In the early stages of the encoder, non-chromosome objects are not differentiated from chromosomal objects.

The third and last images of Figure 4.9 from the left show that in the deepest layer of the encoder, the feature vectors belonging to nuclei and chromosomes are mapped to different principal components, corresponding respectively to blue and red pixels. By stacking convolutions the model aggregates information from a larger subset of the image (the receptive field increases), which is suitable to distinguish chromosomes from nuclei on the basis of their differences in shape. This is a first indication that our model learned to reject debris and nuclei without domain knowledge, which is further confirmed by our performance results.
#### 4.4.4 Visualization of training trajectories

In this section, we discuss the visualizations produced by the method described in Section 4.2.4. In Figure 4.10, we see the latent space of Unet at 6 different epochs. Red points on the scatterplot represent locations containing DCs, blue points fragments and gray points are background locations. Because those snapshot are taken at regular intervals at the end of training, classes are well separated. We see that the decision boundary of the SVM classifier changes from timestep to timestep. Some locations remain well separated in latent space from others, but this is not true for all samples. Regions in the latent space where classes overlap correspond to areas of uncertainty.

While locations with high prediction uncertainty are hard to detect because neural networks tend to be overconfident [5], an ensemble of models can be used to retrieve this information. This fact is also visible in the latent space of Unet, as shown by Figure 4.11. In this Figure, we performed K-Means classification on our bag of features in the latent space for the last considered epoch in Figure 4.10. Once this is done, we plot the trajectories of the barycenters of those clusters over time. As we see, most cluster centers are well separated across epochs. This is reflected by the fact that if we consider the average of SVM classifiers across epochs, most feature vectors are in low-entropy regions. Feature vectors belonging to high-entropy regions correspond to uncertain detections can be filtered as described in Section 4.2.2.

#### 4.4.5 Transfer between training and calibration curve datasets

As explained in Section 4, most images routinely analyzed by biologists do not contain any aberration, even for high doses. In the previous section, we saw that the model described in this Chapter performs very well in terms of object detection metrics. However, we still need to investigate whether a model trained on this dataset could accurate estimate aberration counts on a more realistic dataset.

Initial calibration curve estimates were unsatisfying: our ensemble would overestimate low doses, and underestimate high doses. Two additional details were needed to improve performance. First, we investigated the Cumulative Distribution Function (CDF) of the maximum probabilities predicted by the members of the ensemble for the DC and fragment class. Those CDFs were estimating using the 4 Gy subset of the calibration curve dataset, and are shown in Figure 4.12. Instead of setting a single threshold for all model, we picked a quantile and retrieved the corresponding CDF value for each model. Fur-
<figure><img src="image_65.png" /><figcaption>Figure 4.10 – Snapshot of the training trajectory in feature space. Each point of every scatterplot represents a fixed location in an image (DC, fragment or background). Because of the stochasticity of training, the corresponding feature vector moves in feature space. The contour map showcases the decision boundary of a kernel SVM classifier trained to predict the type of feature vector depending on its location in feature space.</figcaption></figure>
<figure><img src="image_66.png" /><figcaption>Figure 4.11 – Training trajectories of feature barycenters. The scatterplots displayed in Figure 4.10 are clustered with K-Means to simplify visualization. The trajectories of barycenters during training are displayed in this figure. The thickness of the trajectory shows the number of feature points in the barycenter.</figcaption></figure>

thermore, we used domain knowledge to reject spurious DC detection. We considered DC
detections if and only if at least one fragment was present in the same image. As usual
in biological dosimetry, we fitted a linear-quadratic model to the point cloud of average
DC count retrieved on the calibration curve dataset. This led to a large improvement
in our calibration curve estimation, as shown in Figure 4.13. Furthermore, the Metafer
calibration curve is semi-automated: DCs undergo a manual review, because of the very
high FPR of the algorithm, as described in Section 4.4.1. On the contrary, our calibration
curve is obtained in a fully automated setting.
<figure><img src="image_67.png" /><figcaption>Figure 4.12 – Cumulative Distribution Functions (CDF) of the maximum probabilities predicted by every member of the ensemble for the DC (left) and fragment (right) class over all images corresponding to a 4 Gy dose in the calibration curve dataset.</figcaption></figure>

## 4.5 Discussion and future works

In biological dosimetry, estimating the average number of chromosomal aberrations per peripheral blood lymphocyte metaphase is necessary to estimate an ionizing radiation dose. However, human expertise is required and is therefore a bottleneck to scale chromosome counting beyond a few hundred images per patient.

In this Chapter, we evaluated Unet as an aberration detection model for biological dosimetry. Unet outperformed Metafer (a current commercial solution) in terms of Precision and Recall by a wide margin. Our approach is learning-based and differs significantly from the current state of the art in terms of how much domain knowledge of chromosome morphology is incorporated. We demonstrated a high level of performance without the need for significant shape modelling. Feature visualization suggests that the model learns to reject debris and nuclei in an unsupervised manner, without the need for object-specific annotations for monocentric chromosomes or debris. Furthermore, a simple regularization term modelling the intrinsic heatmap sparsity helps performance.

We pushed this performance further by ensembling several checkpoints collected during training. We proposed a visualization of the latent features of Unet during training to explore the relationship between the dynamics of training and this performance improvement. Furthermore, we showed that the variation in performance between different (randomly sampled) ensembles is lower than between single checkpoints of a training run.
<figure><img src="image_68.png" /><figcaption>Figure 4.13 – Calibration curves estimated by the ensemble. Left: calibration curve before setting a threshold per model and using domain knowledge. Right: calibration curve after model-adaptive thresholding and using domain knowledge. To improve readability, we show the 4 curves closest to the manual calibration curve displayed in black. Metafer curve is displayed in red.</figcaption></figure>

This is especially relevant in the context of the deployment of a deep learning model in
an automated fashion in a medical setting. Those improvements can be achieved without
the need for any architectural modifications or extensive hyperparameter calibration.

Our training dataset is not representative of real world data: the average number of
aberration per cell is over 1, which corresponds to an extremely high dose of ionizing
radiation. It is very likely that some of the spurious detections can be attributed to this
training set imbalance. We evaluated our ensemble of Unet in a realistic setting, on a
calibration curve dataset. Using model-adaptive thresholding and the domain knowledge
of co-occurence of DCs and fragments, we reach a very competitive calibration curve,
widely outperforming the Metafer baseline. Furthermore, our ensemble can be used in
a fully automated fashion, while the Metafer solution requires manual review to reduce
the number of false positives. It is therefore possible to build a competitive aberration
detection system even with a large distribution gap between training and inference images.
# CONCLUSION

## Contributions

The main contribution of this thesis is a proof of concept based on deep learning for a chromosomal aberration detection model in cytogenetic biological dosimetry. This model relies on the Unet architecture. It is trained by minimizing a $L_2$ loss and a sparsity-promoting regularization term over a large, labelled dataset. We use the fact that gradient descent samples can be seen as approximate samples for the posterior distribution of the network to retrieve a diverse ensemble of models during training, and we build an interpretable aggregation function for this ensemble. This function relies on two thresholds, a confidence thresholds and an agreement thresholds, which can be set by the end-user to balance between precision and recall.

In our labelled dataset, all images contain at least one dicentric chromosome, which means that aberrations are much more prevalent than in real life. This oversampling prevents the discovery of trivial models during training, where no aberration is every predicted. However, this also means that there is a distribution gap between a real blood sample and our training dataset. Using some domain knowledge, like the co-occurence of DCs and fragments, we reject a large number of spurious detections, and improve the false positive rate.

Using this domain knowledge, the counts provided at each dose match the manually estimated calibration curve very closely. Furthermore, unlike the previous semi-automated tool used at IRSN where semi-automatic, our solution achieves this false positive rate in a fully automated fashion.

Moreover, we provide a PCA-based visualization of the features of our model. This allows one to display the differences between the features corresponding to chromosomal objects and non-chromosomal objects. This rejection of non-chromosomal objects is achieved without the need for specific modelling or labelling. This contrasts with current commercial solutions that need to reject non-chromosomal objects based on their shape or texture during the chromosome detection step. The deep learning paradigm, where rich image features are learned on large datasets instead of modelled explicitly proved rele-
vant in this specific biomedical imaging application, in spite of the all the issues identified
in Section 2.2.

Furthermore, we used PCA to provide a visualization of the training trajectories of Unet outputs. By plotting the trajectories of pixels corresponding to aberrations in feature space, we are able to evaluate how well our model separates aberrations from background during training. We trained several SVM classifiers, each corresponding to a different epoch, and therefore a different feature space. Visualizing the entropy of the distribution predicted by the average classifier allows one to justify our ensembling gains. Low entropy regions in feature space correspond to aberrations receiving a high number of votes, as most members of the ensemble agree on the object class. On the other hand, high entropy regions correspond to aberrations receiving a lower number of votes.

## Perspectives

While our model reaches a high level of performance, both in terms of object detection metrics and calibration curve estimation, it may not be considered as a complete dose estimation pipeline. Because of variations in image quality some metaphases are completely unusable. Chromosomes may be underspread, or the image may be blurry for accurate chromosome classification. Currently, both of our datasets (training and calibration curve estimation) are comprised of metaphases that are good enough for human scoring. A fully automated dose estimation pipeline needs a metaphase selection step. Metaphase selection with conventional computer vision techniques was explored in [89] to guarantee that every metaphase is complete, i.e., that it contains 23 pairs of chromosomes. More recently, a ConvNet was trained to perform metaphase selection, reaching over 99% accuracy [90]. However, the metaphase selection step does not segment chromosome explicitly and relies on a ResNet-based binary classifier instead.

Metaphase selection could be solved with several architectures mentioned in Chapter 1, like Faster R-CNN [9] or ResNet [8]. To ensure that metaphases are complete, Faster R-CNN is appropriate. However, this means that we need bounding box labels for every single chromosome which is time consuming. In [121], Pachitariu *et al.* proposed an active learning biomedical image segmentation software called Cellpose. The underlying segmentation model is trained on biomedical images, and predictions can be edited to re-train the model. This approach could be used to train a chromosome segmentation model without building a large labelled dataset.
Furthermore, while our training dataset is comprised of 5430 labelled images, this dataset is extracted from a larger database containing around 80 000 images. Our training subset is comprised of all images containing at least 1 aberration. Even if this larger archive does not contain any chromosomal aberration, it still brings a much larger repository of MC and cellular debris examples which would help MC-DC discrimination provided that we can design an effective pre-training task. It has been shown that because conventional similarity-based SSL techniques discard spatial information with global average pooling, they are less effective for pre-training dense models [122]. Therefore, designing pretext tasks for dense prediction models can be considered in future works. However, pretraining models on this large dataset would also be significantly more costly. The increased computation time suggests that the exploration of small architecture changes or new hyperparameter configurations is slowed down considerably, which makes designing such a pretext task more challenging.
<figure><img src="image_69.png" /></figure>
# SPARSITY IN OPTICAL FLOW FOR MICROFLUIDICS IMAGING

From June to September of 2022, I stayed at the Engineering Physics department of McMaster University, in Canada. This research stay was supervised by Qiyin Fang and focused on assisting Tianqi Hong's PhD work by designing methods for optical flow estimation in microfluidics cell imaging. This optical flow estimation method uses the same sparse regularizer as the proof of concept presented in Chapter 4. This method is the subject of this Appendix.

## Abstract

Building labelled optical flow datasets is very challenging, which makes the training and evaluation of those on unseen sequences difficult. To overcome this difficulty, optical flow models are often trained on benchmark datasets composed of synthetic image sequences. This means that during inference on unseen image sequences, there will be a distribution gap between the training and testing data. This distribution gap leads to errors if unseen data does not match the training distribution. This is often the case in biological, medical of satellite images, because those applications showcase unusual object sizes, viewpoints and levels of optical resolution or noise. In this context, a pre-trained, supervised model might fail to estimate a flow field where the flow discontinuities match the object boundaries for small objects. Furthermore, slow objects cannot be continuously detected by those pre-trained models, and they tend to be merged in the background.

The main context for this chapter is shadow imaging in microfluidics applications, where the issues mentioned above are especially relevant. Cells are very small with respect to the field of view, and showcase variations in speed. Because of this, the estimated flow discontinuities do not match cell boundaries, and slow cells are not detected. We evaluate boundary mismatch and speed dependence using the agreement between motion
field segmentation and cell masks as a proxy measurement of flow quality. Moreover, we
propose a novel regularizer to adapt unsupervised optical flow models to microfluidics
imaging.

A.1 Introduction

Lensless microscopy systems relying on cost-efficient Commercial-Off-The-Shelf (COTS) hardware are proposed as an alternative to flow cytometry methods [123]. The image sensor used for acquisition is mounted over the transparent microfluidics channel. The channel is backlit with incoherent lighting. Usually, the considered samples (red blood cells, parasites) have enough transparency for shadow imaging. The optical resolution of the device depends on the distance between the channel and the sensor, but cannot exceed twice the pixel size. Nevertheless, lensless devices have several advantages. They are cheaper than a flow cytometer, which is the current gold standard for fluid analysis in the medical field. They also offer a very wide field of view, which improves the throughput of the device. Figure A.1 provides a schematic of the device implemented in [123]. In return, sophisticated image processing is required to reproduce the functionality of flow cytometers, like cell identification or cell counting. Because cell rotation motion can be used to classify cells, estimating a motion field from pairs of images is a crucial step to build a COTS flow cytometer using microfluidics imaging.

Optical flow estimation was traditionally formulated as an optimization problem [124]–
[127]. Variational methods generally penalized the norm of the gradient of the motion
field to favor piecewise smoothness [124]. For several decades, variational methods were
considered to be the most competitive methods for optical flow estimation [127]. However
the errors generated by the pixel-wise data term are often compensated by an undesirable
over-smoothing. An alternative approach consisted in fitting parametric motion models
in the neighborhood of each pixel [128]. The motion field is estimated from a region-
based data term that exploits the consistency of several pixels in the region. Despite their
superior robustness to noise, these methods are usually unable to compete with global
methods in terms of accuracy at low levels of noise [127], in various application fields (fluid
mechanics [129], cell imaging and microscopy [130]–[133]). The reason is the difficulty to
determine the boundaries of regions where the motion can be accurately approximated by
a parametric motion model. However, it has been shown that when appropriate regions
are chosen, local parametric models can yield good performance [134].
<figure><img src="image_70.png" /><figcaption>Figure A.1 – Schematic of the microfluidics device implemented in [123]. A ~ 80 µm channel is machined in PDMS (gray). A fluid (light blue) carrying objects (cells, parasites) is pumped through the channel for analysis. An image sensor (in green) is bounded directly to the channel. The channel is lit from the opposite side to the image sensor.</figcaption></figure>

More recently, convolutional neural networks have been shown to be very effective
models in a wide variety of computer vision tasks, including optical flow estimation. The
first example of supervised flow estimation using deep learning was FlowNet, which uses
a Unet-like architecture [135]; it is trained using labelled datasets, usually built from
synthetic video sequences; in this line of research, the current state-of-the-art method is
Recurrent All-Pairs Field Transforms (RAFT) [136]. Several unsupervised methods were
also published [137]–[139]. Those methods still use the basic components introduced by
[128], but they use convolutional neural networks (or transformers) as non-linear feature
learners. Furthermore, they consider specific training methods that leverage the advan-
tages of deep architectures. Current deep-learning-based unsupervised optical flow esti-
mation provides strong results [138], [140], even in cases where large motions or occlusions
are present.

The remainder of this chapter is organized as follows. In Section A.2, we present
the basics of unsupervised optical flow estimation with convolutional networks, with a
special focus on regularized training losses. In Section A.3, we present the dataset, as
well as a detailed reasoning for our use of segmentation as a proxy measurement when ground truth optical flow is unavailable in the context of microfluidics imaging. Finally, we show qualitative and quantitative results for the competing models also introduced in this section.

## A.2 Unsupervised CNN-based method for sparse optical flow estimation

### A.2.1 Definition of data term and loss function

For a pair of consecutive grayscale images $I_1, I_2 \in [0, 1]^{H \times W}$, let us consider a vector field $f_{12}$, where at each location $(u, v) \in \Omega = \{0, \dots, H\} \times \{0, \dots, W\}$ we have $f_{12}(u, v) = (\Delta_u, \Delta_v)$ so that $I_2(u+\Delta_u, v+\Delta_v) = I_1(u, v)$. This vector field quantifies the displacement of pixels. Note that the displacement is assumed to be small to derive linear equations with respect to image displacement. This is usually correct if the timesteps are small enough, i.e the framerate is high enough.

If the flow $f_{12}$ was perfectly estimated, we would be able to reconstruct $I_2$ from $I_1$ by resampling $I_1$ at the new locations given by the flow field. This resampling step is defined as a warping function:

$$w(I_1, f_{12}) = I_2. \qquad (A.1)$$

The reconstruction of $I_2$ from $I_1$ and the estimated flow field $f_{12}$ provides a natural criterion to evaluate the quality of the estimated flow field, usually called the photometric consistency criterion:

$$\mathcal{D}(I_1, I_2, f_{12}) = \sum_{u,v \in \Omega} \rho(I_2(u, v), w(I_1, f_{12})), \qquad (A.2)$$

where $\rho$ is a robust loss criterion (e.g. $\rho(x) = \sqrt{x^2 + \varepsilon^2}$). In general, variational methods minimize a cost function directly with respect to $f_{12}$ [127]. In deep-learning-based approaches, the major difference is that the training loss is minimized with respect to the estimator parameters. In our case, $f_{12} = g_{\theta}(I_1, I_2)$ where $g_{\theta}$, which is a deep convolutional neural network parameterized by a set of weights $\theta$. Given a dataset of image pairs $\{(I_1, I_2)_1, \dots, (I_1, I_2)_n\}$, the model parameters are obtained by minimizing the following
<figure><img src="image_71.png" /></figure>

<figure><img src="image_72.png" /><figcaption>Figure A.2 – Example of smearing in unsupervised optical flow estimation. Left: input grayscale image, Right: flow field estimated by UFlow without any regularization. Bottom left: HSV wheel corner indicating how the flow map (bottom) should be interpreted. Color indicates flow orientation, while saturation indicates flow norm.</figcaption></figure>

criterion as follows:

$$ \theta^* = \arg \min_{\theta} \sum_{i=1}^{n} \sum_{u,v} \rho(I_2(u, v), w(I_1, g_{\theta}(I_1, I_2)(u, v))). \quad (A.3) $$

Because of the non-linearities of $g_{\theta}$, this training loss is not convex, and finding a global minimum is not feasible. Accordingly, $\theta^*$ usually corresponds to a local minima of the loss function.

### A.2.2 Definition of spatial regularizers

Because of occlusions near object boundaries, matching a pixel between $I_1$ and $I_2$ may be very difficult. A single pixel in $I_1$ might have several low-error (in terms of intensity) matches in $I_2$. Because of this, unsupervised flow estimates often exhibit "smear", that is the flow boundaries do not match object boundaries, as shown in Figure A.2. As the data fitting may not be informative, regularization is used to impose some prior. The following 1st-order smoothness is commonly used in optical flow to adapt to object boundaries:

$$ \mathcal{R}(f_{12}, I_1) = \beta \left[ \sum_{u,v} \| \nabla_{u,v} f_{12} \|_1 e^{-\alpha |\nabla_{u,v} I_1(u,v)|} \right], \quad (A.4) $$
where $\alpha$ modulates edge weighting, and

$$ \nabla_{u,v} f_{12} = \begin{bmatrix} f_{12}(u, v) - f_{12}(u + 1, v) \\ f_{12}(u, v) - f_{12}(u, v + 1) \end{bmatrix}. \qquad (A.5) $$

Note that this is a simple linear transformation of the image so that the computational overhead is minimal. While this criterion forces the alignment of image and flow boundaries, it is unsatisfying because while object edges should also be flow boundaries, not all edges are object edges. Objects with internal edges can exhibit over-regularized flow with this regularizer.

In the context of microfluidics imaging, we know a priori that objects are small and exhibit relatively small motion. The background, which is most of the surface of the image, is expected to exhibit zero flow. Therefore, the flow field is expected to be sparse, that is the number of locations where the flow is non-zero is very low. Therefore, we can consider the following sparse variation regularizer, introduced in [11]:

$$ \mathcal{R}_{SV}(f_{12}) = \lambda \left[ \sum_{(u,v)} \rho^2 \| \nabla_{u,v} f_{12} \|_2^2 + (1 - \rho)^2 f_{12}^2(u, v) \right]. \qquad (A.6) $$

The first term encourages a smooth flow field, while the $f_{12}^2$ term controls the sparsity of the flow. The two terms are balanced with the hyperparameter $\rho$, usually $\rho \approx 0.1$. If $\rho = 1$, (A.6) is nothing else than the popular Tikhonov regularizer. In [11], the square root of this regularization term was taken to ensure convexity, which is not relevant in our training framework, as the non-linearities of $g_\theta$ makes the training objective non-convex.

In our implementation of Sparse Variation, we found that it is necessary to regularize the flow estimate at all layers of the feature pyramid of UFlow (see Section A.3.3 for details of UFlow). However, as the criterion in (A.3) is a sum over all pixels, the hyperparameter $\lambda$ in (4.3) needs to be adjusted at each scale to ensure it scales with the magnitude of the loss. As the height and width are divided by a factor 2 at each layer, we multiply $\lambda$ by a factor 4, and we optimize the average regularization over all layers:

$$ \mathcal{R}_{SVA}(f_{12}) = \frac{1}{L} \sum_{\ell=1}^{L} 4^{\ell} \lambda \mathcal{R}_{SV}(f_{12}^{\ell}), \qquad (A.7) $$

where $f_{12}^\ell$ the estimated flow at layer $\ell$.
<figure><img src="image_73.png" /></figure>

<figure><img src="image_74.png" /><figcaption>Figure A.3 – Single frame of the two sequence dataset used in this chapter. Left: crop of a frame from the red blood cell (RBC) sequence. A red blood cell is contained in the red square. Right: crop of a frame from the polystyrene and *yeast* sequence. A polystyrene bead is localized in the white square and a yeast in the blue square.</figcaption></figure>

# A.3 Experimental results and comparison of supervised and unsupervised methods

## A.3.1 Description of datasets

Our experiments were performed on several sequences acquired through shadow imaging in a microfluidics channel. In the two sequences used for illustration in this Section, each frame is of size $W = 1920$, $H = 1080$. In the first image sequence (*yeast* sequence) composed of 586 frames, the fluid carries polystyrene beads and yeast through the microfluidics channel. In the second image sequence (RBC sequence) comprised of 400 images, a blood sample is imaged, also with shadow imaging. The density of cells is much higher in the RBC sequence, which makes estimating flow boundaries accurately challenging. In Figure A.3, crops from each sequence are shown. As the camera is not moving and the objects are relatively small (around 10 by 10 pixels in a FullHD image), we estimated the background (for background subtraction) using a simple average of the sequence. In what follows, we only consider the background subtracted sequences.

## A.3.2 Evaluation metrics for optical flow estimation

In optical flow benchmark datasets, the ground truths are generally available [127]. If this is not the case, flow fields can be evaluated qualitatively. However, this might not be an objective way to compare several optical flow estimation methods. To provide
a fair comparison, we built segmentation annotation for the *yeast* and RBC sequences with Ilastik interactive segmentation tool [141]. Those annotations are used to build a proxy quality measurement for optical flow, focusing on the difference between flow and object boundaries. This is achieved by building a segmentation map of our optical flow estimation, and using conventional image segmentation metrics like Intersection Over Union (IoU).

To separate cells from the background using an optical flow estimate (i.e., to build a segmentation map), we chose to threshold the $L_2$ norm of the flow field. We translated the histogram of the norm to the [0, 1] interval to ensure that thresholds are comparable from frame to frame. Given that the speed of the cells varies during the sequence (especially for the RBC sequence), we estimated a new threshold for each frame using the Otsu method [80].

Moreover, the speed variation of cells inside a microfluidics channel may complicate the evaluation procedure. While background subtraction ensures that there is no truly stationary cell in the dataset, very slow cells might be estimated as background. Therefore, some segmentation errors are due to slow cells which are missing completely (i.e., every pixel of the cell is predicted to be a static, background pixel), and others are caused by imperfect flow discontinuities. Because of this, we split the segmentation error into two parts and compute two IoUs. For each cell in the true segmentation map, we check if there is a non-zero intersection with the segmented flow map. If this is not the case, the pixels of this cell are labelled as background. Here, $IoU_{boundary}$ is computed between this new true segmentation map and the flow segmentation map. $IoU_{boundary}$ corresponds exclusively to cell boundary estimation error. $IoU_{all}$ is the IoU between the original segmentation map and the segmented flow map. Therefore, $IoU_{boundary} > IoU_{all}$. This decomposition allows one to verify if a model reaches a better IoU because it detects slow cells better, or because it provides a better estimation of cell boundaries.

Finally, we also evaluate the model by comparing cell counts, which gives an indication of how well the model is able to separate close cells. If adjacent cells are insufficiently separated, the global cell count is underestimated. On the other hand, small segmentation artifacts lead to an overestimated cell count.

### A.3.3 Description of competing methods

In our experiments, We compared two different architectures. First, we choose RAFT [136] as the baseline method for supervised optical flow estimation, because of its high
<figure><img src="image_75.png" /><figcaption>Figure A.4 – Comparison of all flow models evaluated in the chapter. First column is a random sample of cells in the sequence, retrieved at the same timestep (first image of the RBC sequence). Second column is the unregularized UFlow model. Third column is UFlow regularized with first-order smoothness. Fourth column is RAFT, and last column is UFlow regularized with Sparse Variation.</figcaption></figure>

performance and availability in the `torchvision` package. We did not modify RAFT, and used the pre-trained weights available in `torchvision`. Second, we considered UFlow [137] as it is a performant and flexible unsupervised optical flow estimation architecture. UFlow [137] consists in a feature pyramid, where both images are repeatedly downsampled, and progressively deeper features are built using convolutional layers. At each layer, a cost volume is computed [135], and an optical flow field is estimated. In [137], the final flow field is estimated at 1/4 of the input resolution and then upsampled to compute the photometric loss. Nevertheless, we made some modifications compared to the original UFlow implementation to better fit our case study. Because dense mention is crucial to accurately estimate cell rotation, and therefore to identify cells, optical flow is estimated at the highest possible resolution, instead of 1/4 of the input resolution. We did not use bidirectional losses, occlusion estimation or self-supervision to keep training as simple as possible. In the end, three models are trained: UFlow without any regularization, UFlow regularized with a Sparse Variation criterion ($\rho = 0.1$ and $\lambda = 2 \times 10^{-10}$) and UFlow regularized with an edge-alignment criterion ($\beta = 0.1$ and $\alpha = 100$). We ran training for 25 epochs on the RBC sequence, using a batch size of 2, and a learning rate of 0.0001.

### A.3.4 Experimental results

**Qualitative evaluation** Figure A.4 shows a qualitative comparison of the different flow estimation methods. In the flow field estimated by RAFT, some cells that exhibit positive
<figure><img src="image_76.png" /><figcaption>Figure A.5 – Cell merging and fading issues in the flow field estimated by RAFT on the RBC sequence. Four consecutive crops (frames 387, 390, 393 and 396) are displayed. Although two cells are present, the boundaries estimation are only correct if the two cells are far enough.</figcaption></figure>

<figure><img src="image_77.png" /><figcaption>Figure A.6 – Crops of frames 278 (top) and 279 (bottom) of the flow field estimated by RAFT on the RBC sequence.</figcaption></figure>

flow in frame $t$ will exhibit zero flow in frame $t + 1$. The cell is effectively "switched off". While a stationary cell should exhibit 0 flow even though the pixel intensities are non-zero, the instantaneous switch background and foreground for specific cells is not realistic. Slowing cells should exhibit a progressively decreasing average flow norm. This also prevents long-term tracking of cells, which is especially problematic because long-term rotation patterns are a crucial feature for cell classification, as shown in [123]. Furthermore, Figure A.5 shows that RAFT sometimes merges cells that are close to each other. This prevents accurate cell counting, and makes cell identification harder. In Figure A.4, we see that 1st-order smoothness tends to produce flow fields with well separated cells. Nevertheless, it's worth noting that this regularizer generates unrealistic square cell borders.

Unsupervised models trained directly on the RBC sequence mitigate the issue of disappearing cells. However, the flow and cell boundary do not match as accurately as with RAFT. Regularization seems to improve this shortcoming.
<figure><img src="image_78.png" /><figcaption>Figure A.7 – IoU as a function of time. Top row is IoUs over time for RBC sequence. Bottom row is IoUs over time for yeast sequence. Left column is IoU for all cells, right column is IoU for detected cells only, as described in Section A.3.2. Four models are compared: RAFT (red), UFlow without any regularization (blue), UFlow with Sparse Variation regularization (green) and UFlow with edge-alignment regularization (black). Shaded areas indicate the ($q_5, q_{95}$) interval of counting performance over the last 5 epochs of training, for all models except RAFT, where only the last epoch weights were available in `torchvision`.</figcaption></figure>

**Quantitative evaluation of boundary estimation**

In Figure A.7, we see the results for the *yeast* sequence on the bottom row. The Sparse Variation model achieves the highest IoU for all cells. Interestingly, if we only consider *detected* cells, the performance gap with the unregularized model is very low which is also true for the RBC sequence, as shown in the top row of Figure A.7. This suggests most of the variation in terms of IoU across time is due to disappearing cells.

The variation of IoU across time is much higher for the RBC sequence than for the
*yeast* sequence, like for counting performance in Figure A.8. As explained before, this performance can be decomposed into an IoU for detected cells, and an IoU for all cells. Thanks to this decomposition, we can confirm that disappearing cells are the largest reason for this drop in IoU across time, as the IoU for detected cells is consistently much higher across time.

<figure><img src="image_79.png" /><figcaption>Figure A.8 – Counting error as a function of time. The error for both sequences is represented on the same figure, but both sequences are not the same length. The RBC sequence (red, black, blue and green lines) is shorter than the *yeast* sequence. Four models are compared: RAFT (red), UFlow without any regularization (blue), UFlow with Sparse Variation regularization (green) and UFlow with edge-alignment regularization (black). Shaded areas indicate the (*q*<sub>5</sub>, *q*<sub>95</sub>) interval of IoU over the last 5 epochs of training, for all models except RAFT, where only the last epoch weights were available in *torchvision*.</figcaption></figure>
**Quantitative evaluation of cell counting** As we see in Figure A.8, cell counting results are heavily speed-dependent, especially in the RBC sequence, although some architectures fare better than others. RAFT produces the highest overestimation of the number of cells. Worse, this overestimation is inconsistent, with very high peaks. Those peaks correspond to frames where RAFT produces high-error outputs, as displayed in Figure A.6. In those frames, the optical flow estimation is over-smoothed, with a very strong background motion. Furthermore, the cells are not visible anymore.

UFlow produces a much more accurate cell count, especially for two models: the unregularized model (in blue in Figure A.7), and the model regularized with Sparse Variation (in green). While the unregularized model produces slightly better results for the RBC sequence, the variation in IoU across training epochs (materialized by the shaded confidence intervals) is much higher.

In the *yeast* sequence, we see that the counting error shows less variation across time, because of the constant speed of the cells in the channel. The Sparse Variation model achieves the lowest counting error, and the lowest variability across training epochs.

## A.4 Conclusion

This chapter shows the limitations of pre-trained, supervised optical flow models for analyzing the motion of small objects in a microfluidics setting. Because of the distribution gap between the training and testing dataset, the estimation of an accurate flow field for slow cells is unsatisfying. Cells tend to disappear unexpectedly, even though cells are not stationary.

Unsupervised optical flow estimation addresses this issue, but the discontinuities of the flow do not match the object boundaries. While edge-alignment regularization is used very often in optical flow to address this issue, it performs worse than our novel Sparse Variation regularization. The model trained with Sparse Variation regularization achieves the best performance (sometimes tied with the unregularized model), and the lowest variation in terms of IoU and cell counting performance over training epochs.
<figure><img src="image_80.png" /></figure>
References

[1] A. Krizhevsky, I. Sutskever, and G. E. Hinton, « ImageNet Classification with Deep Convolutional Neural Networks », in *Advances in Neural Information Processing Systems (NeurIPS)*, vol. 25, 2012.

[2] O. Russakovsky, J. Deng, H. Su, et al., « ImageNet Large Scale Visual Recognition Challenge », *International Journal of Computer Vision*, vol. 115, pp. 211–252, 2015.

[3] T.-Y. Lin, M. Maire, S. Belongie, et al., « Microsoft COCO: Common Objects in Context », in *European Conference on Computer Vision (ECCV)*, 2014, pp. 740–755.

[4] M. Everingham, L. Van Gool, C. K. I. Williams, J. Winn, and A. Zisserman, « The Pascal Visual Object Classes (VOC) Challenge », *International Journal of Computer Vision*, vol. 88, pp. 303–338, 2010.

[5] C. Guo, G. Pleiss, Y. Sun, and K. Q. Weinberger, « On Calibration of Modern Neural Networks », in *International Conference on Machine Learning (ICML)*, 2017, pp. 1321–1330.

[6] W. J. Maddox, P. Izmailov, T. Garipov, D. P. Vetrov, and A. G. Wilson, « A simple baseline for bayesian uncertainty in deep learning », *Advances in Neural Information Processing Systems (NeurIPS)*, vol. 32, 2019.

[7] S. Mandt, M. Hoffman, and D. Blei, « Stochastic Gradient Descent as Approximate Bayesian Inference », *Journal of Machine Learning Research*, vol. 18, 2017.

[8] K. He, X. Zhang, S. Ren, and J. Sun, « Deep Residual Learning for Image Recognition », in *IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)*, 2016, pp. 770–778.

[9] S. Ren, K. He, R. Girshick, and J. Sun, « Faster R-CNN: Towards Real-Time Object Detection with Region Proposal Networks », in *Advances in Neural Information Processing Systems (NeurIPS)*, vol. 28, 2015.

[10] O. Ronneberger, P. Fischer, and T. Brox, « U-Net: Convolutional Networks for Biomedical Image Segmentation », en, in *Medical Image Computing and Computer-Assisted Intervention (MICCAI)*, 2015, pp. 234–241.
[11] S. Prigent, H.-N. Nguyen, L. Leconte, et al., « SPITFIR(e): a supermaneuverable algorithm for fast denoising and deconvolution of 3D fluorescence microscopy images and videos », *Scientific Reports*, vol. 13, p. 1489, 2023.

[12] J. B. Diederik Kingma, « Adam: A Method for Stochastic Optimization », in *International Conference on Learning Representations (ICLR)*, 2015.

[13] P. Izmailov, D. Podoprikhin, T. Garipov, D. Vetrov, and A. Wilson, « Averaging Weights Leads to Wider Optima and Better Generalization », 2018.

[14] I. A. E. Agency, « Cytogenetic Dosimetry: Applications in Preparedness for and Response to Radiation Emergencies », International Atomic Energy Agency, Text, 2011, p. 1.

[15] R. Balestriero, M. Ibrahim, V. Sobal, et al., *A Cookbook of Self-Supervised Learning*, arXiv:2304.12210 [cs], 2023.

[16] C. Cortes and V. Vapnik, « Support-vector networks », *Machine Learning*, vol. 20, pp. 273–297, 1995.

[17] N. Dalal and B. Triggs, « Histograms of Oriented Gradients for Human Detection », in *IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)*, 2005, pp. 886–893.

[18] D. G. Lowe, « Distinctive Image Features from Scale-Invariant Keypoints », *International Journal of Computer Vision*, vol. 60, pp. 91–110, 2004.

[19] M. Abadi, P. Barham, J. Chen, et al., « TensorFlow: a system for large-scale machine learning », in *USENIX conference on Operating Systems Design and Implementation*, 2016, pp. 265–283.

[20] A. Paszke, S. Gross, F. Massa, et al., « PyTorch: An Imperative Style, High-Performance Deep Learning Library », in *Advances in Neural Information Processing Systems (NeurIPS)*, vol. 32, 2019.

[21] G. Cybenko, « Approximation by superpositions of a sigmoidal function », *Mathematics of Control, Signals and Systems*, vol. 2, pp. 303–314, 1989.

[22] A. Botev, H. Ritter, and D. Barber, « Practical Gauss-Newton Optimisation for Deep Learning », in *International Conference on Machine Learning (ICML)*, 2017, pp. 557–565.
[23] H. Ritter, A. Botev, and D. Barber, « A Scalable Laplace Approximation for Neural Networks », in *International Conference on Learning Representations (ICLR)*, 2022.

[24] J. Martens and R. Grosse, « Optimizing neural networks with Kronecker-factored approximate curvature », in *International Conference on Machine Learning (ICML)*, 2015, pp. 2408–2417.

[25] H. Robbins and S. Monro, « A Stochastic Approximation Method », *The Annals of Mathematical Statistics*, vol. 22, pp. 400–407, 1951.

[26] S. L. Smith, B. Dherin, D. Barrett, and S. De, « On the Origin of Implicit Regularization in Stochastic Gradient Descent », in *International Conference on Learning Representations (ICLR)*, 2020.

[27] C. Zhang, S. Bengio, M. Hardt, B. Recht, and O. Vinyals, « Understanding deep learning requires rethinking generalization », in *International Conference on Learning Representations (ICLR)*, 2016.

[28] K. Simonyan, A. Vedaldi, and A. Zisserman, « Deep Inside Convolutional Networks: Visualising Image Classification Models and Saliency Maps », in *International Conference on Learning Representations (ICLR)*, Y. Bengio and Y. LeCun, Eds., 2014.

[29] L. A. Gatys, A. S. Ecker, and M. Bethge, « Image Style Transfer Using Convolutional Neural Networks », in *IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)*, 2016, pp. 2414–2423.

[30] D. Ulyanov, A. Vedaldi, and V. Lempitsky, « Deep Image Prior », *International Journal of Computer Vision*, vol. 128, pp. 1867–1888, 2020.

[31] A. Rosenfeld and J. K. Tsotsos, « Intriguing Properties of Randomly Weighted Networks: Generalizing While Learning Next to Nothing », in *Conference on Computer and Robot Vision (CRV)*, 2019, pp. 9–16.

[32] J. Frankle, D. J. Schwab, and A. S. Morcos, « Training BatchNorm and Only Batch-Norm: On the Expressive Power of Random Features in CNNs », in *International Conference on Learning Representations (ICLR)*, 2020.

[33] S. d'Ascoli, L. Sagun, G. Biroli, and J. Bruna, « Finding the Needle in the Haystack with Convolutions: on the benefits of architectural bias », in *Advances in Neural Information Processing Systems (NeurIPS)*, vol. 32, 2019.
[34] T. Hastie, R. Tibshirani, and J. H. Friedman, *The elements of statistical learning: data mining, inference, and prediction* (Springer series in statistics), 2nd ed. New York, 2009.

[35] I. Goodfellow, J. Pouget-Abadie, M. Mirza, et al., « Generative Adversarial Nets », in *Advances in Neural Information Processing Systems (NeurIPS)*, vol. 27, 2014.

[36] D. P. Kingma and M. Welling, *Auto-Encoding Variational Bayes*, arXiv:1312.6114 [stat.ML], 2022.

[37] T. Karras, M. Aittala, S. Laine, et al., « Alias-Free Generative Adversarial Networks », in *Advances in Neural Information Processing Systems (NeurIPS)*, vol. 34, 2021, pp. 852–863.

[38] M. Arjovsky and L. Bottou, « Towards Principled Methods for Training Generative Adversarial Networks », in *International Conference on Learning Representations (ICLR)*, 2016.

[39] A. Radford, L. Metz, and S. Chintala, « Unsupervised Representation Learning with Deep Convolutional Generative Adversarial Networks », in *International Conference on Learning Representations (ICLR)*, 2016.

[40] I. J. Goodfellow, J. Shlens, and C. Szegedy, « Explaining and Harnessing Adversarial Examples », in *International Conference on Learning Representations (ICLR)*, Y. Bengio and Y. LeCun, Eds., 2015.

[41] M. Mirza and S. Osindero, *Conditional Generative Adversarial Nets*, arXiv:1411.1784 [cs.LG], 2014.

[42] L. Mescheder, A. Geiger, and S. Nowozin, « Which Training Methods for GANs do actually Converge? », arXiv:1801.04406 [cs.LG], 2018.

[43] P. Isola, J.-Y. Zhu, T. Zhou, and A. A. Efros, « Image-to-Image Translation with Conditional Adversarial Networks », in *IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)*, 2017, pp. 5967–5976.

[44] P. Ghosh, M. S. M. Sajjadi, A. Vergari, M. Black, and B. Scholkopf, « From Variational to Deterministic Autoencoders », arXiv:1903.12436 [cs.LG], 2019.

[45] G. E. Hinton and R. Zemel, « Autoencoders, Minimum Description Length and Helmholtz Free Energy », in *Advances in Neural Information Processing Systems (NeurIPS)*, vol. 6, 1993.
[46] I. Higgins, L. Matthey, A. Pal, et al., « Beta-VAE: Learning Basic Visual Concepts with a Constrained Variational Framework », in *International Conference on Learning Representations (ICLR)*, 2016.

[47] T. Miyato, T. Kataoka, M. Koyama, and Y. Yoshida, « Spectral Normalization for Generative Adversarial Networks », in *International Conference on Learning Representations (ICLR)*, 2018.

[48] K. He, H. Fan, Y. Wu, S. Xie, and R. Girshick, « Momentum Contrast for Unsupervised Visual Representation Learning », in *IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)*, 2020, pp. 9726–9735.

[49] T. Chen, S. Kornblith, M. Norouzi, and G. Hinton, « A Simple Framework for Contrastive Learning of Visual Representations », in *International Conference on Machine Learning (ICML)*, 2020, pp. 1597–1607.

[50] Z. Wu, Y. Xiong, S. X. Yu, and D. Lin, « Unsupervised Feature Learning via Non-parametric Instance Discrimination », in *2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition*, Salt Lake City, UT, 2018, pp. 3733–3742.

[51] Y. LeCun, C. Cortes, and C. Burges, *MNIST handwritten digit database*. [Online]. Available: http://yann.lecun.com/exdb/mnist/.

[52] M. Caron, P. Bojanowski, A. Joulin, and M. Douze, « Deep Clustering for Un-supervised Learning of Visual Features », in *European Conference on Computer Vision (ECCV)*, V. Ferrari, M. Hebert, C. Sminchisescu, and Y. Weiss, Eds., 2018, pp. 139–156.

[53] M. Caron, I. Misra, J. Mairal, P. Goyal, P. Bojanowski, and A. Joulin, « Un-supervised Learning of Visual Features by Contrasting Cluster Assignments », in *Advances in Neural Information Processing Systems (NeurIPS)*, vol. 33, 2020, pp. 9912–9924.

[54] P. A. Knight, « The Sinkhorn–Knopp Algorithm: Convergence and Applications », *SIAM Journal on Matrix Analysis and Applications*, vol. 30, pp. 261–275, 2008.

[55] X. Chen and K. He, « Exploring Simple Siamese Representation Learning », in *IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)*, 2021, pp. 15750–15758.
[56] E. Daxberger, A. Kristiadi, A. Immer, R. Eschenhagen, M. Bauer, and P. Hennig, « Laplace Redux - Effortless Bayesian Deep Learning », in *Advances in Neural Information Processing Systems*, vol. 34, 2021, pp. 20089–20103.

[57] B. Lakshminarayanan, A. Pritzel, and C. Blundell, « Simple and Scalable Predictive Uncertainty Estimation using Deep Ensembles », in *Advances in Neural Information Processing Systems (NeurIPS)*, vol. 30, 2017.

[58] G. Huang, Y. Li, G. Pleiss, Z. Liu, J. E. Hopcroft, and K. Q. Weinberger, « Snapshot Ensembles: Train 1, Get M for Free », in *International Conference on Learning Representations (ICLR)*, 2016.

[59] T. Garipov, P. Izmailov, D. Podoprikhin, D. P. Vetrov, and A. G. Wilson, « Loss Surfaces, Mode Connectivity, and Fast Ensembling of DNNs », in *Advances in Neural Information Processing Systems (NeurIPS)*, vol. 31, 2018.

[60] D. Ruppert, « Efficient Estimations from a Slowly Convergent Robbins-Monro Process », Tech. Rep., 1988.

[61] B. Polyak, « New stochastic approximation type procedures », *Avtomatica i Telemekhanika*, vol. 7, pp. 98–107, 1990.

[62] S. Hochreiter and J. Schmidhuber, « Flat Minima », *Neural computation*, vol. 9, pp. 1–42, 1997.

[63] N. S. Keskar, D. Mudigere, J. Nocedal, M. Smelyanskiy, and P. T. P. Tang, « On Large-Batch Training for Deep Learning: Generalization Gap and Sharp Minima », in *International Conference on Learning Representations (ICLR)*, 2016.

[64] J. A. Hoeting, D. Madigan, A. E. Raftery, and C. T. Volinsky, « Bayesian model averaging: a tutorial », *Statistical Science*, vol. 14, 4, pp. 382–401, 1999.

[65] R. Girshick, « Fast R-CNN », in *IEEE International Conference on Computer Vision (ICCV)*, 2015, pp. 1440–1448.

[66] J. R. R. Uijlings, K. E. A. Van De Sande, T. Gevers, and A. W. M. Smeulders, « Selective Search for Object Recognition », *International Journal of Computer Vision*, vol. 104, pp. 154–171, 2013.

[67] X. Zhou, D. Wang, and P. Krähenbühl, *Objects as Points*, arXiv:1904.07850 [cs.CV], 2019.
[68] T.-Y. Lin, P. Goyal, R. Girshick, K. He, and P. Dollár, « Focal Loss for Dense Object Detection », *IEEE Transactions on Pattern Analysis and Machine Intelligence*, vol. 42, pp. 318-327, 2020.

[69] H. Law and J. Deng, « CornerNet: Detecting Objects as Paired Keypoints », *International Journal of Computer Vision*, vol. 128, pp. 642-656, 2020.

[70] R. Bayley, A. Carothers, X. Chen, et al., « Radiation dosimetry by automatic image analysis of dicentric chromosomes », *Mutation Research/Environmental Mutagenesis and Related Subjects*, vol. 253, pp. 223-235, 1991.

[71] L. Roy, M. Delbos, N. Paillole, V. Durand, and P. Voisin, « Comparaison de systèmes d’analyse d’images cytologiques en dosimétrie biologique », *Radioprotection*, vol. 38, pp. 323-340, 2003.

[72] L. Roy, I. Sorokine-Durm, and P. Voisin, « Comparison between fluorescence in situ hybridization and conventional cytogenetics for dicentric scoring: a first-step validation for the use of FISH in biological dosimetry », *International Journal of Radiation Biology*, vol. 70, pp. 665-669, 1996.

[73] R. Huber, U. Kulka, T. Lörch, et al., « Technical report: application of the Metafer2 fluorescence scanning system for the analysis of radiation-induced chromosome aberrations measured by FISH-chromosome painting », *Mutation Research/Genetic Toxicology and Environmental Mutagenesis*, vol. 492, pp. 51-57, 2001.

[74] B. C. Shirley, J. H. M. Knoll, J. Moquet, et al., « Estimating partial-body ionizing radiation exposure by automated cytogenetic biodosimetry », *International Journal of Radiation Biology*, vol. 96, 11, pp. 1492-1503, Nov. 2020.

[75] A. Subasinghe, J. Samarabandu, Y. Li, et al., « Centromere detection of human metaphase chromosome images using a candidate based method », *F1000Research*, 2016.

[76] P. Finnon, D. C. Lloyd, and A. A. Edwards, « An assessment of the metaphase finding capability of the Cytoscan 110 », *Mutation Research*, vol. 164, pp. 101-108, 1986.

[77] J. Piper and J. Sprey, « Adaptive classifiers for dicentric chromosomes », *Journal of Radiation Research*, vol. 33, pp. 159-170, 1992.

[78] P. Finnon and D. Lloyd, « A preliminary evaluation of the Edinburgh dicentric hunter », *Journal of Radiation Research*, vol. 33 Suppl, pp. 215-221, 1992.
[79] Y. Li, A. Wickramasinghe, A. A. Subasinghe, et al., « Towards large scale auto-
mated interpretation of cytogenetic biodosimetry data », in *IEEE International
Conference on Information and Automation for Sustainability*, 2012, pp. 30–35.

[80] N. Otsu, « A Threshold Selection Method from Gray-Level Histograms », *IEEE Transactions on Systems, Man, and Cybernetics*, vol. 9, pp. 62–66, 1979.

[81] A. Subasinghe Arachchige, J. Samarabandu, J. Knoll, W. Khan, and P. Rogan, « An image processing algorithm for accurate extraction of the centerline from human metaphase chromosomes », in *IEEE International Conference on Image Processing*, 2010, pp. 3613–3616.

[82] A. S. Arachchige, J. Samarabandu, J. H. M. Knoll, and P. K. Rogan, « Intensity Integrated Laplacian-Based Thickness Measurement for Detecting Human Metaphase Chromosome Centromere Location », *IEEE Transactions on Biomedical Engineering*, vol. 60, pp. 2005–2013, 2013.

[83] X. Bai, L. Latecki, and W.-y. Liu, « Skeleton Pruning by Contour Partitioning with Discrete Curve Evolution », *IEEE Transactions on Pattern Analysis and Machine Intelligence*, vol. 29, pp. 449–462, 2007.

[84] J. Liu, Y. Li, R. Wilkins, F. Flegal, J. H. Knoll, and P. K. Rogan, « Accurate cytogenetic biodosimetry through automated dicentric chromosome curation and metaphase cell selection », *F1000Research*, vol. 6, p. 1396, 2017.

[85] Y. LI, J. Knoll, R. Wilkins, F. Norton, and P. Rogan, « Automated discrimination of dicentric and monocentric chromosomes by machine learning-based image processing », *Microscopy Research and Technique*, vol. 79, 2016.

[86] P. K. Rogan, Y. Li, A. Wickramasinghe, et al., « Automating dicentric chromosome detection from cytogenetic biodosimetry data », *Radiation Protection Dosimetry*, vol. 159, pp. 95–104, 2014.

[87] A. S. A., J. Samarabandu, J. Knoll, W. Khan, and P. Rogan, « An Accurate Image Processing Algorithm for Detecting FISH Probe Locations Relative to Chromosome Landmarks on DAPI Stained Metaphase Chromosome Images », in *Canadian Conference on Computer and Robot Vision*, 2010, pp. 223–230.

[88] P. K. Rogan, Y. Li, R. C. Wilkins, F. N. Flegal, and J. H. M. Knoll, « Radiation Dose Estimation by Automated Cytogenetic Biodosimetry », *Radiation Protection Dosimetry*, vol. 172, pp. 207–217, 2016.
[89] Y. Li, B. C. Shirley, R. C. Wilkins, F. Norton, J. H. M. Knoll, and P. K. Rogan, « Radiation Dose Estimation By Completely Automated Interpretation Of The Dicentric Chromosome Assay », *Radiation Protection Dosimetry*, 2019.

[90] S. Jang, S.-G. Shin, M.-J. Lee, et al., « Feasibility Study on Automatic Interpretation of Radiation Dose Using Deep Learning Technique for Dicentric Chromosome Assay », *Radiation Research*, vol. 195, pp. 163-172, 2021.

[91] X. Shen, T. Ma, C. Li, Z. Wen, J. Zheng, and Z. Zhou, « High-precision automatic identification method for dicentric chromosome images using two-stage convolutional neural network », *Scientific Reports*, vol. 13, p. 2124, Feb. 2023.

[92] P. Voisin, M. Benderitter, M. Claraz, et al., « The cytogenetic dosimetry of recent accidental overexposure », *Cellular and Molecular Biology*, vol. 47, pp. 557-564, 2001.

[93] A. Vaurijoux, G. Gruel, F. Pouzoulet, et al., « Strategy for population triage based on dicentric analysis », *Radiation Research*, vol. 171, pp. 541-548, 2009.

[94] A. Vaurijoux, E. Gregoire, S. Roch-Lefevre, et al., « Detection of Partial-Body Exposure to Ionizing Radiation by the Automatic Detection of Dicentrics », *Radiation Research*, vol. 178, pp. 357-364, 2012.

[95] G. Gruel, E. Grégoire, S. Lecas, et al., « Biological dosimetry by automated dicentric scoring in a simulated emergency », *Radiation research*, vol. 179, 2013.

[96] H. Romm, L. Ainsbury, S. Barnard, et al., « Automatic scoring of dicentric chromosomes as a tool in large scale radiation accidents », *Mutation research*, vol. 756, May 2013.

[97] T. L. Ryan, M. B. Escalona, T. L. Smith, J. Albanese, C. J. Iddins, and A. S. Balajee, « Optimization and validation of automated dicentric chromosome analysis for radiological/nuclear triage applications », *Mutation Research/Genetic Toxicology and Environmental Mutagenesis*, vol. 847, p. 503-87, Nov. 2019.

[98] D. Endesfelder, U. Kulka, J. Einbeck, and U. Oestreicher, « Improving the accuracy of dose estimates from automatically scored dicentric chromosomes by accounting for chromosome number », *International Journal of Radiation Biology*, vol. 96, 12, pp. 1571-1584, Dec. 2020.
[99] A. Defazio, F. Bach, and S. Lacoste-Julien, « SAGA: a fast incremental gradient method with support for non-strongly convex composite objectives », in *Advances in Neural Information Processing Systems (NeurIPS)*, 2014, pp. 1646–1654.

[100] F. Pedregosa, G. Varoquaux, A. Gramfort, et al., « Scikit-learn: Machine Learning in Python », *Journal of Machine Learning Research*, vol. 12, pp. 2825–2830, 2011.

[101] R. Padilla, S. L. Netto, and E. A. B. da Silva, « A Survey on Performance Metrics for Object-Detection Algorithms », in *International Conference on Systems, Signals and Image Processing (IWSSIP)*, 2020, pp. 237–242.

[102] A. Bulat and G. Tzimiropoulos, « Binarized Convolutional Landmark Localizers for Human Pose Estimation and Face Alignment with Limited Resources », in *IEEE International Conference on Computer Vision (ICCV)*, 2017, pp. 3726–3734.

[103] A. Bulat and G. Tzimiropoulos, « How far are we from solving the 2D and 3D Face Alignment problem? (and a dataset of 230,000 3D facial landmarks) », in *IEEE International Conference on Computer Vision (ICCV)*, 2017, pp. 1021–1030.

[104] V. Lempitsky and A. Zisserman, « Learning To Count Objects in Images », in *Advances in Neural Information Processing Systems (NeurIPS)*, vol. 23, 2010.

[105] J. P. Cohen, G. Boucher, C. A. Glastonbury, H. Z. Lo, and Y. Bengio, « Count-ception: Counting by Fully Convolutional Redundant Counting », 2017, pp. 18–26.

[106] P. Chattopadhyay, R. Vedantam, R. R. Selvaraju, D. Batra, and D. Parikh, « Counting Everyday Objects in Everyday Scenes », in *IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)*, 2017, pp. 4428–4437.

[107] I. H. Laradji, N. Rostamzadeh, P. O. Pinheiro, D. Vazquez, and M. Schmidt, « Where Are the Blobs: Counting by Localization with Point Supervision », in *European Conference on Computer Vision (ECCV)*, vol. 11206, 2018, pp. 560–576.

[108] F. Wang, G. Wei, Q. Liu, J. Ou, x. wei xian, and H. Lv, « Boost Neural Networks by Checkpoints », in *Advances in Neural Information Processing Systems (NeurIPS)*, vol. 34, 2021, pp. 19719–19729.

[109] R. E. Schapire, « A brief introduction to boosting », in *IJCAI International Joint Conference on Artificial Intelligence*, vol. 2, 1999, pp. 1401–1406.
[110] D. Marmanis, J. D. Wegner, S. Galliani, K. Schindler, M. Datcu, and U. Stilla, « Semantic Segmentation Of Aerial Images With An Ensemble Of CNNs », *ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences*, vol. 3, pp. 473-480, 2016.

[111] K. Shaga Devan, H. A. Kestler, C. Read, and P. Walther, « Weighted average ensemble-based semantic segmentation in biological electron microscopy images », *Histochemistry and Cell Biology*, vol. 158, pp. 447-462, 2022.

[112] L. Nanni, C. Fantozzi, A. Loreggia, and A. Lumini, « Ensembles of Convolutional Neural Networks and Transformers for Polyp Segmentation », *Sensors*, vol. 23, p. 4688, 2023.

[113] A. Neubeck and L. Van Gool, « Efficient Non-Maximum Suppression », in *International Conference on Pattern Recognition (ICPR)*, 2006, pp. 850–855.

[114] R. Solovyev, W. Wang, and T. Gabruseva, « Weighted boxes fusion: Ensembling boxes from different object detection models », *Image and Vision Computing*, vol. 107, p. 104 117, 2021.

[115] A. G. Wilson and P. Izmailov, « Bayesian Deep Learning and a Probabilistic Perspective of Generalization », in *Advances in Neural Information Processing Systems (NeurIPS)*, vol. 33, 2020, pp. 4697–4708.

[116] G. Franchi, A. Bursuc, E. Aldea, S. Dubuisson, and I. Bloch, « TRADI: Tracking Deep Neural Network Weight Distributions », in *European Conference on Computer Vision (ECCV)*, Series Title: Lecture Notes in Computer Science, vol. 12362, 2020, pp. 105–121. DOI: 10.1007/978-3-030-58520-4_7.

[117] M. Li and C. Scheidegger, *Comparing Deep Neural Nets with UMAP Tour*, 2021.

[118] L. McInnes, J. Healy, N. Saul, and L. Großberger, « UMAP: Uniform Manifold Approximation and Projection », *Journal of Open Source Software*, vol. 3, 29, p. 861, 2018.

[119] J. C. Gower, « Generalized procrustes analysis », *Psychometrika*, vol. 40, pp. 33–51, 1975.

[120] A. Araujo, W. Norris, and J. Sim, « Computing Receptive Fields of Convolutional Neural Networks », *Distill*, 2019.
[121] M. Pachitariu and C. Stringer, « Cellpose 2.0: how to train your own model », *Nature Methods*, vol. 19, 12, pp. 1634-1641, Dec. 2022, Number: 12 Publisher: Nature Publishing Group.

[122] W. Zhang, J. Pang, K. Chen, and C. C. Loy, « Dense Siamese Network for Dense Unsupervised Learning », in *European Conference on Computer Vision (ECCV)*, 2022, pp. 464-480.

[123] J. Kun, M. Smieja, B. Xiong, L. Soleymani, and Q. Fang, « The Use of Motion Analysis as Particle Biomarkers in Lensless Optofluidic Projection Imaging for Point of Care Urine Analysis », *Scientific Reports*, vol. 9, p. 17255, 2019.

[124] B. K. P. Horn and B. G. Schunck, « Determining optical flow », *Artificial Intelligence*, vol. 17, pp. 185-203, 1981.

[125] T. Brox, A. Bruhn, N. Papenberg, and J. Weickert, « High Accuracy Optical Flow Estimation Based on a Theory for Warping », in *European Conference on Computer Vision (ECCV)*, T. Pajdla and J. Matas, Eds., 2004, pp. 25-36.

[126] E. Mémin and P. Pérez, « Hierarchical Estimation and Segmentation of Dense Motion Fields », *International Journal of Computer Vision*, vol. 46, pp. 129-155, 2002.

[127] D. Fortun, P. Bouthemy, and C. Kervrann, « Optical Flow Modeling and Computation: A Survey », *Computer Vision and Image Understanding*, vol. 134, 2015.

[128] B. D. Lucas and T. Kanade, « An iterative image registration technique with an application to stereo vision », in *International Joint Conference on Artificial Intelligence (IJCAI)*, vol. 2, 1981, pp. 674-679.

[129] E. Mémin and P. Pérez, « Dense estimation and object-based segmentation of the optical flow with robust techniques », *IEEE Transactions on Image Processing*, vol. 7, pp. 703-719, 1998.

[130] J. Delpiano, J. Jara, J. Scheer, O. A. Ramírez, J. Ruiz-del-Solar, and S. Härtel, « Performance of optical flow techniques for motion analysis of fluorescent point signals in confocal microscopy », *Machine Vision and Applications*, vol. 23, pp. 675-689, 2012.
[131] D. Fortun, P. Bouthemy, P. Paul-Gilloteaux, and C. Kervrann, « Aggregation of patch-based estimations for illumination-invariant optical flow in live cell imaging », in *IEEE International Symposium on Biomedical Imaging*, 2013, pp. 660–663.

[132] D. Fortun, N. Debroux, and C. Kervrann, « Spatially-Variant Kernel for Optical Flow Under Low Signal-to-Noise Ratios Application to Microscopy », in *IEEE International Conference on Computer Vision Workshops (ICCVW)*, 2017, pp. 42–48.

[133] S. Manandhar, P. Bouthemy, E. Welf, G. Danuser, P. Roudot, and C. Kervrann, « 3D flow field estimation and assessment for live cell fluorescence microscopy », *Bioinformatics*, vol. 36, 5, pp. 1317–1325, 2020.

[134] D. Fortun, P. Bouthemy, and C. Kervrann, « Aggregation of local parametric candidates with exemplar-based occlusion handling for optical flow », *Computer Vision and Image Understanding*, Light Field for Computer Vision, vol. 145, pp. 81–94, 2016.

[135] A. Dosovitskiy, P. Fischer, E. Ilg, et al., « FlowNet: Learning Optical Flow with Convolutional Networks », in *IEEE International Conference on Computer Vision (ICCV)*, 2015, pp. 2758–2766.

[136] Z. Teed and J. Deng, « RAFT: Recurrent All-Pairs Field Transforms for Optical Flow », in *European Conference on Computer Vision (ECCV)*, vol. 12347, 2020, pp. 402–419.

[137] R. Jonschkowski, A. Stone, J. T. Barron, A. Gordon, K. Konolige, and A. Angelova, « What Matters in Unsupervised Optical Flow », in *European Conference on Computer Vision (ECCV)*, 2020, pp. 557–572.

[138] L. Liu, J. Zhang, R. He, et al., « Learning by Analogy: Reliable Supervision From Transformations for Unsupervised Optical Flow Estimation », in *IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)*, 2020, pp. 6488–6497.

[139] Z. Huang, X. Shi, C. Zhang, et al., « FlowFormer: A Transformer Architecture for Optical Flow », in *European Conference on Computer Vision (ECCV)*, vol. 13677, 2022, pp. 668–685.
[140] A. Stone, D. Maurer, A. Ayvaci, A. Angelova, and R. Jonschkowski, « SMURF: Self-Teaching Multi-Frame Unsupervised RAFT with Full-Image Warping », in *IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)*, 2021, pp. 3886–3895.

[141] S. Berg, D. Kutra, T. Kroeger, et al., « Ilastik: interactive machine learning for (bio)image analysis », *Nature Methods*, vol. 16, pp. 1226–1232, 2019.
<figure><img src="image_81.png" /></figure>
**Titre :** Apprentissage machine et réseaux de convolution pour une expertise augmentée en dosimétrie biologique

**Mot clés :** Apprentissage profond, Agrégation de modèles, Incertitude, Détection d'objets, Imagerie médicale

**Résumé :** La dosimétrie biologique est la branche de la physique de la santé qui se préoccupe de l'estimation de doses de rayonnement ionisants à partir de biomarqueurs. Dans le procédé standard défini par l'AIEA, la dose est calculée en estimant la fréquence d'apparition de chromosomes dicentriques lors de la métaphase des lymphocytes périphériques. La variabilité morphologique des chromosomes, ainsi que celle des conditions d'acquisition des images rend ce problème de détection d'objets complexe. De plus, l'estimation fiable de cette fréquence nécessite le traitement d'un grand nombre d'image. Etant donné les limites du comptage

humain (faible nombre de personnes qualifiées, charge cognitive), l'automatisation est une nécessité dans le contexte d'exposition de masse. Dans ce contexte, l'objectif de cette thèse est de tirer parti des progrès récents en vision par ordinateur (et plus spécifiquement en détection d'objets) apportés par l'apprentissage profond. La contribution principale de ce travail est une preuve de concept pour un modèle de détection de chromosomes dicentriques. Ce système repose sur l'agrégation de modèles pour parvenir à un haut niveau de performance, ainsi qu'a une bonne quantification de son incertitude, une exigence essentielle dans un contexte médical.

**Title:** Machine learning and convolutional networks for automated biological dosimetry

**Keywords:** Deep learning, Model Aggregation, Uncertainty, Object detection, Medical imaging

**Abstract:** Biological dosimetry is the branch of health physics dealing with the estimation of ionizing radiation doses from biomarkers. The current gold standard (defined by the IAEA) relies on estimating how frequently dicentric chromosomes appear in peripheral blood lymphocytes. Variations in acquisition conditions and chromosome morphology makes this a challenging object detection problem. Furthermore, the need for an accurate estimation of the average number of dicentric per cell means that a large number of image has to be processed. Human counting is intrinsically

limited, as cognitive load is high and the number of specialist insufficient in the context of a large-scale exposition. The main goal of this PhD is to use recent developments in computer vision brought by deep learning, especially for object detection. The main contribution of this thesis is a proof of concept for a dicentric chromosome detection model. This model aggregates several Unet models to reach a high level of performance and quantify its prediction uncertainty, which is a stringent requirement in a medical setting.