# Machine learning and convolutional networks for automated biological dosimetry

Antonin Deschemps

► To cite this version:

Antonin Deschemps. Machine learning and convolutional networks for automated biological dosimetry. Machine Learning [cs.LG]. Université de Rennes, 2023. English. NNT: 2023URENS104. tel-04597585

HAL Id: tel-04597585

https://theses.hal.science/tel-04597585

Submitted on 3 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
## THÈSE DE DOCTORAT DE

L'UNIVERSITÉ DE RENNES

ÉCOLE DOCTORALE N° 601
Mathématiques, Télécommunications, Informatique, Signal, Systèmes, Électronique
Spécialité : Signal, Image, Vision

Par
Antonin DESCHEMPS

## Apprentissage machine et réseaux de convolution pour une expertise augmentée en dosimétrie biologique

Thèse présentée et soutenue à Rennes, le 19 décembre 2023
Unité de recherche : Centre Inria de l'université de Rennes

### Rapporteurs avant soutenance :

Daniel RACOCEANU Professeur, Sorbonne Université

David ROUSSEAU Professeur, Université d'Angers

### Composition du Jury :

Attention, en cas d'absence d'un des membres du Jury le jour de la soutenance, la composition du jury doit être revue pour s'assurer qu'elle est conforme et devra être répercutée sur la couverture de thèse

Président : Prénom NOM Fonction et établissement d'exercice (à préciser après la soutenance)

Rapporteurs : Daniel RACOCEANU Professeur, Sorbonne Université

David ROUSSEAU Professeur, Université d'Angers

Examinatrices : Maria VAKALOPOULOU Maitre de conférence, Centrale Supélec

Ruth WILKINS Adjunct Research Professor, Carleton University

Elisa FROMONT Professeure, Université de Rennes

Dir. de thèse : Charles KERVRANN Directeur de recherche, Inria Rennes

Co-dir. de thèse : Mohamedamine BENADJAOUD Chercheur, IRSN / SERAMED
<figure><img src="image_1.png" /></figure>
# TABLE OF CONTENTS

Notations and acronyms 17

Résumé en français 21

Introduction 31

1 Deep Learning for image analysis 41

1.1 Deep Learning architectures 42

1.1.1 ImageNet and deep learning 42

1.1.2 Deep neural networks 42

1.1.3 Convolutional Neural Networks 47

1.2 Unsupervised learning for image analysis 52

1.2.1 Generative models 53

1.2.2 Similarity-based methods 59

1.3 Approximate Bayesian deep learning 65

1.3.1 SGD as an Orstein-Uhlenbeck process 66

1.3.2 Ensembling neural networks for performance and uncertainty estimation 68

1.4 Object detection with deep learning models 73

1.4.1 Object detection with Region-CNN 73

1.4.2 NMS-free object detection: Object as Points 77

1.5 Conclusion 78

2 Image analysis for the automation of biological dosimetry 79

2.1 Basics of biological dosimetry 79

2.2 Challenges of automated aberration counting 82

2.3 An historic example of automated dicentric scoring systems 84

2.4 ADCI: implementing an ADS framework 87

2.4.1 Chromosome feature extraction 87

2.4.2 Chromosome classification 88
TABLE OF CONTENTS

<table><tr><td>2.4.3</td><td>A pipeline for dose estimation</td><td>90</td></tr><tr><td>2.5</td><td>Deep learning for biological dosimetry.</td><td>94</td></tr><tr><td>2.6</td><td>DCScore: evaluating an ADS in a semi-automatic regime</td><td>97</td></tr><tr><td>2.6.1</td><td>Evaluating DCScore in realistic scenarios</td><td>97</td></tr><tr><td>2.6.2</td><td>Exploring DCScore shortcomings</td><td>99</td></tr><tr><td>2.7</td><td>Conclusion</td><td>102</td></tr></table>

3 Two-stage chromosomal aberration detection with patch classification 105

<table><tr><td>3.1</td><td>Datasets</td><td>106</td></tr><tr><td>3.1.1</td><td>Patch dataset</td><td>106</td></tr><tr><td>3.1.2</td><td>Skeleton datasets</td><td>106</td></tr><tr><td>3.2</td><td>Chromosome classification</td><td>109</td></tr><tr><td>3.2.1</td><td>Resnet-based classifier</td><td>109</td></tr><tr><td>3.2.2</td><td>Chromosome patches autoencoder</td><td>111</td></tr><tr><td>3.2.3</td><td>Latent space classifier</td><td>111</td></tr><tr><td>3.2.4</td><td>Performance results and commentary</td><td>111</td></tr><tr><td>3.3</td><td>Simulation of chromosome patches</td><td>115</td></tr><tr><td>3.3.1</td><td>Simulating chromosomes with VAEs</td><td>115</td></tr><tr><td>3.3.2</td><td>Simulating chromosomes with pix2pix</td><td>118</td></tr><tr><td>3.4</td><td>Training Faster R-CNN on a synthetic dataset</td><td>118</td></tr><tr><td>3.5</td><td>Conclusion</td><td>119</td></tr></table>

4 End-to-end chromosomal aberration detection in metaphase images 121

<table><tr><td>4.1</td><td>Related works</td><td>122</td></tr><tr><td>4.1.1</td><td>Key-point regression in deep learning</td><td>122</td></tr><tr><td>4.1.2</td><td>Object detection and counting</td><td>122</td></tr><tr><td>4.1.3</td><td>Model aggregation</td><td>123</td></tr><tr><td>4.2</td><td>Methods</td><td>124</td></tr><tr><td>4.2.1</td><td>Keypoint regression with heatmap regression models</td><td>124</td></tr><tr><td>4.2.2</td><td>Implicit ensembling of neural networks</td><td>127</td></tr><tr><td>4.2.3</td><td>Setting of model parameters</td><td>128</td></tr><tr><td>4.2.4</td><td>Visualization of the training dynamics of single model</td><td>129</td></tr><tr><td>4.3</td><td>Materials</td><td>132</td></tr><tr><td>4.3.1</td><td>Data description</td><td>132</td></tr><tr><td>4.3.2</td><td>Evaluation metrics</td><td>134</td></tr></table>
<table><tbody><tr><td>4.4</td><td>Experimental results</td><td>136</td></tr><tr><td>4.4.1</td><td>Performance of single model</td><td>136</td></tr><tr><td>4.4.2</td><td>Performance of model ensemble</td><td>139</td></tr><tr><td>4.4.3</td><td>Robustness to non-chromosome objects in metaphase images</td><td>140</td></tr><tr><td>4.4.4</td><td>Visualization of training trajectories</td><td>142</td></tr><tr><td>4.4.5</td><td>Transfer between training and calibration curve datasets</td><td>142</td></tr><tr><td>4.5</td><td>Discussion and future works</td><td>145</td></tr><tr><td colspan="2"><strong>Conclusion</strong></td><td><strong>147</strong></td></tr><tr><td colspan="2"><strong>A Sparsity in optical flow for microfluidics imaging</strong></td><td><strong>151</strong></td></tr><tr><td>A.1</td><td>Introduction</td><td>152</td></tr><tr><td>A.2</td><td>Unsupervised CNN-based method for sparse optical flow estimation</td><td>154</td></tr><tr><td>A.2.1</td><td>Definition of data term and loss function</td><td>154</td></tr><tr><td>A.2.2</td><td>Definition of spatial regularizers</td><td>155</td></tr><tr><td>A.3</td><td>Experimental results and comparison of supervised and unsupervised methods</td><td>157</td></tr><tr><td>A.3.1</td><td>Description of datasets</td><td>157</td></tr><tr><td>A.3.2</td><td>Evaluation metrics for optical flow estimation</td><td>157</td></tr><tr><td>A.3.3</td><td>Description of competing methods</td><td>158</td></tr><tr><td>A.3.4</td><td>Experimental results</td><td>159</td></tr><tr><td>A.4</td><td>Conclusion</td><td>163</td></tr><tr><td colspan="2"><strong>Bibliography</strong></td><td><strong>165</strong></td></tr></tbody></table>
LIST OF FIGURES

1 Chronologie d'une estimation de dose de rayonnement ionisant en dosimétrie biologique. Un échantillon sanguin est mis en culture pendant 48h, et la mitose est inhibée durant la métaphase avec du Demecolcine. Les microscopes modernes permettent une acquisition automatisée des métaphases, comme montré en a). Les experts biologistes comptent les chromosomes dicentriques (entouré en rouge) et identifient également les fragments (en vert) et les anneaux-centriques (en bleu), voir b). Finalement, le nombre moyen d'aberrations par lymphocyte permet d'estimer la dose, comme montré en c). 22

2 A) Figure montrant les variations possibles en terme d'illumination, de longueur de chromosomes et de présence de noyaux et débris cellulaires B) Figure montrant les variations possibles en terme de longueur de chromosome, de position du centromere et de repliement des chromatides 24

3 Prédictions d'un ensemble de quatres modèles différents. Le premier modèle ne détecte pas le fragment de chromosome (en vert) en bas à gauche de l'image. Les trois autres modèles de l'ensemble le détectent correctement. Les quatres modèles prédisent le chromosome dicentrique près du centre de l'image correctement. 28

4 Chronology of a biological dosimetry-based estimation of an ionizing radiation dose. First, a blood sample is taken and cells are grown for 48 hours, Demecolcine is used to stop cells in the metaphase stage of mitosis. Metaphase acquisition is automated with modern microscopy systems, as shown in a), and specialists count DCs (circled in red), fragments (circled in green) and ring-centric chromosomes (circled in blue) as displayed in b). The proportion of chromosomal aberration is linked to the ionizing radiation dose through a linear-quadratic relationship, see c). 32
<table>
  <tbody>
    <tr>
      <td>1.12</td>
      <td>SimSiam architecture. Note that the dimension of z and p are identical.</td>
      <td>65</td>
    </tr>
    <tr>
      <td>1.13</td>
      <td>In [58], parameter space is explored using a cyclical learning rate schedule. An ensemble is sampled from saved checkpoints.</td>
      <td>69</td>
    </tr>
    <tr>
      <td>1.14</td>
      <td>L<sub>2</sub>-regularized cross-entropy train loss (left) and test error (middle) as a function of the interpolation parameter t between local minimas. Straight lines between local minimas showcase a large increase in loss, while a Bezier interpolation between both local minimas shows the presence of a low-loss path between those models. "Polychain" (a polygonal chain between minimas) shows similar behavior.</td>
      <td>70</td>
    </tr>
    <tr>
      <td>1.15</td>
      <td>Comparison of the sampling of the posterior distribution for several methods. An ensemble of several DNNs trained separately [57] effectively provides several true samples of the posterior distribution. However, the number of samples is low. SWAG can cover a local minima of the loss function (a mode of the posterior distribution) much more densely, but can only cover one mode. Multi-SWAG covers several modes of the posterior distribution by fitting several Gaussian distributions during training, to form a mixture model of the posterior distribution.</td>
      <td>72</td>
    </tr>
    <tr>
      <td>1.16</td>
      <td>For a proposed region, features from a ConvNet are retrieved. An adaptive pooling mechanism called ROI pooling ensures that features have the same spatial size, no matter the size of the proposed region. Pooled features are used for classification and bounding box regression. (Source: [65])</td>
      <td>74</td>
    </tr>
    <tr>
      <td>1.17</td>
      <td>For every pixel of the feature map produced by the backbone, the RPN produces a set of regression vector modifying the pre-defined anchors, and a set of objectness scores which quantify how likely each anchor is to contain an object. (Source: [9])</td>
      <td>75</td>
    </tr>
    <tr>
      <td>1.18</td>
      <td>Architecture of Faster R-CNN. A ConvNet produces a feature map for an input image. This feature map is used to produce region proposals. Those region proposals are selected, and an object detection network is trained as presented in Figure 1.16. (Source: [9])</td>
      <td>76</td>
    </tr>
    <tr>
      <td>1.19</td>
      <td>Output of the CenterNet model. The center of the bounding box is predicted, a specific heatmap is predicted for every class in the dataset. For every location of the heatmap, a (height, width) tuple is predicted, representing the size of the bounding box at this location. (Source: [67])</td>
      <td>77</td>
    </tr>
    <tr>
      <td>2.1</td>
      <td>Basic morphology of a metaphase chromosome.</td>
      <td>80</td>
    </tr>
  </tbody>
</table>
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). 137

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. 139

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. 140

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. 143

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. 144
A.6 Crops of frames 278 (top) and 279 (bottom) of the flow field estimated by RAFT on the RBC sequence. 160

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 indicates 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`. 161

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 indicates the $(q_5, q_{95})$ 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`. 162
# NOTATIONS AND ACRONYMS

## Notations

### General

- $\mathbb{R}$: set of real numbers

- $1_K$: vector where every component has value 1 of size $K$

- $x^T$: transpose of $x$

- diag($x$): square matrix with diagonal $x$

- $\|x\|_2$: $L_2$-norm of $x$

- $\theta$: parameters of a model

- $\mathbb{E}$: expectation of random variable

- $\mathbb{KL}(X|Y)$: Kullback-Leibler divergence between distributions $X$ and $Y$

- $\nabla_{\theta}$: gradient with respect to variable $\theta$ (usually model parameters)

- $\mathcal{H}_{\theta}$ Hessian with respect to $\theta$

- $\mathcal{D}$: training dataset

- $\mathcal{J}$: batch of training data (random subset of $\mathcal{D}$)

- $|A|$: number of elements in set $A$, for example $|\mathcal{J}|$ is the batch size

- $\sigma$: activation function of a neural network

- $\frac{\partial f}{\partial x}$: partial derivative of $f$ with respect to variable $x$

- $\mathcal{L}$: Risk criterion

- $H, W$: Height and Width of an image

### Unet ensembling

- $\mathcal{R}$: Sparse Variation regularization term

- $\rho$: Sparse Variation sparsity parameter

- $\lambda$: regularization strength parameter

- $M$: number of models in ensemble

- $(u, v) \in \Omega$: pixel coordinates ($\Omega$) is the set of all pixel coordinates, with $|\Omega| = H \times W$

- $\phi_{\theta}$: Unet model with parameters $\theta$
— $y_L$: image downsampled by a factor $L$, i.e of size $H/L, W/L$

— $\nabla_{u,v}$: intensity gradients of image (with respect to image axes)

— $\Lambda_{\theta_k}(x_i)$: binarized prediction of $\phi_{\theta_k}$ on image $x_i$

— $\Lambda_{GT}(y_i)$: binarized ground truth

— $I_i$: intersection between binarized prediction and ground truth

— $T_A$: agreement threshold

— $T_C$: single-model confidence threshold

— $S_i(u, v)$: sum of the binary predictions of ensemble members at location $u, v$ for image $x_i$

— $D_i(u, v)$: binary aggregated decision at location $u, v$ for image $x_i$

Acronyms

— IAEA: International Atomic Energy Agency / Agence Internationale pour l'Energie Atomique

— DC: Dicentric Chromosome

— MC: Monocentric Chromosome

— LRAcc: Laboratory for Radiobiology of Accidental Exposures / Laboratoire pour la Radiobiologie des Expositions Accidentelles

— IRSN: Radioprotection and Nuclear Safety Institute / Institut pour la Radioprotection et la Sûreté Nucléaire

— SERPICO: Space-TimE RePresentations, Imaging and cellular dynamics of Molecular COmplexes

— Inria: National Institute for Research in Informatics and Automatics / Institut National pour la Recherche en Informatique et Automatique

— SSL: Self-Supervised Learning

— ADCI: Automated Dicentric Chromosome Identifier

— ADS: Automated Dicentric Scoring

— DL: Deep Learning

— CNN: Convolutional Neural Network

— R-CNN: Region CNN

— ConvNet: Convolutional Neural Network

— FCN: Fully Connected Network

— DNN: Deep Neural Network
- SVM: Support Vector Machine

- SGD: Stochastic Gradient Descent

- MCMC: Monte Carlo Markov Chain

- GAN: Generative Adversarial Networks

- VAE: Variational AutoEncoders

- IoU: Intersection over Union

- ROC: Receiver Operating Characteristic

- AUC: Area Under the Curve

- MCMC: Monte Carlo Markov Chain

- UMAP: Uniform Manifold Approximation and Projection
<figure><img src="image_2.png" /></figure>
# RÉSUMÉ EN FRANÇAIS

## Contexte de la thèse

### La dosimétrie biologique

La dosimétrie est la discipline qui vise à estimer les doses de radiation ionisantes reçues par les être vivants ou les objets. Dans le cas où un dosimètre est disponible, l'estimation peut en général être réalisée de manière rapide et fiable par lecture directe. Si le contexte d'exposition est connu, la dose peut être reconstruite par simulation Monte-Carlo. Cependant, le scénario d'exposition peut être inconnu dans le cas d'une exposition accidentelle. Dans ce contexte, la dosimétrie biologique propose des méthodes alternatives pour estimer ces doses. Les rayonnements ionisants ayant des effets délétères sur les cellules du corps humain, les traces résiduelles peuvent être exploitées pour estimer la dose initialement reçue.

Actuellement, la méthodologie recommandée par l'Agence Internationale pour l'Energie Atomique (AIEA) en matière de dosimétrie biologique cytogénétique est le comptage de chromosomes dicentriques dans les lymphocytes périphériques. Plus exactement, le biomarqueur d'intérêt est le nombre moyen de chromosomes dicentriques par cellule. Même si les rayonnements ionisants produisent d'autres types d'aberration chromosomique, comme les fragments de chromosome ou les chromosomes en anneau, le comptage de ces derniers ne fait pas partie du protocole défini par l'AIEA. Au cours d'un examen, un échantillon sanguin est collecté, et les cellules sanguines sont mises en culture pendant 48 heures. La mitose est inhibée durant la métaphase avec du Demelcocine<fnref n="1" />. Ensuite, l'échantillon sanguin est fixé sur une lame, coloré en Giemsa et examiné avec un microscope optique. Le Giemsa est un colorant couramment utilisé en histologie, composé de bleu de méthylène et d'éosine. Il adhère spécifiquement au groupes phosphate de l'ADN. Par conséquent, c'est un outil couramment utilisé en cytogénétique pour la visualisa-

<footnote n="1">1. Le Demecolcine (également connu sous le nom de colcemid) est un médicament utilisé en chimiothérapie. Le Demelcocine inhibe la mitose pendant la métaphase en limitant la formation des microtubules. Ce médicament permet de synchroniser les cellules cancéreuses en métaphase afin de maximiser leur radio-sensibilité. Source: https://en.wikipedia.org/wiki/Demecolcine.</footnote>
<figure><img src="image_3.png" /></figure>

Figure 1 – Chronologie d'une estimation de dose de rayonnement ionisant en dosimétrie biologique. Un échantillon sanguin est mis en culture pendant 48h, et la mitose est inhibée durant la métaphase avec du Demecolcine. Les microscopes modernes permettent une acquisition automatisée des métaphases, comme montré en a). Les experts biologistes comptent les chromosomes dicentriques (entouré en rouge) et identifient également les fragments (en vert) et les anneaux-centriques (en bleu), voir b). Finalement, le nombre moyen d'aberrations par lymphocyte permet d'estimer la dose, comme montré en c).

Durant l'étape de métaphase, les chromosomes des lymphocytes sont visibles. Pour cette raison, une image des 23 paires de chromosomes appartenant à un lymphocyte en métaphase sera appelée une "métaphase" par les spécialistes du domaine, un terme que nous reprenons dans la suite du manuscrit. Un expert biologiste examine ensuite la métaphase afin de compter les chromosomes dicentriques. Ces aberrations sont rares, même pour une exposition à une dose élevée. Par exemple, en l'absence d'exposition, environ une cellule sur mille contient un chromosome dicentrique. Une fois ce nombre moyen d'aberrations calculé, une courbe de calibration permet de déterminer la dose correspondante. La Figure 1 présente un résumé du protocole utilisé en dosimétrie biologique par cytogénétique.

Une formation de plusieurs mois est nécessaire pour habiliter une personne à compter fiablement les aberrations chromosomiques. Les chromosomes sont des objets déformables, et certains chromosomes monocentriques peuvent être déformés de telle sorte qu'ils ressemblent à des chromosomes dicentriques. Si une population d'individus importante est exposée à une source de radiation, le comptage manuel des aberrations pour chaque individu est impossible. Le développement d'une solution automatisée est l'objet principal de cette thèse pluridisciplinaire, réalisée dans le cadre d'une collaboration entre deux laboratoires: le LRAcc (IRSN) et SERPICO (Inria). Le Laboratoire pour la Radiobiologie des expositions Accidentelles (LRAcc) est l'un des rares laboratoires français détenant une expertise en dosimétrie biologique. Il fait partie de l'Institut pour la Radioprotection et la
Sûreté Nucléaire (IRSN), et est considéré comme l'expert de référence pour la dosimétrie biologique en conditions accidentelles. L'expertise en matière d'imagerie biologique est fournie par l'équipe Space-timE RePresentation, Imaging and cellular dynamics of molecular COmplexes (SERPICO) de l'Institut National pour la Recherche en Informatique et Automatique (Inria) de Rennes.

Cahier des charges pour un outil de dosimétrie biologique automatisée

La conception d'un système de détection automatique de chromosomes dicentriques présente de nombreuses difficultés. Par exemple, l'adhésion du colorant Giemsa au chromosomes présente une certaine variabilité, ce qui affectera le contraste entre les chromosomes et l'arrière plan. Même si la mitose est inhibée en métaphase par le Demelcocine, les cellules ne sont pas toutes figées au même instant, ce qui induit une certaine variabilité inter-cellulaire. Cela se traduit par des variations de longueur, d'épaisseur et de texture des chromosomes (voir Figure 2 B)). Lorsque l'échantillon sanguin est fixé à la lame, certaines cellules sont proches du rétroéclairage du microscope et d'autres sont plus éloignées, ce qui induit des variations d'illumination dans les images. L'étalement de l'échantillon sur la lame est par ailleurs un procédé aléatoire, sans garantie de produire des métaphases suffisamment étalées pour pouvoir compter les chromosomes et identifier les dicentriques (voir Figure 2 A)). Enfin, la déformabilité des chromosomes produit une variabilité morphologique extrême qui complique la classification automatique des chromosomes. Un système de détection automatique des aberrations se doit d'être robuste au facteurs de variation visible sur la Figure 2.

Finalement, le taux de base d'aberrations impose un taux de faux positifs très faible. Tout taux de faux positifs dépassant 1 pour 1000 conduira nécessairement à une surestimation de la dose. Par conséquent, le système doit être capable de mesurer l'incertitude de prédiction: les faux positifs doivent être interprétés comme des détections incertaines. Sur la base de ces constats, cette thèse vise à utiliser les avancées récentes de l'apprentissage statistique et les réseaux de convolution profond pour élaborer un système de détection automatique d'aberrations plus performant que les systèmes commerciaux actuels.
<figure><img src="image_4.png" /><figcaption>Figure 2 – A) Figure montrant les variations possibles en terme d’illumination, de longueur de chromosomes et de présence de noyaux et débris cellulaires B) Figure montrant les variations possibles en terme de longueur de chromosome, de position du centromere et de repliement des chromatides</figcaption></figure>
L'apprentissage profond en vision par ordinateur

Dans un article fondateur pour la vision par ordinateur, Krizhevsky et al. [1] on dé-
montré il y a une dizaine d'années qu'il était possible d'obtenir des performances élevées
sur ImageNet [2] (un jeu de données de classification) avec un réseau de neurones convo-
lutionnel (appelé AlexNet) dont les paramètres sont appris en minimisant une entropie
croisée via un algorithme de descente de gradient stochastique. Depuis, les réseaux convo-
lutifs ont fait l'objet de recherches extrêmement actives dans la communauté scientifique
de la vision par ordinateur pour résoudre des tâches comme la classification, la détection
d'objets, la segmentation d'image, l'estimation du flot optique, etc.

Les grandes bases de données annotées telles que ImageNet [2], Common Objects in Context (COCO) [3] ou encore Pascal VOC [4] se sont révélées crucales pour la recherche sur les réseaux convolutionnels. En effet, ces modèles nécessitent en général d'importantes quantités de données pour atteindre un haut niveau de performance. La performance progresse presque continûment avec la taille de la base de données, ce qui n'était pas nécessairement une caractéristique des modèles disponibles auparavant. Finalement, ces bases de données permettent d'évaluer différentes architectures de manière objective, et facilite les comparaisons entre les méthodes proposées par la communauté scientifique.

Cependant, ces grandes bases de données annotées sont plus rares dans le contexte
de l'imagerie biomédicale pour plusieurs raisons. Les annotations sont particulièrement
chronophages (et donc coûteuses) car les objets à détecter dans les images de microscopie
ou de scanner sont souvent ambigus. Construire des grandes bases de données annotées
peut également poser des problèmes d'ordre juridique relatifs à la préservation de la vie
privée. La diversité des modalités d'imagerie et des contextes biologiques rend également
la construction d'une base de données "universelle" (similaire à ImageNet) impossible. Par
conséquent, pré-entraîner les réseaux convolutifs sur des données sans annotation est une
préoccupation majeure en apprentissage machine. Ce besoin est au coeur d'un champ de
recherche appelé "apprentissage auto-supervisé", ou *Self Supervised Learning* (SSL). En
apprentissage auto-supervisé, l'objectif est de construire des critères d'apprentissage qui
ne dépendent d'aucune annotation. Une fois qu'un modèle est pré-entraîné sur une grande
base de données sans annotation (en général plus facile à obtenir), il est possible d'utiliser
ce modèle comme initialisation pour un apprentissage sur une base de données supervisée.
Ce principe est appelée apprentissage par transfert (*transfer learning*).

Un autre problème majeur des les réseaux convolutifs est la quantification des incertitudes. En effet, pour une tâche de classification, la probabilité prédite ne correspond pas
au taux d'erreur observé sur cette classe [5]. Par exemple, un modèle prédisant une classe avec une confiance moyenne de 99% aura en général un taux d'erreur supérieur a 1% sur cette classe. En d'autre mots, la probabilité prédite ne peut pas être utilisée comme indicateur de confiance du modèle. Cette caractéristique des modèles profonds est particulièrement problématique dans le contexte de l'imagerie médicale, ou les conséquences d'une erreur du modèle peuvent être très importantes. Construire des indicateurs de confiance fiables est donc un domaine de recherche actif. Le paradigme Bayésien est une manière naturelle de construire ces indicateurs de confiance, mais l'échantillonnage de la distribution *a posteriori* des poids avec les méthodes habituellement utilisées (e.g. méthode de Monte Carlo) n'est pas une option réaliste pour les réseaux convolutionnels profonds, en raison du grand nombre de paramètres. Une approche alternative est celle proposée par l'apprentissage profond Bayésien approximatif. La distribution *a posteriori* est supposée Gaussienne, et plusieurs techniques ont été développées pour estimer les paramètres sous-jacents. Par exemple, dans [6] les itérés de la descente de gradient sont considérés comme des échantillons de la distribution *a posteriori* (un résultat démontré dans [7]), et utilisés pour estimer la moyenne et la matrice de covariance de la distribution *a posteriori* approximative. Cette distribution peut être échantillonnée pour construire un ensemble de modèles afin de quantifier l'incertitude de la prédiction.

## Contributions de cette thèse

La première partie de cette thèse (Chapitres 1 et 2) est consacrée à une présentation (probablement incomplète) de l'état de l'art de l'apprentissage profond en vision par ordinateur et en dosimétrie biologique. Dans le Chapitre 1, nous présentons les éléments de base de l'apprentissage profond, des réseaux convolutionnels, de l'apprentissage auto-supervisé et de l'apprentissage profond Bayésien approximatif. Dans le Chapitre 2 nous présentons un état de l'art de la dosimétrie biologique automatisée. Dans ce contexte, le logiciel Automated Dicentric Chromosome Identifier (ADCI) fait l'objet d'une attention particulière, car c'est la solution commerciale dont le fonctionnement est le mieux documenté, à travers une série d'articles publiés de 2012 à 2019. Nous présentons également l'évaluation de DCScore (une autre solution commerciale de détection de chromosomes centriques) dans des contextes variés. Finalement, des avancées très récentes en matière de dosimétrie automatisée, reposant sur l'apprentissage profond sont présentées. En décrivant le fonctionnement de ces solutions, ce chapitre permet également de présenter les points
de difficulté majeurs dans la conception d'un tel système, dont les prémices sont présentés dans le Chapitre 3.

Les premières tentatives pour réaliser un tel système de détection d'aberrations chromosomiques reposent sur un schéma en deux étapes: la détection de tous les chromosomes, puis la classification de ces derniers en deux classes (chromosomes monocentriques et chromosomes dicentriques). Dans ce système analogue à ADCI, nous utilisons un réseau convolutionnel profond appelé ResNet [8] comme classifieur monocentrique vs dicentrique. Ce classifieur est entraîné sur une base de données de petites images annotées (appelées patch). Ce classifieur permet d'atteindre des niveaux de performance élevés sur cette tâche de classification binaire sans grande difficulté. Cependant, nous ne sommes pas parvenus à reproduire ce niveau de performance avec des patches extraits par des modèles de détection d'objets simples, comme un clustering K-Means. Nous avons identifié deux raisons principales à cet échec: i) la base de données n'était pas suffisamment représentative, ii) le modèle de détection d'objets n'était pas suffisamment performant. Les patches à classer contenaient parfois plusieurs chromosomes, ou des fragments de chromosomes entiers qu'il n'était pas possible de classer correctement. La détection imparfaite des chromosomes dégrade *in fine* assez vite la performance de classification

Nous avons également tenté d'utiliser cette base de données de patches pour apprendre à simuler des chromosomes afin de pré-entraîner des modèles de détection d'objets plus sophistiqués (basés sur l'apprentissage profond) sur des données synthétiques. Construire un simulateur de chromosomes performant s'est révélé particulièrement difficile en définitive, et nous n'avons pas réussi à concevoir une solution véritablement satisfaisante. Afin de tout de même pouvoir évaluer le comportement de modèles de détection d'objets sur des données simulées, nous avons construit des métaphases synthétiques annotées (pour lesquelles les chromosomes dicentriques sont localisés par des boites englobantes) à partir de "squelettes" de chromosomes. Même si ces métaphases synthétiques ne sont pas particulièrement réalistes, nous les avons utilisées pour évaluer Faster R-CNN [9] dans un contexte de détection d'objets rares (tel que la dosimétrie biologique). La rareté des chromosomes dicentriques induit une dépendance forte entre la performance du modèle et le nombre d'images utilisées pour entraîner le modèle. En effet, pour les doses les plus faibles, une base de données d'apprentissage comprenant plusieurs milliers d'images ne contient que quelques exemples de chromosomes dicentriques différents. Nous avons également constaté qu'un modèle entraîné sur des données associées à un scénario de dose forte (i.e., avec un nombre moyen élevé de chromosomes dicentriques par métaphase) s'avérait moins
<figure><img src="image_5.png" /><figcaption>Figure 3 – Prédictions d'un ensemble de quatres modèles différents. Le premier modèle ne détecte pas le fragment de chromosome (en vert) en bas à gauche de l'image. Les trois autres modèles de l'ensemble le détectent correctement. Les quatres modèles prédissent le chromosome dicentrique près du centre de l'image correctement.</figcaption></figure>

performant lorsqu'il est évalué sur des données correspondant à un scénario de dose plus faible. Sur-représenter les chromosomes dicentriques pour faciliter l'apprentissage tend à augmenter le nombre de faux positifs pour les faibles doses. La simulation nous a permis d'identifier ces problèmes et d'élaborer une autre stratégie plus efficace.

La contribution principale de cette thèse est une preuve de concept en dosimétrie autotomatisée qui exploite une grande base de données annotée, pour laquelle les centres de chromosomes dicentriques et de fragments sont localisés dans les métaphases. Cette preuve de concept repose sur l'architecture Unet [10], un terme de régularisation encourageant une prédiction parcimonieuse [11] et l'algorithme d'optimisation Adam [12]. Notre approache utilise également les résultats récents d'Izmailov et al. [13] et de Mandt et al. [7] en matière d'agrégation d'itérés de la descente de gradient collectés pendant l'apprentissage. Les paramètres du modèle sont collectés à intervalles réguliers au cours de l'apprentissage, et nous construisons un ensemble en tirant au hasard plusieurs vecteurs de paramètres. Chacun de ces modèles atteint un haut niveau de performance, mais l'ensemble est également suffisamment diversifié (comme le montre la Figure 3) pour qu'il soit possible d'être plus performant en agrégant les prédictions de plusieurs modèles qu'en sélectionnant un seul modèle. Les aberrations sont acceptées ou rejetées sur la base de l'accord entre les différents modèles de l'ensemble. Finalement, en utilisant certaines connaissances a priori sur la dosimétrie biologique, nous parvenons à dépasser les problèmes de transfert entre les scénarios de dose identifiés sur les données simulées, et nous démontrons que notre modèle parvient à estimer des courbes de calibration très satisfaisantes sur un jeu de données de test non-annoté de plus de 21 000 images.
# Organisation du manuscript

Le Chapitre 1 présente quelques articles de référence de l'apprentissage profond pour la vision par ordinateur. Ce Chapitre introduit les bases des réseaux de neurones ainsi que des réseaux de convolutions, de la rétropropagation et des méthodes du premier ordre utilisées pour entraîner cette classe de modèles. Les principales architectures en apprentissage auto-supervisé sont également présentées, telles que les autoencodeurs et les approches basées sur la consistence. Finalement, ce chapitre décrit les grands principes de l'apprentissage Bayésien profond approximatif, car c'est une composante majeure de la preuve de concept présentée dans le Chapitre 4.

Le Chapitre 2 est un état de l'art en dosimétrie biologique automatisée utilisant des images colorées au Giemsa. Nous passons en revue deux logiciels conçus pour compter automatiquement les chromosomes dicentriques. Tout d'abord, ADCI est une solution commerciale développée par Cytognomix, une entreprise canadienne. Cette solution est un pipeline "détection-classification" qui n'utilise pas de réseaux convolutionnels. Les articles présentant les composants de cette solution donnent une idée précise de son mode de fonctionnement, ce qui est rare dans le contexte des solutions commerciales en imagerie biomédicale. Ensuite, nous présentons une série d'articles issus principalement du LRAcc qui évaluent DCScore, une autre solution commerciale de détection de chromosomes dicentriques proposée par l'entreprise allemande MetaSystems. Même si DCScore n'est pas suffisamment performant pour être utilisé de manière complètement automatisée, ce logiciel peut être utilisé de manière semi-automatique, les chromosomes dicentriques détectés font l'objet d'une revue manuelle. Finalement, nous présentons les solutions automatisées existantes qui reposent sur le principe de l'apprentissage profond.

Le Chapitre 3 décrit nos travaux portant sur la reconnaissance de chromosomes individuels. Trois approches principales sont décrites: la classification de chromosomes, la synthèse de chromosomes et l'évaluation de modèle de détection d'objets sur métaphases synthétiques. Entraîner un classifieur de chromosomes est relativement aisé avec l'architecture ResNet [8]. Cependant, la performance de ce classifieur sur des vraies métaphases dépend de la représentativité de la base d'apprentissage, ainsi que de la qualité de notre détecteur de chromosomes, et construire un pipeline efficace est une tâche difficile. Finalement, même si nous ne sommes pas parvenus a construire un simulateur de chromosome performant, l'entraînement de Faster R-CNN [9] sur une base de données de "squelette" de chromosomes nous a permis d'étudier le comportement de ce modèle en fonction du volume de
données.

Dans le Chapitre 4 nous présentons notre détecteur de chromosomes dicentrique, qui repose sur une architecture Unet, une méthode d'optimisation stochastique du premier ordre (Adam), un terme de régularisation encourageant des prédictions parcimonieuses et une technique d'agrégation de modèles. L'agrégation permet de quantifier en particulier le niveau de confiance associé à une détection, et de rejeter les détections peu confiantes. Notre méthode d'agrégation repose sur des paramètres interprétables, et atteint un haut niveau de performance en considérant à la fois les métriques de détection d'objets et l'estimation de courbes de calibration.

Communication scientifique

Notre preuve de concept pour un détecteur automatique de chromosomes dicentriques a été décrite dans un article intitulé "Ensembling Unets, sparse representation and low dimensional visualization for rare chromosomal aberration detection in light microscopy images" (https://doi.org/10.1101/2023.09.11.557124), en cours de soumission à une revue internationale. Cette preuve de concept a également été présentée lors de l'International Conference on Radiation Research à Montréal en août 2023.

De juin à septembre 2022, j'ai réalisé un séjour scientifique dans le département En-
gineering Physics de l'université McMaster au Canada, sous la direction de Qiyin Fang.
Ce séjour avait pour objectif d'implémenter des méthodes d'estimation du flot optique
pour l'imagerie cellulaire microfluidique, dans le cadre d'une thèse (Tianqi Hong). Cette
méthode d'estimation du flot optique réutilise le régulariseur parcimonieux utilisé pour
la preuve de concept du Chapitre 4. Ce séjour a débouché sur la rédaction d'un pre-print
qui fait l'objet de l'Annexe A de cette thèse.
# INTRODUCTION

## Context and motivations

### Biological dosimetry

Dosimetry is the branch of health physics dealing with estimating ionizing radiation doses received by living organisms. If a dosimeter is available, the dose can be read directly with this tool. Even if a dosimeter is not available, the dose can be reconstructed using Monte-Carlo simulations if the exposition context is known. However, if the exposition is accidental the exposition context is usually unknown. In this case, biological dosimetry provides alternative tools to estimate doses of ionizing radiation. As radiation has damaging effects on human cells, several biomarkers can be used to estimate this dose only using information available after the exposition.

The current gold standard in cytogenetic biological dosimetry, defined by the International Atomic Energy Agency (IAEA) is Dicentric Chromosome (DC) scoring [14] in peripheral blood lymphocytes. More precisely, the biomarker of interest is the average number of DC per cell, also called dicentric yield. While radiation induces additional aberrations like chromatid fragments (usually abbreviated as "fragment") or ring-centric chromosomes (usually abbreviated as "rings"), they are currently not considered in the protocol defined by the IAEA.

To estimate a dose, a blood sample is taken from a patient and blood cells are grown. After 48 hours, cell mitosis is inhibited with Demecolcine<fnref n="2" />, the blood sample is stained with Giemsa, mounted on a slide, and examined with a brightfield microscope. Giemsa is a well known histological stain composed of methylene blue and eosin. It binds specifically to the phosphate groups of DNA. Because of this, it is commonly used in cytogenetics to visualize chromosomes. With Giemsa staining, the chromosomes of lymphocytes are visible during metaphase. Because of this, images of chromosomes belonging to a metaphase cell

<footnote n="2">2. "Demecolcine (also known as colcemid) is a drug used in chemotherapy [...] During cell division, demecolcine inhibits mitosis at metaphase by inhibiting spindle formation. Medically, demecolcine has been used to improve the results of cancer radiotherapy by synchronising tumour cells at metaphase, the radiosensitive stage of the cell cycle". Source: https://en.wikipedia.org/wiki/Demecolcine.</footnote>
<figure><img src="image_6.png" /></figure>

Figure 4 – Chronology of a biological dosimetry-based estimation of an ionizing radiation dose. First, a blood sample is taken and cells are grown for 48 hours, Demecolcine is used to stop cells in the metaphase stage of mitosis. Metaphase acquisition is automated with modern microscopy systems, as shown in a), and specialists count DCs (circled in red), fragments (circled in green) and ring-centric chromosomes (circled in blue) as displayed in b). The proportion of chromosomal aberration is linked to the ionizing radiation dose through a linear-quadratic relationship, see c).

are called "metaphases" by biological dosimetry experts. A trained specialist examines the
metaphase and detects DCs.

While any cell contains 23 chromosome pairs, *i.e.* 46 chromosomes, most chromosomes
are healthy Monocentric Chromosomes (MCs), even at high doses, and DCs are rare
events. The base rate of aberration (when there is no exposition) is around 1 DC per
1000 lymphocytes. Once the dicentric yield is computed, a calibration curve provides the
underlying dose. Figure 2.2 provides an overview of the chronologized steps of biological
dosimetry.

Training is needed to become proficient at aberration detection, as separating DCs
from normal MCs is difficult. Chromosomes are deformable objects and some MCs may
twist and become indistinguishable from a DC for an untrained expert. In the context
of a large scale exposition, manual aberration scoring is widely acknowledged to be a
bottleneck and automated solutions are required.

This thesis aims to explore recent advances in deep learning for computer vision to
improve on currently available commercial Automated Dicentric Scoring (ADS) systems.
This PhD has been a multidisciplinary endeavour and was completed between two labs.
The Laboratory for Radiobiology of Accidental Expositions (LRAcc) is one of the main
french laboratories researching biological dosimetry. It is part of the Institute for Radio-
protection and Nuclear Safety (IRSN). LRAcc is focused on dose estimation in irradiation
settings. The interest of LRAcc in automating DC detection led to the funding and su-
pervision of this thesis.

As the LRAcc does not have any specific expertise in biomedical imaging, it was
provided by the Space-timE RePresentation, Imaging and cellular dynamics of molecular
COmplexes (SERPICO) team at the National Institute for Research in Informatics and
Automation (Inria).

Challenges and requirements of an automated DC detector

It is well established that building an ADS system is challenging. The Giemsa staining
process is inhomogeneous and staining does not bind to every chromosome evenly. This
leads to variations in intensities values of image pixels. While Demelcocine is used to
stop mitosis, cells are stopped at various degrees of metaphase which leads to variation
in chromosome length, width and texture. As the blood sample is spread on the glass
slide, some cells are very close to the center of the slide, while some others are closer to
the edges. This leads to variations in illumination, and therefore to variations in pixel
intensity values.

Furthermore, chromosome spread is inherently random and there is no way to guaran-
tee that metaphases are sufficiently spread for chromosome counting and DC detection,
as is shown in Figure 5 A). Finally, the deformability of chromosomes generates extreme
morphological variety (as shown in Figure 5 B)), which complicates automated chromo-
some classification. Any effective ADS system needs to be robust to the variations in
the chemical process, the changes in illumination and the morphological variations of
chromosomes, as shown in Figure 5.

Finally, the base rate of chromosomal aberration imposes extremely strict constraints
on the rate of false positives. Any system that produces more than one erroneous detection
every thousand cell overestimates the dose of a patient that was not irradiated. Because
of this, model uncertainty should be quantified: false positives should be low-confidence
detections. Accordingly, we investigated recent advances in deep learning to improve on
the performance of currently available commercial ADS solutions.

Deep learning for computer vision

Since Krizhevsky et al. [1] demonstrated state-of-the-art performance on the ImageNet image classification dataset [2] using a Convolutional Network (ConvNet) model called AlexNet trained with Stochastic Gradient Descent (SGD), ConvNets have been
<figure><img src="image_7.png" /><figcaption>Figure 5 – A) Figure showcasing possible variations in illumination, chromosome length, chromosome debris and presence of nuclei and debris B) Figure showcasing variations in chromosome length, centromere location and folding of chromatides</figcaption></figure>
an extremely active area of research in computer vision. This class of model has been
used for tasks like object detection, object segmentation, image classification, optical flow
estimation, etc.

Large annotated benchmark datasets like ImageNet [2], Common Objects in Context (COCO) [3] and Pascal VOC [4] have made research on ConvNets for computer vision considerably easier. First, a large dataset is usually needed to reach a high level of performance with ConvNets. It is worth noting that performance keeps improving as the dataset grows, which was not always a given with previous image analysis models. Second, those datasets are used as standardized benchmarks, which makes objective comparison between methods possible.

However, biomedical imaging lacks benchmark datasets like those available in com-
puter vision for several reasons. Labels are much more expensive as they usually require
the expertise of trained medical experts and the relevant objects are harder to locate and
classify in microscopy or scanner images. Furthermore, curating such a dataset may run
into privacy issues. Finally, the wide variety of imaging modalities and biological settings
makes building an "universal" dataset unfeasible.

Therefore, the ability to pretrain ConvNets on un-labelled data is crucial. This re-
quirement has motivated works in Self-Supervised Learning (SSL) [15]. SSL focuses on
designing training criterion that do not use any label. Once a model is trained on a large
unlabelled dataset it can be fine-tuned on a smaller labelled dataset, a process called
"transfer learning".

Another issue facing ConvNets is uncertainty quantification. ConvNets trained on clas-
sification tasks are usually overconfident [5] which means that the predicted probability
does not match the error rate. A ConvNet that predicts a class with a probability of 99%
has an error rate greater than 1% on this class. In other words, predicted probability can-
not be used as a metric for model confidence. This unwanted characteristic of deep models
is especially problematic in medical imaging, where the consequences of a classification
error can be catastrophic.

Accordingly, building accurate confidence metrics for ConvNets becomes an active area
of research. The Bayesian framework is a natural way of modelling model uncertainty, but
sampling the posterior distribution over the networks weight with conventional methods
(like Monte Carlo Markov Chain) is usually computationally intractable, because of the
very large number of parameters in modern ConvNets.

Alternative solutions have been investigated. In approximate Bayesian deep learning,
the posterior distribution over the weights is usually assumed to be Gaussian, and several techniques have been developed to estimate the parameters of this distribution. A noteworthy example is the one provided by Maddox *et al.* in [6]. In this paper successive iterates of the network weights during SGD are taken as approximate samples of the posterior distribution over the weights, as suggested by Mandt *et al.* in [7]. Those samples are used to estimate the mean and the variance of the approximate posterior distribution. This approximate posterior distribution can be used to sample an ensemble of models which can be used for uncertainty quantification.

## Contribution of this thesis

The first part of this thesis (Chapter 1 and 2) provides an introduction to deep learning and automated biological dosimetry. Chapter 1 includes the basics of deep learning, ConvNets, unsupervised learning, and approximate Bayesian deep learning. In Chapter 2, we give an overview of automated biological dosimetry. Most importantly, we describe the Automated Dicentric Chromosome Identifier (ADCI), a current commercially available automated dosimetry solution. The components of this solution are described in several papers published between 2012 and 2019.

We present a few papers written by researchers at LRAcc that evaluate DCScore in various contexts. Finally, we mention two papers implementing ADS systems based on deep learning. By describing current ADS solutions, Chapter 2 provides an overview of the major technical difficulties encountered in the design of automated biodosimetry software.

In our first attempts at designing an ADS system, we tried to reproduce a two-step, detection-classification ADS pipeline with a ResNet-based [8] MC-DC classifier. We used a dataset of labelled chromosome patches to train the ResNet binary classifier used in the second step.

While reaching a very high level of performance on this classification step was relatively easy to achieve, reproducing this performance on real world data (i.e., chromosomes extracted from arbitrary metaphases) proved to be much harder. We identified two main reasons for this failure: our dataset was not representative enough, and our object detection model was not good enough, leading to ambiguous patches that could not be classified accurately. In other words, the performance of the detection and classification steps were strongly dependent.

Furthermore, we attempted to synthesize images of chromosomes, so that we could
<figure><img src="image_8.png" /><figcaption>Figure 6 – Prediction diversity for a set of four different models. The first model does not predict the fragment in the bottom left corner, but the three next model do. All four models predict the DC near the center of the image correctly.</figcaption></figure>

pre-train object detection models on simulated data. Designing a performant chromosome simulator proved to be a challenging task and we failed to find a satisfying solution.

In order to investigate the performance of Faster R-CNN, a relatively modern, deep learning-based object detection model, we simulated labelled metaphase images (metaphases where DCs are located with bounding boxes) using chromosome "skeletons". Because of the under-representation of DCs model performance exhibited a strong dependence to the number of training images.

In a low dose scenario, a dataset containing several thousands images might only contain a few examples of DCs. Furthermore, models trained on data corresponding to a high dose scenario (with a high average number of aberrations per cell) performed worse on "low-dose" data. This means that oversampling images containing large number of aberrations to circumvent their under-representation leads to an increase in false positives.

While we did not build a viable ADS system with our chromosome dataset and did not learn to synthesize realistic data, the insights of this simulation study have been helpful to design an alternative strategy in the next chapter.

The main contribution of this thesis is a proof of concept for DC detection and dose estimation. It is trained on large labelled dataset which was built over the course of several months at LRAcc. In this dataset, the centers of all chromosomal aberrations (including fragments and rings) are located in metaphases. Our solution is based on the Unet architecture [10], a sparsity-promoting regularizing term [11] and gradient-based training with Adam [12].

Furthermore, we benefit from recent results in DNN aggregation by Izmailov et al. [13], based on theoretical insight from Mandt et al. [7]. We build an ensemble from a random
sample of the checkpoints collected at the end of each epoch. While each model of this
ensemble reaches a high level of performance, it is also highly diverse as shown in Figure
6.

This can be used to improve performance beyond the single-model baseline by ac-
cepting or rejecting aberration detections based on the agreement between the models
of the ensemble. Finally, using a priori knowledge on biological dosimetry, we overcome
the issues related to the transfer between different dose scenarios already identified on
simulated data. We show that our model is able to estimate very satisfying calibration
curves on a separate unlabelled testing dataset comprised of over 21 000 images.

Organization of the manuscript

In Chapter 1, we provide a review of relevant references in the deep learning litera-
ture for computer vision. The basics of neural networks, ConvNets, backpropagation, and
gradient-based training are introduced.

A few important references of the similarity-based unsupervised learning literature are
provided as examples of the dominant approach in the current state of the art. Finally, the
approximate Bayesian deep learning literature is given an introduction, as it is a major
component of the proof of concept introduced in Chapter 3.

In Chapter 2, the field of automated biological dosimetry in Giemsa-stained images is
introduced. As biological dosimetry remains a niche application in computer vision there
are only a few competing ADS solutions and to the best of our knowledge, all of them
being commercial software. First, a thorough overview of the Automated Dicentric Chro-
mosome Identifier (ADCI) software is given. ADCI is developed by a Canadian company
(Cytognomix). The papers covering ADCI provide detailed insight into the design of a
"pipeline" (detect, then classify) ADS system.

Second, we go over another set of papers written by researchers at LRAcc evaluating
DCScore, another proprietary ADS solution designed by MetaSystems. The performance
of DCScore in terms of DC detection is insufficient for fully automated use, but it can
still be of significant help in a semi-automated setting (where detected DCs are reviewed
manually). Finally, we review several recent solutions using deep learning for DC detection.

Chapter 3, describes our contributions at the chromosome level. Three main ap-
proaches are described: chromosome classification, chromosome synthesis and training
object detection models on synthetic metaphases. In Chapter 4, our proof of concept for
a ConvNet-based ADS system is introduced. Training an effective chromosome classifier is relatively easy with ResNet [8]. However, the performance of this classifier on true chromosomes depends on the distribution gap between real metaphases and our training dataset and on the quality of our chromosome detection algorithm. Building an effective pipeline from this classifier proved to be a difficult task. We did not succeed in building a performant chromosome simulator, but training Faster R-CNN [9] provided insight into the relationship between training dataset size and performance in the context of biological dosimetry.

Chapter 4 introduces our proof of concept for DC detection. This solution uses well-known building blocks like Unet, stochastic first-order optimization, a sparsity-promoting regularizer, and model ensembling. Aggregation makes rejection of low confidence detections easy, and relies on interpretable parameters. This ensemble-based approach demonstrates strong performance, both in terms of object detection metrics and in terms of calibration curve estimation.

## Scientific communication

The proof of concept for an ADS system was described in an article titled "Ensembling Unets, sparse representation and low dimensional visualization for rare chromosomal aberration detection in light microscopy images" (https://doi.org/10.1101/2023.09.11.557124), which is currently being submitted a journal. This proof of concept was also presented at the *International Conference on Radiation Research* in Montreal (August 2023).

From June to September of 2022, I stayed in the Engineering Physics department at 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 research stay led to a pre-print which is presented in Annex A.
<figure><img src="image_9.png" /></figure>
# DEEP LEARNING FOR IMAGE ANALYSIS

Deep learning has been at the heart of significant progress in machine learning over the past decade. Image analysis is one of the fields where the influence of Deep Neural Networks (DNNs) has been the most significant. As the application of modern deep learning to chromosomal aberration detection was the driving force behind this PhD project, we present useful references for the rest of the thesis in this chapter.

First, we give a brief introduction of deep learning in the context of image analysis. After introducing Fully Connected Networks (FCNs) and their shortcomings for image data, we give the rationale behind an ubiquitous alternative in image analysis: Convolutional Networks (ConvNets). We introduce backpropagation, and the challenges associated with training Deep Neural Networks (DNNs) using Stochastic Gradient Descent (SGD).

As labelled data is often sparse in biomedical imaging, we provide an introduction to unsupervised deep learning for images. Two main approaches are highlighted: generative models and similarity-based methods.

Providing an accurate modelling of model uncertainty is a key requirement for an automated dosimetry system. Therefore, we provide a review of a few major references in the approximate Bayesian deep learning literature. In those references, the authors use a Gaussian approximation for the posterior distribution over the weights of the networks. Several ways of estimating the mean and covariance of this Gaussian distribution are provided.

Finally, we give two examples of ConvNet-based object detection models, Faster R-CNN and CenterNet. First, we present Faster R-CNN as an example of a two-stage object detector. We focus on the technical complexity of this model. As an alternative, we highlight the simplicity of one-stage object detectors using CenterNet as an example.
1.1 Deep Learning architectures

1.1.1 ImageNet and deep learning

The ImageNet database [2] was published in 2009 to investigate large scale image clas-
sification problems. The ImageNet Large Scale Visual Recognition Challenge (ILSVRC)
ran yearly, to incentivize and publicize improvements in image classification research.
Initially, the best performing methods were popular machine learning techniques at the
time, like Support Vector Machines (SVMs) [16]. Those methods usually relied on con-
vex optimization which means that their convergence could be theoretically guaranteed.
Before the early 2010s, most methods used hand-crafted features, like Histogram of Gra-
dients (HOG) [17] or Scale Invariant Feature Transforms (SIFT) [18]. In a seminal paper
published in 2012, Krizhevsky et al. showed that ConvNets trained with SGD brought a
considerable improvement on the current performance state-of-the-art [1]. This sparked a
renewal of neural networks methods in machine learning, which were previously thought
to be very hard to train, owing to the non-convexity of their loss functions.

The development of ConvNets was enabled by the existence of the ImageNet database but also by improvements in computing, both on the software and hardware front. With an appropriate redefinition of the image that trades some memory for speed, convolution can be seen as a matrix multiplication between a kernel and an image, which is an embarrassingly parallel problem. The Compute Unified Device Architecture (CUDA) framework made it relatively easy to write fast 2D convolution implementations running on Nvidia Graphical Processing Units (GPUs) with hundreds of cores. While the initial implementation of AlexNet relied on custom CUDA kernels (which required a working knowledge of C++), higher level frameworks quickly made deep learning much more accessible by packaging so-called "primitives" (convolutions, Fourier transforms, etc) into simple Application Programming Interfaces (APIs) like PyTorch and Tensorflow [19], [20].

1.1.2 Deep neural networks

DNNs are flexible models composed of stacks of linear and non-linear transformations
of the input data. For a two layer neural network, with some activation function σ and
an arbitrarily shaped input x, the output of this neural network can be written as:

$$
\hat{y} = f_{\theta}(x) = \sigma(W_2^T \sigma(W_1^T x + b_1) + b_2) \quad (1.1)
$$
with $\theta = W_1, W_2, b_1, b_2$ being trainable weights and biases respectively. For $x \in \mathbb{R}^k$, $W_1$ is a matrix of size $k \times m$ where $m$ is the feature dimension in the first layer, while $b_1$ is of size $m$. Here, $m$ is a user-defined parameter. This architecture is arbitrarily flexible: one can add layers, make them larger or smaller or build connections between different layers. The number of successive layers is usually refered as the depth of the neural network. Neural networks are often described as deep, without any agreed-upon standard on the number of layers needed for this description. Because $\sigma$ is a non-linearity (for example, $\sigma(x) = \frac{1}{1+e^{-x}}$ or $\sigma(x) = \max(0, x)$), extremely complex, non-linear relationships between inputs and outputs can be modelled. In fact, a number of universal approximation theorems exist. A famous example was proven by George Cybenko [21]:

**Theorem 1** Let $I_n$ be the $n$-dimensional unit cube $[0, 1]^n$. The space of continuous functions on $I_n$ is denoted as $C(I_n)$. Let $M(I_n)$ be the space of finite signed regular Borel measures on $I_n$. We say that $\sigma$ is discriminatory if for a measure $\mu \in M(I_n)$

$$ \int_{I_n} \sigma(y^T x + \theta) d\mu(x) = 0, \qquad (1.2) $$

for all $y \in \mathbb{R}^n$ and $\theta \in \mathbb{R}$ implies that $\mu = 0$. Let $\sigma$ be any continuous discriminatory function. Then finite sums of the form

$$ G(x) = \sum_{j=1}^{N} \alpha_j \sigma(y_j^T x + \theta_j) \qquad (1.3) $$

are dense in $C(I_n)$. In other words, given any $f \in C(I_n)$ and $\varepsilon > 0$ there is a sum $G(x)$ for which

$$ |G(x) - f(x)| < \varepsilon, \forall x \in I_n. \qquad (1.4) $$

## Training DNNs with first order methods

While Theorem 1 suggest that DNNs are a very flexible class of functions, it tells us nothing about how to approximate a desired function with a neural network i.e., learn the parameters of the model. This is be done by optimizing some loss function $\mathcal{L}$ that depends on data and model parameters. Cross entropy (for a classification problem) or $L_2$ distance (for a regression problem) are common examples of loss functions. For a set of $N$ input and output $\mathcal{D} = \{(x_1, y_1)\dots(x_N, y_N)\}$ sampled from some joint probability distribution $(X, Y)$, we can approximate the true risk $\mathbb{E}_{(X,Y)}(\mathcal{L}(X, Y))$ with the empirical
risk, which is also our training criterion:

$$ \mathbb{E}_{(X,Y)}(\mathcal{L}(X, Y)) \approx \sum_{i=1}^{N} \mathcal{L}(f_{\theta}(x_i), y_i). \qquad (1.5) $$

This empirical risk can be minimized with respect to the model parameters $\theta$, which is called training the model:

$$ \theta^* = \arg \min_{\theta} \sum_{i=1}^{N} \mathcal{L}(f_{\theta}(x_i), y_i). \qquad (1.6) $$

This simple idea has several issues:

- There is no unique solution to Equation 1.6, and no analytical description of the set of solutions. Therefore, $\theta^*$ is be estimated using iterative methods that rely on the efficient computation of $\nabla_{\theta} \sum_{i} \mathcal{L}(f_{\theta}(x_i), y_i)$, as detailed in Section 1.1.2. Note that $\nabla_{\theta} \mathcal{L}(f_{\theta}(x), y) \in \mathbb{R}^k$, with $k$ being the dimension of $\theta$. Therefore, $\mathcal{H}_{\theta} \mathcal{L}(f_{\theta}(x), y) = \nabla_{\theta} [\nabla_{\theta} \mathcal{L}(f_{\theta}(x), y)] \in \mathbb{R}^{k \times k}$. As modern neural networks may feature billions of trainable parameters, it is not even be possible to store the Hessian matrix in memory let alone compute it in any reasonable time. Note that despite this fact, a number of fast approximations of the Hessian of DNNs have been suggested, see the papers of Ritter et al. [22], [23] or Martens et al. [24]. However, the overwhelming majority of deep learning literature still relies on first order methods, like SGD [25] or Adam [12].

- The number of training samples may be extremely large. Internal Google datasets like JFT-3B contain billions of images: computing the gradient of the empirical risk over the complete dataset is not feasible because of memory constraints. Instead, this gradient is approximated on a random subset $\mathcal{J}$ of the training data, i.e $\nabla_{\theta} \sum_{x,y \in \mathcal{D}} \mathcal{L}(f_{\theta}(x), y) \approx \nabla_{\theta} \sum_{x,y \in \mathcal{J}} \mathcal{L}(f_{\theta}(x), y)$, with $|\mathcal{J}| < |\mathcal{D}|$. Thus, model parameters are updated with a noisy estimate of the gradient.

- This training criterion is non-convex; there is no guarantee of finding the global minimum and there is a large number of local minima, usually with comparable loss values. Because of our noisy estimation of the gradient and the fact that we usually do not have any curvature information the exploration of the parameter space is a stochastic process: two successive training runs of the same model with identical initialization but different batch samplings will reach two different parameter values at the end of training. While this was initially seen as a insurmountable issue before
the revival of neural networks in the early 2010s, it has since been suggested that the randomness in parameter space exploration may actually be a key positive factor in deep learning performance. This hypothesis has been named implicit regularization in the literature, see [26], [27], and is discussed further in Section 1.3.

— Neural networks architecture are extremely flexible: the number of trainable pa- rameters may be similar to the number of training samples, or it may even be much larger. This means that there is a large set of values for θ that achieve a training loss of zero. In fact, surprisingly, DNNs can fit the training data perfectly and generalize very well on unseen test data, see [27]. As explained in the previ- ous item, randomness in batch sampling helps to regularize training. Contrary to previous assumptions, "bad" minimas that do not generalize are rare because of implicit regularization but also thanks to architectural bias embedded in convolu- tional architectures as we will see in Section 1.1.3.

Computing gradients in network architectures

To train a DNN with first order methods, we need an efficient way to compute ∇$_{θ}$ L ( f$_{θ}$ ( x ) , y ) . Let us consider the same 2-layer FCN as before:

$$
f_{\theta}(x) = \sigma(W_2\sigma(W_1x + b_1) + b_2), \quad (1.7)
$$

with θ = { W$_{1}$, W$_{2}$, b$_{1}$, b$_{2}$ } , where W$_{1}$, W$_{2}$ , and b$_{1}$, b$_{2}$ denote weights and biases respectively. Let us also define two intermediate quantities at layer l , the output a$_{l}$ of the activation function σ (called the "activation" of the layer):

$$
a_l = \sigma(z_l), \tag{1.8}
$$

and z$_{l}$ an affine transform of the activations at the previous layer, i.e

$$
z_l = W_l a_{l-1} + b_l, \tag{1.9}
$$

Activations are defined in a recursive fashion: the activations a$_{l}$ at layer l are a function of the activations at layer l − 1, composed of an affine transform W$_{l}$a$_{l}$$_{−}$$_{1}$ + b$_{l}$ and a non-linear activation function σ applied element-wise. This recursive definition makes representing this class of functions as a computational graph possible, as displayed in Figure 1.1. This
<figure><img src="image_10.png" /><figcaption>Figure 1.1 – Backpropagation algorithm for a 2 layer neural network. The derivative of every function node (in blue) is computed with respect to the previous one, for all nodes on the path to a leaf (in red). Leaves are trainable parameters.</figcaption></figure>

graph has two types of nodes, either representing function evaluation (in blue) or model parameters (in red) and data (in yellow). Terminal nodes (i.e., nodes that are not used as input for any other node) are called leaves. Usually, leaves are model parameters. Every function node has an associated derivative with respect to its input. During the forward pass, the output of the network with respect to its input is computed. Activations at each function nodes are retained, because they are needed for gradient computation.

During the backward pass, the derivatives of the criterion with respect to trainable parameters are computed for every leaves of the graph. The derivative of this composition of differentiable functions is computed with the chain rule. Derivatives of the loss with respect to a specific leaf are computed by multiplying the derivatives computed at every node along the path to this leaf. For example, the derivative of the loss function with respect to W₂ is

$$ \frac{\partial \mathcal{L}(f_{\theta}(x), y)}{\partial W_2} = \frac{\partial \mathcal{L}(f_{\theta}(x), y)}{\partial f_{\theta}(x)} \times \frac{\partial \sigma(z_2)}{\partial z_2} \times \frac{\partial z_2}{\partial W_2}, \quad (1.10) $$

as seen in Figure 1.1. Note that $\frac{\partial z_2}{\partial W_2} = a_1$: the gradient at note i depends on the activations (i.e the output) of the previous node, which is why all activations at every layer are retained during the forward pass. This immensely flexible representation of chained computations is at the core of autodifferentiation, and offers numerous possibilities:

— we can differentiate the loss function with respect to any terminal node (i.e any
<figure><img src="image_11.png" /><figcaption>Figure 1.2 – 2D convolution with a (3 × 3) kernel.</figcaption></figure>

node that is not an input of another node), the difference between input data and
a trainable parameter is only semantic. This means that we can differentiate with
respect to training data, and even update the input of the model with respect to
some criterion. This has been used to find images maximizing the activation of a
specific neuron [28], or for style transfer [29] ;

- while neural networks can implement very complex functions, the computation of derivatives is always analytic, without any approximation beyond the limitations of representing real numbers as floats. This is extremely efficient (only one evaluation of the function is needed), and enables the use of lower precision floats (float32, or even float16) because the additional precision is not as if derivatives were evaluated numerically, which saves memory and computation time.

- In PyTorch, the computational graph is defined lazily: for each sample $x_i$, a node is added to the graph if and only if a new (differentiable) function is called on the output of a function node of the current computation graph. Nodes can be added and removed based on user-defined conditions, with an automatic adaptation of the backward pass. This means that the computation graph can be different for every sample of the dataset, or modified based on external input. This flexibility is an advantage of lazily defined graphs over compiled graphs, like the one implemented in Tensorflow. However, this flexibility costs some performance.

1.1.3 Convolutional Neural Networks

Fully-connected neural networks architectures do not scale well to images: for an image $x$ of size $H \times W$, a linear transform $W_{\theta}$ of $x$ that does not change the resolution of the input image has $(H \times W)^2$ parameters if we let $W_{\theta}$ be a dense matrix. For $H = W = 256$, $W_{\theta}$ would contain $256^4$ float32, which is not feasible in practice. Implicitely, for any pixel
<figure><img src="image_12.png" /><figcaption>Figure 1.3 – 2D convolution with a multi-channel kernel.</figcaption></figure>

of the image, this large matrix has a trainable coefficient modelling an interaction with
any other pixel. ConvNets make two simplifying assumptions:

- pixels usually exhibit high correlation with their closest neighbors but lower correlation with pixels located far away. Instead of keeping the full matrix, long-range interactions can be removed so that only interactions between neighboring pixels are modelled which means that $W_{\theta}$ is block-sparse ;

- a good correlation model with neighboring pixels should be location-independent: instead of learning a local correlation at every location on the image grid, the assumption is that a single correlation model is valid for the whole image ;

In other words, we learn the parameters of a convolution kernel, a set of parameters computing a local average at every location of the image. For any given image, local content interacts with the convolution filter based on their similarity, as shown in Figure 1.2. Given the immense variety of natural images, if one wants to detect some type of object reliably, a single convolutional filter is not enough: a filter providing a strong response for some orientation and lighting condition might not be activated reliably by the same object in other conditions. Furthermore, objects in natural images exhibit scale variation, so that a filter might be activated by the object if it is sufficiently distant, but it might not provide useful features if the object is closer. For those reasons, ConvNets rely on a large number of convolution filters, stacked in successive "blocks", to provide an expressive and hierarchical representation of the image, as shown in Figure 1.3. Obviously, convolution did not appear in the deep learning literature first, and has been a core component of image analysis for decades. For example, horizontal and vertical gradients can be computed using convolution filters, as shown in Figure 1.4.
<figure><img src="image_13.png" /><figcaption>Figure 1.4 – Left: input image, middle: horizontal gradient, right: vertical gradient.</figcaption></figure>

Efficient implementation of 2D convolution

A naive implementation of convolution might rely on two **for** loops, iterating over the height and width of the input image and multiplying the local intensities of the image with the convolution kernel. However, this fundamental operation can be parallelized very effectively, as described in this section.

Let us define a batch of $\mathcal{J}$ images, with spatial dimensions $(H, W)$ and $N_i$ channels. For some convolution kernel with spatial size $K_H$, $K_W$, a stride $S$ (the step size of the convolution kernel) and a padding of $P$ pixels, the output of the convolution operation has the following spatial size:

$$H_o = \left( \frac{H - K_H + 2P}{S} \right) + 1, W_o = \left( \frac{W - K_W + 2P}{S} \right) + 1. \qquad (1.11)$$

Using this, the batch of images image can be divided into $N_p$ overlapping patches in an unfolded view, which is an array with dimensions $[\mathcal{J}, K_H \times K_W \times N_i, N_p]$. Here, $N_p = H_o \times W_o$ with $H_o$ and $W_o$ being the spatial size of the convolutions output. Note that this representation is redundant because patches are overlapping: some pixels are represented multiple times. This means the unfolded representation takes up more memory than the initial, "folded" representation.

For a 2D convolution taking an image with $N_i$ channels as input and outputting another image with $N_o$ channels, the convolution kernel (i.e the trainable parameters) is an array of size $[N_o, N_i, K_H, K_W]$. This convolution kernel can be flattened to an array of size $[N_o, N_i \times K_H \times K_W]$, so that the unfolded batch of size $[\mathcal{J}, K_H \times K_W \times N_i, N_p]$ and the flattened convolution kernel can be multiplied together. The output batch is of size
<figure><img src="image_14.png" /><figcaption>Figure 1.5 – Unfold-GEMM-Fold implementation of convolution showcased on a 1 channel image.</figcaption></figure>

$[\mathcal{J}, N_o, N_p]$, i.e $[\mathcal{J}, N_o, H_o \times W_o]$, which can be reshaped into $[\mathcal{J}, N_o, H_o, W_o]$.

The local averaging effect of the convolution can clearly be seen in the output shape: the $K_H \times K_W$ term has disappeared, as the output at every spatial location of the input corresponds to the sum of the element-wise product between the patch and the kernel. Furthermore, this representation makes the fact that convolution is a linear operation obvious. In turns, this explains why computing derivatives in ConvNets is easy, as they are mostly stack of linear transforms and element-wise non-linearities.

As shown in Figure 1.5, 2D convolution can be implemented with General Matrix Multiply (GEMM), which makes parallel implementations easy. Because matrix multiplication is an embarrassingly parallel problem, this makes accelerating this computation of GPUs very easy, but increases memory consumption. This increase is the largest limitation in scaling up convolutional models.
## Architectural bias

Every ConvNet is a FCN, but even for small images where FCNs might be feasible, ConvNets are usually much easier to train. This is because the assumptions made in section 1.1.3 introduce architectural bias which makes modelling images easier. To use an alternate wording, the statistical learning community would say that ConvNets are a good prior model to describe images. This bias has been extensively studied, and has led to a number of surprising discoveries:

- Completely random ConvNets can be effective image priors for image restoration tasks, as shown in Deep Image prior [30]. For some noise vector $z$ and a corrupted image $x$, finding $\theta^* = \arg\min_{\theta} \|f_{\theta}(x) - z\|_2^2$ with early stopping yields good performance on inverse problems, even though the model is trained to reconstruct noise;

- Training only the batch normalization parameters, or only a very small subset of the model parameters is sufficient to reach non-trivial levels of performance, as shown in [31], [32];

- ConvNets can be used as a prior for a FCN. In [33], d'Ascoli et al. find that the weights of a trained convolutional neural network can be used as an initialization for a FCN. FCNs initialized that way reach a higher level of performance than purely convolutional architectures.

## Receptive field

In most cases, convolutional architectures use downsampling to aggregate information from different spatial locations in the image, either by using large convolutional kernels without padding, or by pooling pixels, retaining only the maximum or the average of the sub-region of the image. This spatial aggregation means that as we reach deeper layers of the network, a single pixel in feature space aggregates information from a large subset of the input in image space.

More formally, for a fully convolutional neural network $f_{\theta}$, $a_l$ is the output activation volume at layer $l$. As we are in a convolutional architecture, we have $a_l \in \mathbb{R}^{H/K \times W/K \times N_i}$ with $W, H$ the input image height and width, $K$ some downsampling factor that depends on the stride and padding of the successive convolutions and $N_i$ a pre-defined number of convolution kernels (also called channels). For a given feature pixel $a_l[u, v]$ for $u, v \in \{0...H/K - 1\} \times \{0...W/K - 1\}$ (corresponding to a vector of size $N_i$), the receptive
<figure><img src="image_15.png" /><figcaption>Figure 1.6 – Receptive field in ConvNets.</figcaption></figure>

field of this pixel corresponds to the specific subset of the input image $x[U, V]$ (with $U, V$ a set of image coordinates) so that $a_l[u, v] = f_\theta(x[U, V])$, as shown in Figure 1.6. There is a tradeoff between the receptive field and spatial resolution of a ConvNet: as the receptive field grows larger and the model aggregates information from distant locations of the input image, the spatial resolution is reduced: a single pixel models a large part of the input image. Therefore, accurately modelling objects at different scales in the same image requires specific architectural choices.

## 1.2 Unsupervised learning for image analysis

The early successes of deep learning depended on the availability of large annotated datasets like ImageNet [2], COCO (Common Objects in Context) [3] or Pascal VOC [4]. However, curating large annotated datasets is expensive, especially when it requires specialized knowledge, which is extremely common in biomedical imaging. While the deep learning community quickly noticed that pre-training on ImageNet was usually an effective way to achieve good results with limited amounts of data this solution is less effective when there is a large domain gap between ImageNet and the application domain like in satellite or biomedical imaging. Many imaging modalities found in those applications are not even available in ImageNet, like hyperspectral or volumetric imaging.

Therefore, unsupervised training methods that do not rely on annotated data are needed. Those training techniques should learn useful features from an unlabelled dataset $\mathcal{D} = \{x_1, \dots, x_N\}$. In other words, a pretext task is needed, which can be solved without
any labels. This pretext task is solved by training a feature extractor which can then be
reused for supervised downstream tasks. Designing effective pretext tasks is a longstanding
challenge, as different downstream tasks may requires different features. In this section,
we focus on two research directions for unsupervised learning in image analysis: generative
models, and similarity-based approaches.

A generative model is trained by learning a data distribution from a set of samples.
The curse of dimensionality makes it difficult for images, but a number of successful
methods have been proposed, like Generative Adversarial Networks (GANs) or Variational
Autoencoders (VAEs).

In similarity-based approaches, a similarity criterion is optimized between different
augmented views of the same image. Augmentations can be noise, blurring, flipping, color
shifting ... This way, the network learns features that are invariant to those perturbations.
The minimization of those criterions usually have trivial solutions, where the same fea-
tures are predicted for all images, which is called representation collapse. Preventing this
collapse is a central issue in similarity-based methods.

This section is not intended as an exhaustive survey of unsupervised feature learning
for images as thousands of pages could be written on the subject, even by restricting
ourselves to the deep learning literature. We focus on a few impactful architectures to
highlight the similarity, differences and tradeoffs between those training criterions.

1.2.1 Generative models

A number of approaches have tried approaching the unsupervised learning problem
as a distribution learning problem. As images are very high dimensional, non-parametric
distribution learning techniques suffer from the curse of dimensionality: sample density is
extremely low in image space [34]. In the deep learning era, two noteworthy approaches are
GANs [35] and VAEs [36]. While adversarial models have managed to produce impressive
samples [37], the adversarial criterion remains difficult to optimize [38]. On the other
hand, VAEs enjoy a more principled training criterion but fail to match GANs on sample
quality.

Generative Adversarial Networks

GANs learn to transform low dimensional Gaussian samples from a latent distribution
$p_z(z)$ into images with an adversarial training scheme between a generator $G$ and a dis-
criminator $D$. The discriminator tries to predict wether or not any given image is sampled
from the learned generator distribution $p_g(x)$ or from the true data distribution $p_r(x)$.

In a seminal paper, Goodfellow et al. [35] showed that generating images using this
aforementioned adversarial game was possible. A generator function $G(z, \theta) : \mathbb{R}^l \to \mathbb{R}^d$
maps Gaussian samples to data space, with $l < d$. A discriminator $D(x, \phi) : \mathbb{R}^d \to \mathbb{R}$
solves a binary classification task: it predicts wether the input is a sample of the generator
distribution $p_g(x)$ or the real data distribution $p(x)$. Therefore, the discriminator outputs
a single scalar.

The generator $G$ is trained to minimize $\log(1 - D(G(z)))$, i.e. the probability that the discriminator correctly identifies the output of the generator as fake. On the other hand, the discriminator is trained to maximize the probability of correctly classifying real and fake samples. The training objective of a standard GAN, as proposed in Goodfellow et al. is:

$$
\min_G \max_D V(D, G) = \mathbb{E}_{x \sim p_r(x)} [\log D(x)] + \mathbb{E}_{z \sim p_z(z)} [\log(1 - D(G(z)))]. \quad (1.12)
$$

The min-max formulation makes optimizing this criterion difficult. Optimizing the dis-
criminator to convergence before updating the discriminator would be too costly, so most
GAN training schemes alternate updates of the generator and discriminator. Training
can also collapse rather easily, especially in early phases of training. When the generator
provides poor samples, classifying generator output is easy, so that D(G(z)) is very close
to 0. The generator gradient update (where m is the batch size)

$$
\nabla_{\phi} \frac{1}{m} \sum_{i=1}^{m} \log (1 - D(G(z_i))), \qquad (1.13)
$$

is very small. To prevent this, Goodfellow et al. suggest maximizing $D(G(z))$ instead,
but collapse of training because the discriminator is too strong remains an issue in GAN
training. In fact, [38] proved that the gradient updates of the generator tend to 0 as the
discriminator improves, under some conditions on $p_g(x)$ and $p_r(x)$. In [39], the authors
replace the FCNs used in [40], [41] with ConvNets, which improves training stability and
performance. Furthermore, they show that the discriminator of a GAN can be used as a
feature extractor to solve downstream task like image classification. Additional study of
GAN training stability was provided in [42]. In Mirza et al. [41], the training criterion of
(1.12) is rewritten for class-conditional generation:

$$
\min_G \max_D V(D, G) = \mathbb{E}_{x \sim p_r(x)} [\log D(x, y)] + \mathbb{E}_{z \sim p_z(z)} [\log(1 - D(G(z, y))], \quad (1.14)
$$
<figure><img src="image_16.png" /><figcaption>Figure 1.7 – Image-to-image translations capabilities of pix2pix [43].</figcaption></figure>

where *y* denotes labels. While initially, Mirza et al. showed that GANs can be conditioned with class information, more sophisticated conditioned conditional GANs were later demonstrated. For example, in [43] Isola et al. show that GANs can be conditioned with images, for image-to-image translation, as shown in Figure 1.7. This architecture is called pix2pix. The authors use the following criterion:

$$ \min_G \max_D V(D, G) = \mathbb{E}_{x \sim p_r(x)} [\log D(x, y)] + \mathbb{E}_{z \sim p_z(z)} [\log(1 - D(G(z, y)))] + \mathbb{E}_{x, y, z} [\|x - G(z, y)\|_1]. \quad (1.15) $$

The $\mathbb{E}_{x,y,z} [\|x - G(z, y)\|_1]$ term improves sample quality by forcing the generator to produce outputs that match the true input image. The generator *G* is conditioned both on the noise *z* and the label *y*. However, the authors note that there is little stochasticity in the decoder output: sampling new noise vectors does not lead to a diverse output of the generator. Variants of pix2pix where the generator output does not depend on any noise produce almost identical results to the model trained with the criterion in (1.15).

## Deterministic and probabilistic autoencoders

This section introduces autoencoders, first as a compression model and then as a latent variable model of a true data distribution, as introduced in the famous paper by Kingma et al. [36]. The similarities and differences between the two approaches are highlighted, as presented in [44].

Autoencoders rely on a simple idea: learning how to compress an image, and reconstruct it from a low-dimensional latent representation requires learning good image features.
tures. Image compression was frequently used as a pretext task in early deep learning experiments [45]. This class of models can be described in the following way: an input image $x$, with dimension $d$ is compressed to a latent representation $z \in \mathbb{R}^k$, $k < d$ by a parameterized function $f_\phi(x)$. This latent representation is decompressed into a reconstruction of the input image $\hat{x} = g_\theta(z) = g_\theta(f_\phi(x))$. The parameters $(\theta, \phi)$ can be learned by minimizing a reconstruction error, for example the squared euclidean norm:

$$ \theta^*, \phi^* = \arg \min_{\phi, \theta} \|x - g_\theta(f_\phi(x))\|_2^2. \qquad (1.16) $$

If $f_\phi$ and $g_\theta$ are defined as linear transforms and the reconstruction error is the $L_2$ norm, the optimal encoder and decoder correspond to the first $k$ vectors of the Principal Component Analysis (PCA) basis. Using DNNs, deep autoencoders can be trained with backpropagation and first order optimization. Kingma et al. refined deep autoencoders further in [36], by parameterizing a probability model for a random variable $x$ with DNNs. In this paper, the authors suggest modelling $x$ as a function of a latent unobserved variable $z$:

$$ X \sim p_\theta(x|z) = p(x|g_\theta(z)), \qquad (1.17) $$

where $g_\theta$ is a deterministic function mapping the latent space to the data space, usually called a decoder. In the same way, the probability distribution over the latent space is defined as

$$ z \sim q_\phi(z|x) = q(z|f_\phi(x)), \qquad (1.18) $$

where $f_\phi$ is a deterministic function mapping data space to latent space, usually called an encoder. The likelihood (or model evidence) of this model is

$$ p(x) = \int p(x, z)dz, \qquad (1.19) $$

which requires evaluating the joint distribution over the whole dataset. This integral is usually intractable and an alternative formulation is needed. One can show that the log-evidence can be decomposed as

$$ \log p(x) = \mathrm{ELBO}(\theta, \phi) + \mathrm{KL}[q_\phi(z|x)|p(z|x)], \qquad (1.20) $$
where $p(z|x)$ is the true model posterior, defined as

$$p(z|x) = \frac{p(x|z)p(z)}{p(x)}, \qquad (1.21)$$

which again depends on the intractable model evidence. The Evidence Lower BOund (ELBO) can be rewritten as

$$\mathrm{ELBO}(\theta, \phi) = \mathbb{E}_{q_{\phi}(z|x)} [\log p_{\theta}(x|z)] - \mathrm{KL}[q_{\phi}(z|x)|p(z)], \qquad (1.22)$$

where $p(z)$ is a prior on the latent space. Because of the Jensen inequality, the Kullback-Leibler (KL) divergence between the true and parametric posterior is always positive, and the model evidence $p(x)$ is a constant. Therefore, maximizing the ELBO with respect to the parameters $\phi, \theta$ is equivalent to minimizing the KL divergence between the true and parametric posterior. To optimize the lower bound of (1.22), specific probability models for $q_{\phi}$ and $p_{\theta}$ are needed. Usually, in the context of VAEs, Gaussian probability models are chosen:

$$q_{\phi}(z|x) = \mathcal{N}(\mu_{\phi}(x), \sigma_{\phi}(x)), \quad p_{\theta}(x|z) = \mathcal{N}(\mu_{\theta}(z), \sigma_{\theta}(z)), \qquad (1.23)$$

where $\mu_{\phi}, \sigma_{\phi}$ are functions outputting a mean and a variance for the latent probability model, so that $f_{\phi}(x) = (\mu_{\phi}(x), \mu_{\sigma}(x))$ and $g_{\theta}(z) = (\mu_{\theta}(z), \sigma_{\theta}(z))$. The first expectation term in the ELBO must be estimated with samples of the latent distribution. For computational reasons, a single sample is usually used, which makes the decoder effectively deterministic. Therefore, we have

$$\mathbb{E}_{q_{\phi}(z|x)} [\log p_{\theta}(X|z)] = \|x - \mu_{\theta}(\mu_{\phi}(x))\|_{2}^{2}. \qquad (1.24)$$

Let us focus more closely on $z$. Equation 1.23 tells us that the latent probability model is a Gaussian variable. Therefore, we have

$$p_{\phi}(z|x) = \mu_{\phi}(x) + \sigma_{\phi}(x) \odot \varepsilon, \quad \varepsilon \sim \mathcal{N}(0, I), \qquad (1.25)$$

Where $\odot$ is point-wise multiplication (Hadamard product). Therefore, as pointed out by Ghosh et al. in [44], in practice, VAEs are noise-regularized autoencoders. The deterministic latent representation is corrupted with Gaussian noise, and the variance of this noise depends on the input image. Finally, the second term of (1.22) is a KL divergence between
two Gaussian distributions, which can be computed analytically:

$$ \mathbb{KL}(q_{\phi}(z|x)|p(z)) = \frac{1}{2} \left( \|\mu_{\theta}(x)\|_{2}^{2} + d + \sum_{i=1}^{d} (\sigma_{\phi}(x)_{i} - \log \sigma_{\phi}(x)) \right). \quad (1.26) $$

This divergence is minimized if the trained encoder outputs a centered Gaussian with unit variance. This is a regularization term that tends to smooth the latent space of the autoencoder.

In the initial definition of the ELBO, the reconstruction term of (1.24) and the regularization term of (1.26) are equally balanced. However, in practice, training this model can be difficult. Often, some kind of annealing schedule is needed, where the training starts without any regularization, and progressively increases the weight of the KL divergence term of (1.26) in the ELBO. The training criterion of (1.22) can be rewritten to depend on an hyperparameter $\beta$ which balances reconstruction and regularization:

$$ ELBO(\theta, \phi, \beta) = \mathbb{E}_{q_{\phi}(z|x)} [\log p_{\theta}(x|z)] - \beta \mathbb{KL}[q_{\phi}(z|x)|p(z)]. \quad (1.27) $$

Here, $\beta$ can be increased beyond 1, which over-regularizes the model, trading off reconstruction performance for a smoother latent space. A surprising effect of this over-regularization is that it tends to produce disentangled latent spaces, as shown in [46].

In a disentangled latent space, every dimension of the latent space controls a specific type of variation in output space. This property has not received a rigorous definition, because defining semantic variations identified by humans in images is close to impossible. However, visual examples can provide a better intuition of this property, as shown in Figure 1.8.

## Regularized autoencoders

While VAEs provide a complete probability model for images, where traditional density estimation methods (like kernel density estimation) usually fail, practical implementations of this probability model are usually fairly simplistic. Noticing that VAEs are effectively noise-regularized autoencoders, Ghosh et al. compare VAEs to Regularized Autoencoders (RAEs) in [44]. The following criterion is optimized:

$$ \begin{aligned} \mathcal{L}_{RAE} &= \mathcal{L}_{REC} + \beta \mathcal{L}_{Z}^{RAE} + \lambda \mathcal{L}_{REG} \\ &= \|x - g_{\theta}(f_{\phi}(x))\|_{2}^{2} + \beta \|f_{\phi}(x)\|_{2}^{2} + \lambda \mathcal{L}_{REG}. \end{aligned} $$
<figure><img src="image_17.png" /><figcaption>Figure 1.8 – Latent factors learned by β-VAE on celebA: specific dimensions of the latent space capture factors of semantic variations in images. Note that some semantic factors are not disentangled, like age and gender in b). (Source: [46])</figcaption></figure>

We re-use the deterministic autoencoder notations introduced in Section 1.2.1. The $\beta \mathcal{L}_Z^{RAE}$ is used to constrain the latent variable near the origin of the latent space, to prevent uncontrolled "growth" of the latent space which could lead to numerical issues.

The decoder $g_\theta$ is architecturally close to a GAN generator, as it transforms a latent representation into an image. Therefore, common GAN generator regularizer can be used. Tikhonov regularization on the decoder weights is equivalent to weight decay, which is commonly used. Two noteworthy examples for GAN generators are gradient penalty regularization, $\mathcal{L}_{REG} = ||\nabla_{\phi,\theta}g_\theta(f_\phi(x))||_{2}$, and spectral normalization [47], where weight matrices are divided by an estimate of their largest singular values.

Finally, the authors of [44] note that this deterministic formulation prevents any sampling from the latent space, which is a major advantage of VAEs. To remedy this, they fit a Gaussian mixture model on the latent space once model training is finished. For a random sample $z_{new}$ of this mixture model, $g_\theta(z_{new})$ provides a decoded view.

## 1.2.2 Similarity-based methods

While generative pretext tasks can be easy to train (although some of them are fairly unstable, like GANs), there is no guarantee that the features learned during training will be useful for discriminative downstream tasks like classification. Therefore, pretext tasks that are better suited to discriminative tasks have been designed. In this section, we review a broad class of methods that can be summarized as "similarity based": some kind of similarity metric between augmented views of the same batch of images is optimized during training.
<figure><img src="image_18.png" /><figcaption>Figure 1.9 – Conceptual comparison of three contrastive loss mechanisms (empirical comparisons are in Figure 3 and Table 3). Here we illustrate one pair of query and key. The three mechanisms differ in how the keys are maintained and how the key encoder is updated. (a): The encoders for computing the query and key representations are updated end-to-end by back-propagation (the two encoders can be different). (b): The key representations are sampled from a memory bank [50]. (c): MoCo encodes the new keys on-the-fly by a momentum-updated encoder, and maintains a queue (not illustrated in this figure) of keys. (Source: [48])</figcaption></figure>

Contrastive learning

In contrastive learning, the pretext task is discriminative: the model is trained by
minimizing a classification loss. For a batch of images, two views are produced with
random augmentations for each image. The loss is computed over all pairs of images. The
model should attribute high probability to pairs of matching views, and low probability
to other pairs of images. In all contrastive approaches, selection of negative samples is a
major topic: if the set of negative samples is extremely small, trivial features might be
enough to solve the classification problem. On the other hand, enlarging the set of negative
samples brings computational challenges. We review two approaches in detail: simCLR
(Simple Contrastive Learning of Representations) and MoCo (Momentum Contrast) [48],
[49].

In MoCo [48], the authors frame the contrastive learning problem as a dictionary
lookup task, which can be used as a general framework to highlight the difference between
the approaches suggested in the literature. For any image $x \in \mathcal{D} = \{x_1, \cdots, x_N\}$, a query
encoder $f_q$ produces a query $q \in \mathbb{R}^d$. This query is matched to a set of keys $\{k_1, \cdots, k_m\}$.
The dictionary contains a single positive key $k_+$ (which corresponds to another augmented
<figure><img src="image_19.png" /></figure>

Learning criterion

$$
\begin{align*}
\text{criterion for a pair of images:} \quad l(i, j) &= -\log \frac{e^{s_{i,j}/\tau}}{\sum_{k \neq i} e^{s_{i,k}/\tau}} \\
&= \frac{e^{\square}}{\sum_{\square} e^{\square}}
\end{align*}
$$

$$
\text{total learning criterion:} \quad \mathcal{L} = \frac{1}{2N} \sum_k [l(2k-1, 2k) + l(2k, 2k-1)]
$$

Figure 1.10 – In simCLR, $f_{\theta}$ is trained to maximize similarity between different views of every key in the batch (i.e matching keys), and to minimize similarity between non-matching keys

view of the same image), while all the other negative keys $k_-$ correspond to images
$x_j, j \neq i$. for a single key (i.e a single image), the training criterion is:

$$
\mathcal{L}_q = \frac{e^{q \cdot k_+}}{\sum_{k_-} e^{q \cdot k_-}}. \qquad (1.28)
$$

This way, the similarity between views of the same image is maximized, while the sim-
ilarity between different images is minimized. The fact that the model simultaneously
maximize similarity between positive pairs and dissimilarity between negative pairs is the
key ingredient to prevent representation collapse. Contrastive methods for unsupervised
learning mostly differ in how the key encoder and the set of negative keys are defined, as
explained in Figure 1.2.2. In MoCo [48], the query encoder is a moving average of the key
encoder (with the same architecture), i.e

$$
\theta_k = m\theta_k + (1-m)\theta_q, \tag{1.29}
$$

with $\theta_k, \theta_q$ being respectively the parameters of the key and query encoder. In SimCLR [49], The key and query encoder are shared, so that $\theta_k = \theta_q$. In MoCo, the set of keys is continuously updated, but kept much larger than a single batch to avoid making the discriminative task trivial. The set of keys is a queue: the oldest keys are removed at every batch step, while the new batch represents the newest member of the dictionary.

SimCLR [49] avoids queue mechanisms by using a very large batch sizes. This leads to a
very simple architecture, but requires very large batches (around 4k images on ImageNet),
which quickly becomes intractable with large images (like in biomedical imaging). In this
architecture, there is only one single shared model that is trained using a contrastive
loss. There is no dictionary, and for a positive key, the matching key corresponds to an
augmented view of the same image.

Clustering-based approaches

In supervised image classification, we have a class label $y_i$ for every image $x_i$. Even if this label is not available, we could still optimize a classification loss if we had some proxy label that corresponds to semantic characteristics but does not rely on human annotations. This pseudo-label should be retrieved in an automatic fashion, which is hard because this is effectively equivalent to the supervised classification problem.

For datasets like MNIST featuring well-defined factors of variation [51], simple feature
learning methods like PCA might be enough to build relatively well-defined clusters based
on pixel values in an unsupervised fashion. However, this is not possible for rich and
complex datasets like ImageNet [2].

In [52], Caron et al. proposed an iterative strategy that relies on bootstrapped labels and circumvents the need for human annotations. Initially, image embeddings are computed using an untrained ConvNet. Then, pseudo-labels are built from cluster assignments computed with K-Means. Finally, the ConvNet is trained with a conventional cross-entropy loss using those pseudo-labels. After every training epoch, cluster assign-
ments are re-computed, using the embeddings provided by the classification network.

To prevent a trivial solution where all images are mapped to a single cluster, the centroid of an empty cluster is recomputed from the centroid of a non-empty cluster with an additional random perturbation, and K-Means assignments are re-computed. This effectively splits the non-empty cluster into two. Furthermore, to prevent the model from predicting a single representation, the classification loss of a single sample is weighted based on the number of images in the cluster, i.e, if we have K different clusters of sizes {$n_1, \cdots, n_K$}

$$ \mathcal{L}(X, Y) = \sum_{i} \lambda_{i} y_{i} \log s(f_{\theta}(x_{i})), \qquad (1.30) $$

with $\lambda_i = \frac{1}{n_p}$ if sample $i$ belongs to cluster $p$. Therefore, high classification performance on large clusters brings a lower reduction in loss, which prevents representation collapse.

While this approach is effective, cluster assignments are computed offline, which means that it does not scale very effectively to large datasets. After every epoch of training, clusters have to be re-computed with K-Means. In [53], Caron et al. optimize a similar classification loss based on pseudo-labels. However, those pseudo labels are computed in a different fashion: instead of computing assignments using K-Means on feature vectors, assignments are computed using the Sinkhorn-Knopp algorithm [54], as explained below.

A set of B feature vectors $Z = \{z_1, \cdots, z_B\}$ needs to be mapped to a set of K prototypes $C = \{c_1, \cdots, c_K\}$. Feature vectors should be mapped to prototypes on the basis of their similarity. This mapping $Q \in \mathbb{R}_+^{K \times B}$ is the solution of the following optimization problem:

$$ \arg \max_{Q \in \mathcal{Q}} \mathrm{Tr}(Q^T C^T Z) + \varepsilon H(Q) \qquad (1.31) $$

with $H(Q) = \sum_{i,j} Q_{ij} \log Q_{ij}$, $\varepsilon$ taking a low value ($\varepsilon \approx 0.001$), and

$$ Q = \left\{ Q \in \mathbb{R}_+^{K \times B}, \ Q1_B = \frac{1}{K}1_K, \ Q^T1_K = \frac{1}{B}1_B \right\}. \qquad (1.32) $$

The constraints in (1.32) prevent solutions where all images are mapped to the same prototype, by ensuring that the distribution of feature vectors over the different prototypes is roughly uniform.

The rest of the paper uses conventional ideas in unsupervised learning. We have two augmented views $x_s, x_t$ of the same batch of images. A ConvNet $f_\theta$ builds two sets of features $z_s = f_\theta(x_s), z_t = f_\theta(z_s)$ from those augmented views. For each view, those features
<figure><img src="image_20.png" /><figcaption>Figure 1.11 – Swapping Assignments between Views (SWaV). (Source: [53])</figcaption></figure>

are mapped to prototypes using the Sinkhorn-Knopp algorithm, yielding sets of codes $q_s, q_t$. Finally, the ConvNet $f_\theta$ is trained with the following criterion:

$$ \mathcal{L}(z_s, z_t) = - \left( \sum_k q_s^k \log p_t^k + \sum_k q_t^k \log p_s^k \right) \qquad (1.33) $$

with

$$ p_s^k = \frac{\exp\left(\frac{1}{\tau} z_s^T c_k\right)}{\sum_{k'} \exp\left(\frac{1}{\tau} z_s^T c_{k'}\right)}. \qquad (1.34) $$

The model is trained by predicting the prototype assignments of one view from the features of the other views. This ensures that $f_\theta$ produces features that are invariant to the augmentations used to produce the views $x_s, x_t$. This process is summed up in Figure 1.11.

## Similarity-maximization approaches

In [55], Chen et al. point out that a number of successful unsupervised learning architectures that avoid representational collapse have been designed (see Sections 1.2.2, 1.2.2). However, there has been little interest in the simplest architecture that can avoid this issue. For example, one could ask if negative samples are absolutely necessary to prevent representation collapse. Chen et al. suggest that this is not the case, and an architecture
<figure><img src="image_21.png" /><figcaption>Figure 1.12 – SimSiam architecture. Note that the dimension of z and p are identical.</figcaption></figure>

that relies only on maximizing similarity between positive samples is feasible. SimSiam [55] relies on a encoder network $f_{\theta}(x) : [0, 1]^{H \times W \times C} \rightarrow \mathbb{R}^d$ and a projector $p_{\theta}(z) : \mathbb{R}^d \rightarrow \mathbb{R}^d$, as shown in Figure 1.12. This neural network is trained using the following criterion:

$$ \mathcal{L} = \frac{1}{2} D(p_1, \operatorname{sg}(z_2)) + \frac{1}{2} D(p_2, \operatorname{sg}(z_1)), \qquad (1.35) $$

where sg denotes the stop-gradient operator. The similarity criterion is cosine similarity, i.e.,

$$ D(u, v) = \frac{u^T v}{\|u\|_2 \|v\|_2}. \qquad (1.36) $$

The key ingredient to prevent representation collapse is the stop-gradient operator. This means although $z = f_{\theta}(x)$, no gradient is computed with respect to $\theta$ for this term of the loss criterion. Without this stop-gradient operator, training quickly collapses to a degenerate solution where the same representation is predicted for every image. This technique does not rely on sampling negative examples like in [49], and there is no need to maintain a key encoder like in [48].

## 1.3 Approximate Bayesian deep learning

Because DNNs are often perceived to be "black boxes", there has been particular interests in modelling the sensitivity of the model output with respect to its parameters, and Bayesian approaches are a natural fit for this task. However, computing a posterior distribution over model parameters raises particular challenges, as modern DNNs have
been scaled to hundred of billions of parameters.

Bayesian deep learning techniques using tractable variants of Monte Carlo sampling have been designed, but those methods are usually slower than gradient-based methods. Furthermore, they require an accurate implementation of a completely different optimization procedure than what is usually used in modern DNNs. This complexity and computational cost has made them a relatively niche research direction [56].

Because of this, this section focuses on approximate Bayesian deep learning. First, we briefly describe some theoretical work describing approximate Bayesian deep learning based on SGD iterates. Then, we review some articles that use a set of models to improve performance or provide uncertainty estimation, either by combining their predictions or their weights. Finally, we review a several papers that focus on estimating approximate posterior distributions over the weights of a neural network. The posterior weight distribution is assumed to be a high-dimensional Gaussian. Several ways of fitting the mean and variance of this Gaussian distribution are proposed, either based on gradient descent iterates, or estimates of the Hessian matrix of the network.

### 1.3.1 SGD as an Orstein-Uhlenbeck process

In a well-known paper, Mandt et al. [7] used tools of the stochastic differential equations (SDEs) literature to study the stationary distribution of SGD. Let us consider a loss function over a dataset of size $N$ that depends on model parameters $\theta$:

$$ \mathcal{L}(\theta) = \frac{1}{N} \sum_{i=1}^{N} l_i(\theta), \quad g(\theta) = \nabla_{\theta} \mathcal{L}(\theta), \qquad (1.37) $$

where each $l_i$ is the contribution to the total loss of a single sample, i.e., $l_i(\theta) = l_i(x, \theta)$. In the case of a MAP estimate, the objective is the sum of a log-likelihood and a prior, so that

$$ l_i(\theta) = \log p(x_i|\theta) - \frac{1}{N} \log p(\theta). \qquad (1.38) $$

For computational reasons, the full gradient $g(\theta)$ is usually not computed, but estimated on a random sample $\mathcal{J}$ of the full dataset, i.e

$$ \hat{g}_S(\theta) = \nabla_{\theta} \mathcal{L}_{\mathcal{J}}(\theta), \qquad (1.39) $$
where $\mathcal{L}_{\mathcal{J}}(\theta) = \frac{1}{|\mathcal{J}|} \sum_{i=1}^{|\mathcal{J}|} l_i(\theta)$. The authors [7] make a set of assumptions to view SGD as an Orstein-Uhlenbeck process:

- As the batch is a random sample of the full dataset, the stochastic gradient (and therefore the gradient noise $g(\theta) - g_{\mathcal{J}}(\theta)$) is a sum of identically distributed independent contributions. According to the central limit theorem, it follows that

$$ \hat{g}_{\mathcal{J}}(\theta) \approx g(\theta) + \frac{1}{\sqrt{|\mathcal{J}|}} \nabla g(\theta), \quad \nabla g(\theta) \sim \mathcal{N}(0, C(\theta)), \qquad (1.40) $$

- the covariance matrix of the gradient noise is a symmetric positive-semidefinite matrix which is approximately constant with respect to $\theta$, i.e

$$ C(\theta) \approx C = BB^T, \qquad (1.41) $$

- the difference between two successive iterates is

$$ \Delta\theta(t) = -\varepsilon g(\theta(t)) + \frac{\varepsilon}{\sqrt{|\mathcal{J}|}} B \Delta W(t), \quad W \sim \mathcal{N}(0, I). \qquad (1.42) $$

Furthermore, the authors assume that this finite-difference equation can approximate the following continuous-time SDE:

$$ d\theta(t) = -\varepsilon g(\theta)dt + \frac{\varepsilon}{\sqrt{|\mathcal{J}|}} B dW(t), \qquad (1.43) $$

- the stationary distribution of SGD iterates is constrained to a region where the loss can be approximated in the following fashion:

$$ \mathcal{L}(\theta) \approx \frac{1}{2} \theta^T A \theta, \qquad (1.44) $$

where $A$ is positive definite matrix, the Hessian of the loss at the optimum.

Thanks to those assumption, SGD can be defined as an Orstein-Uhlenbeck process and the Gaussian stationary distribution $q(\theta)$ can be computed analytically:

$$ q(\theta) \propto \exp \left[ -\frac{1}{2} \theta^T \Sigma^{-1} \theta \right]. \qquad (1.45) $$

The authors note that those assumptions do not hold in the general case for DNNs. During early phases of training, the iterates are not located in the vicinity local mininum. It is
only where the iterates are "close" to a local minima that the assumptions used above prove correct. If we use a very long training schedule with large enough gradient steps, several local minima might be explored during training.

### 1.3.2 Ensembling neural networks for performance and uncertainty estimation

In this subsection, we focus on papers [57]–[59] that use the intuitions detailed above, even if they do not explicitly estimate the parameters of the approximate posterior distribution. The authors of those papers combine several models, either by combining their predictions or their weights.

In [57], Lakshminarayanan et al. evaluate the performance of ensembles of neural networks for two specific goals: higher performance, and better model calibration. Calibration describes the relationship between model confidence and error rate. DNNs are usually overconfident: their predicted probability does not match their error rate. In other words, a classifier predicting a class with a probability of 0.99 does not achieve a 1% error rate on this class.

Models in the ensemble are trained in a conventional way, aside from the use of adversarial training [40]. For a set of weights $\theta$ and an image-label pair $(x, y)$, an additional "adversarial" training example $x'$ can be generated as

$$x' = x + \varepsilon\mathrm{sign}(\nabla_x \mathcal{L}(\theta, x, y)), \qquad (1.46)$$

Where $\varepsilon$ is a small random value. This additional training example has the same label $y$ as the original one but is displaced along the gradient of the loss with respect to the original input $x$. As the model has to predict the same label for the original training sample and the adversarial example, this technique regularizes training and makes the loss surface smoother and the classifier more robust.

In [57] The authors note that conventional models usually need additional randomness to ensure that the ensemble has enough diversity to provide an increase in performance. This randomness can be injected by training the model on random samples of the data (bagging), or by choosing a base learner with a stochastic training procedure (like random selection of features for random forest). However, the batch sampling of DNNs usually
<figure><img src="image_22.png" /><figcaption>Figure 1: Left: Illustration of SGD optimization with a typical learning rate schedule. The model converges to a minimum at the end of training. Right: Illustration of Snapshot Ensembling. The model undergoes several learning rate annealing cycles, converging to and escaping from multiple local minima. We take a snapshot at each minimum for test-time ensembling.</figcaption></figure>

Figure 1.13 – In [58], parameter space is explored using a cyclical learning rate schedule. An ensemble is sampled from saved checkpoints.

provides enough randomness. Models are aggregated as an uniform mixture:

$$p(y|x) = \frac{1}{M} \sum_{m=1}^{M} p(y|x, \theta_m), \qquad (1.47)$$

where *M* is the number of models in the ensemble. For classification, this aggregation method provides an improvement in performance, and improved model calibration. While the authors demonstrate the potential improvements brought by ensembling DNNs, training *M* different models remain costly. In [58], the authors find that a specific training schedule can be used to extract a diverse enough ensembles from checkpoint of a single training run. To ensure a thorough exploration of the parameter space, the authors use a cyclic learning rate schedule. When the learning rate is low, the model converges to a local minima, and the modelling suggested in Section 1.3.1 is correct. When the learning rate gets higher, gradient steps are larger and the model parameters escape the current local minima. Figure 1.13 illustrates this phenomenon.

While this method is computationally cheaper than training *M* models, the period of the cyclical learning rate remains high (around 10 epochs), so the training cost is increased. In [59], Garipov et al. suggest an alternative. Noticing that most local minima are connected by curves where the loss stays low, they suggest retrieving model weights
<figure><img src="image_23.png" /><figcaption>Figure 1.14 – $L_2$-regularized cross-entropy train loss (left) and test error (middle) as a function of the interpolation parameter $t$ between local minima. Straight lines between local minima showcase a large increase in loss, while a Bezier interpolation between both local minima shows the presence of a low-loss path between those models. "Polychain" (a polygonal chain between minima) shows similar behavior.</figcaption></figure>

along those curves. Ensembles can be built by averaging the predictions of the model retrieved along curves linking a pair of local minima. The low loss path along this curve is shown in Figure 1.14.

The authors also confirm the insight of [58], but show that a strong ensemble can be retrieved with a cyclical learning that has a much shorter period (2 to 4 epochs), with improved performance.

All methods mentioned previously perform model aggregation in prediction space. However, models can also be aggregated in weight space. This is what Izmailov et al. demonstrate in [13]. While using an Exponential Moving Average (EMA) of SGD iterates is fairly common to stabilize training, this paper uses a uniform average of SGD iterates, called Stochastic Weight Averaging (SWA). An early example of uniform weight averaging for regularization is Polyak-Ruppert averaging [60], [61]. In this paper, the strength of this idea is demonstrated even in the context of non-convex optimization for DNNs. The main difference between SWA and Polyak-Ruppert averaging is that SWA uses a constant high learning rate instead of a decaying learning rate, to ensure parameter space exploration. The authors hypothesize that because once a local minima is reached, SGD samples the posterior weight distribution (see Section 1.3.1), the barycenter of those weights is more likely to lie near the "center" of the local minima. In this location, it is less likely that a small step in weight space will lead to a large increase in loss. In other words, this minima is "flatter", which has long been hypothesized to correlate with improved generalization [62], [63].

While all of those papers use the insights provided by Mandt et al. in [7] and the literature dealing with the limit distribution of SGD in non-convex optimization in general,
none of them actually estimate a posterior distribution over the weights of the model. We cover a few references dealing with this problem in the next section.

## Approximating a posterior distribution with SGD iterates

In [6], Maddox et al. approximate the posterior distribution over the model weights $\theta$ as a Gaussian distribution, and call this method Stochastic Weight Average Gaussian (SWAG). The mean of this distribution is estimated with the running mean of SGD iterates, i.e the SWA solution $\theta_{SWA}$ as described in [13]:

$$ \theta_{SWA} = \frac{1}{T} \sum_{i=1}^{T} \theta_i. \qquad (1.48) $$

Several covariance approximations are possible, although computational costs are an issue, as the size of the covariance matrix scales quadratically with the number of parameters. Because of this, a diagonal covariance can be used:

$$ \Sigma_{\text{diag}} = \text{diag}(\bar{\theta}^2 - \theta_{\text{SWA}}^2), \qquad (1.49) $$

with $\bar{\theta}^2 = \frac{1}{T} \sum_{i=1}^{T} \theta_i^2$. This diagonal covariance model can be enriched with an additional low-rank covariance matrix. This approximation has the following form:

$$ \Sigma_{lr} = \frac{1}{K-1} \sum_{i=1}^{K} (\theta_i - \theta_{SWA})(\theta_i - \theta_{SWA})^T = \frac{1}{K-1} DD^T. \qquad (1.50) $$

Note that as $\theta_{SWA}$ is not available until the end of training, it is approximated during training by $\bar{\theta}_i$, which is the running mean of the first $i$ SGD iterates. $D$ is a deviation matrix with $K$ columns of value $\theta_i - \bar{\theta}_i$. The number of columns $K$ of $D$ (the number of considered SGD steps) is an hyperparameter. Both covariance approximations can be combined, so that the final posterior approximation is $N(\theta_{SWA}, \frac{1}{2}(\Sigma_{\text{diag}} + \Sigma_{lr}))$. One can sample from this approximation with the following equation:

$$ \tilde{\theta} = \theta_{SWA} + \frac{1}{\sqrt{2}} \sum_{\text{diag}}^{1/2} z_1 + \frac{1}{\sqrt{2(K-1)}} D z_2, \qquad (1.51) $$

where $z_1 \sim N(0, I_d)$ and $z_2 \sim N(0, I_K)$ where $d$ is the number of parameters in the network. This approximate posterior can be used to approximate a Bayesian Model Average (BMA) [64]:
<figure><img src="image_24.png" /><figcaption>Figure 1.15 – Comparison of the sampling of the posterior distribution for several methods. An ensemble of several DNNs trained separately [57] effectively provides several true samples of the posterior distribution. However, the number of samples is low. SWAG can cover a local minima of the loss function (a mode of the posterior distribution) much more densely, but can only cover one mode. Multi-SWAG covers several modes of the posterior distribution by fitting several Gaussian distributions during training, to form a mixture model of the posterior distribution.</figcaption></figure>

$$p(y_t|\mathcal{D}, x_t) = \int p(y_t|\theta, x_t)p(\theta|\mathcal{D})d\theta, \qquad (1.52)$$

where $\mathcal{D}$ is the data distribution, $(x_t, y_t)$ is a test input-output pair and $p(\theta|\mathcal{D})$ is the posterior weight distribution. This integral is approximated with samples of the approximate posterior distribution estimated above:

$$p(y_t|\mathcal{D}, x_t) \approx \frac{1}{S} \sum_{s=1}^{S} p(y_t|\theta_s, x_t), \quad \theta_s \sim p(\theta|\mathcal{D}). \qquad (1.53)$$

Ensembles generated from the posterior weight distribution bring improvements in terms of performance and model calibration, and do so at a lower compute cost than ensembles composed of separately trained models. However, true posteriors of neural networks are very unlikely to be unimodal. The stationary distribution proposed in Mandt et al. [7] is only relevant in the vicinity of a local minima. Therefore, while the approximation proposed above provides performance improvements, it does not model this multimodality. In Maddox et al. [6], the authors improve SWAG by representing the posterior distribution
as a mixture of Gaussians, centered in several different basins of attraction of the loss surface. Figure 1.15 provides a schematic summary of the idea.

## 1.4 Object detection with deep learning models

Object detection is most commonly defined as fitting axis-aligned bounding boxes around objects of interest. This task is relevant in a number of fields, like autonomous driving or medical imaging. Bounding boxes have an associated class label, so that any object of interest is localized and classified. In two-stage object detectors, the detection task is divided into three subtasks. First, there is a region proposal step, where a model predicts bounding boxes around foreground objects. Then, a bounding box regression model refines those proposals, using the content of the region proposal. Finally, a classification networks predicts a class label for the bounding box.

One-stage detectors do not have any region-proposal. This is done to improve inference time, but usually lowers detection performance. In this section, we will review two object detection models and explain how they solved the most common problems associated with this task. First, we look at Faster R-CNN as a well-known example of a two-stage object detector. Then, we study CenterNet as an example of a one-stage object detector, as it is very simple to implement which makes it especially interesting.

### 1.4.1 Object detection with Region-CNN

#### Fast R-CNN

In "Fast R-CNN", Girshick et al. [65] describe the basic components of their object detection model. Figure 1.16 gives a summary of the architecture. Fast R-CNN region proposal step is separated from the object detection model. At the time, Selective Search was used for region proposal [66]. Therefore, the complete architecture could not be trained end-to-end.

For a given image, a feature map is produced by a ConvNet which is usually called a backbone. Given that region proposals might have different size, the subsets of the feature map corresponding to those region proposals will also have different sizes. To ensure that features corresponding to those region proposals all have the same size, a Region of Interest (ROI) pooling layer is used. This adaptive pooling layer aggregates a feature map to a constant output size, no matter the size of the input. Once this is
<figure><img src="image_25.png" /><figcaption>Figure 1.16 – For a proposed region, features from a ConvNet are retrieved. An adaptive pooling mechanism called ROI pooling ensures that features have the same spatial size, no matter the size of the proposed region. Pooled features are used for classification and bounding box regression. (Source: [65])</figcaption></figure>

done, two small neural networks are used for bounding box regression and bounding box classification. The object detector is trained with the following criterion:

$$L(p, u, t^u, v) = L_{cls}(p, u) + \lambda L_{loc}(t^u, v) \quad (1.54)$$

with $p$ the predicted probability distribution over the classes, and $u$ the one-hot vector describing the true class label. $L_{cls}$ is a cross entropy loss. $t^u = (t_x^u, t_y^u, t_h^u, t_w^u)$ is the predicted bounding box regression vector and $v$ the true bounding box vector. The first two coordinates of the bounding box vector are the coordinates of the top right corner. The last two coordinates indicate the height and width of the bounding box. The bounding box regression loss is defined as the following robust loss:

$$L_{loc}(t^u, v) = \sum_{i \in \{x,y,h,w\}} \mathrm{smooth}_{L1}(t_i^u - v_i) \quad (1.55)$$

with smooth$_{L1}$ defined as:

$$\mathrm{smooth}_{L1}(r) = \begin{cases} \frac{1}{2}r^2 & |r| < 1 \\ |r| - \frac{1}{2} & |r| > 1. \end{cases} \quad (1.56)$$

Note that this bounding box regression loss is only active if the considered region does
<figure><img src="image_26.png" /><figcaption>Figure 1.17 – For every pixel of the feature map produced by the backbone, the RPN produces a set of regression vector modifying the pre-defined anchors, and a set of objectness scores which quantify how likely each anchor is to contain an object. (Source: [9])</figcaption></figure>

not belong to the background. During training, background examples (regions with a small intersection with a ground truth bounding box) are sampled to ensure that the classification network accurately classifies the background.

## Faster R-CNN

In [9] Ren et al. suggested training a model for regression proposal, called a Region Proposal Network (RPN) so that the entire object detection model could be trained end-to-end. This network uses the same backbone as the object detection model (see Section 1.4.1 for an explanation). Figure 1.18 shows the architecture of Faster R-CNN. For every pixel of the feature map, there are *k* (where *k* is a tunable hyperparameter) anchor boxes, usually with different aspect ratios or size, as shown in Figure 1.17. This is mostly because the resolution of the output feature map is relatively low, so that a single anchor per pixel would be insufficient to accurately locate every object in the image. Those anchors act as pre-defined region proposals. The RPN refines those anchors by predicting a regression vector (like the object detection model). Furthermore, the RPN predicts an objectness
<figure><img src="image_27.png" /><figcaption>Figure 1.18 – Architecture of Faster R-CNN. A ConvNet produces a feature map for an input image. This feature map is used to produce region proposals. Those region proposals are selected, and an object detection network is trained as presented in Figure 1.16. (Source: [9])</figcaption></figure>

score for every anchor. This score corresponds to the likelihood that the current bounding box contains an object, instead of background.

The RPN is trained like the object detection model, by minimizing the sum of a classification and regression loss. Here, the classification task is the rejection of background. While this idea makes Faster R-CNN considerably faster than Fast R-CNN, there an issue. For 2 adjacent pixels of the feature map at the output of the backbone, 2 anchors can be adjusted to cover the exact same subset of the image. In other words, several regions with a high objectness score can be proposed for a single object.

To deal with this issue, the set of region proposals is reduced by thresholding the objectness score. Furthermore, overlapping boxes in the set of object proposals are merged using Non-Maximum Suppression (NMS). If two boxes have high overlap (defined by their Intersection Over Union, or IoU), only the one with the highest objectness score is kept. This method is effective, but requires setting another hyperparameter, and the enumer-
<figure><img src="image_28.png" /><figcaption>Figure 1.19 – Output of the CenterNet model. The center of the bounding box is predicted, a specific heatmap is predicted for every class in the dataset. For every location of the heatmap, a (height, width) tuple is predicted, representing the size of the bounding box at this location. (Source: [67])</figcaption></figure>

ation of region proposals in two-stage detectors is computationally costly. However, this class of object detection models usually reaches a high level of performance. In the next section, we discuss a one-stage object detector that does not require region proposal or NMS.

1.4.2 NMS-free object detection: Object as Points

In [67] Zhou *et al.* proposed a simple one-stage object detector that does not need NMS
nor any pre-defined anchors. Instead, the object detection model (named CenterNet) is
implemented as a ConvNet with three output branches. For an image of size *H* × *W*, a
first branch outputs a heatmap of size *H/K* × *W/K* × *C* (with *K* a downsampling factor,
and *C* the number of classes) predicting the center of bounding boxes. A second branch
predicts the size of the bounding boxes as a (height, width) tuple for every location *u*, *v*
in the output of the heatmap branch, with *u*, *v* ∈ {1, ..., *H/K*} × {1, ..., *W/K*}. Figure
1.19 gives a visual depiction of the output of CenterNet. Bounding box height and width
are predicted for every location of the image, even when this location does not contain
any object. Bounding box height and width are only retrieved where the output of the
bounding box center prediction branch exceeds a pre-defined threshold.

Finally, the heatmap output has resolution *H/K* × *W/K*, which is lower than the
resolution of the input image. Therefore, a third output branch is used to predict an
offset for every location *u*, *v* ∈ {1, ..., *H/K*} × {1, ..., *W/K*} so that the model can
predict bounding box centers with sub-pixel accuracy.
Training

The heatmap branch is trained by minimizing a modified cross entropy loss called Focal Loss [68], [69] between its output and a heatmap of Gaussian "spots" laid over the center of bounding boxes. While one could predict a binary image where bounding box centers would take value 1 and the background value 0, bounding box centers would be too under-represented. Therefore, Zhou *et al.* [67] choose to blur the heatmap with a Gaussian kernel.

Local maximas of the heatmap prediction are retrieved. Those location corresponds
to the predictions of the model. Bounding box size and offset are retrieved in the outputs
of the two other branches. Bounding box and offset errors are computed only in those
locations. Both branches are trained by minimizing a L1 distance between prediction and
ground truth.

1.5 Conclusion

Modern deep learning has received considerable attention in image analysis. Extremely flexible architecture can be trained end-to-end with simple criterions, using backpropagation and stochastic gradient descent. In spite of the non-convexity of commonly used training criterions because of the non-linearities in this class of models, finding "good" local minimas is relatively easy because of implicit regularization.

Bringing the benefits of this class of high-performing architectures to biomedical imaging is more difficult, as labelled data is rare. SSL lets us pre-train models without any labels, either by learning to generate images or by learning invariance to common image augmentation techniques (noise, cropping, flipping, etc).

Because of the black box nature of modern DNNs, Bayesian deep learning has been an active field of research. In approximate Bayesian deep learning, the posterior distribution over the weights of the model is usually modelled as Gaussian. The parameters of this distribution can be estimated with approximate samples, like SGD iterates.

Finally, DNNs can be used to build object detection models. Those models are usually divided in two categories: one stage and two stage detectors. Two stage detectors are more complex, but reach a higher level of performance. In Chapter 4, we explore the potential of one-stage object detection models for chromosomal aberration detection.
# IMAGE ANALYSIS FOR THE AUTOMATION
OF BIOLOGICAL DOSIMETRY

In the first part of this chapter we introduce the basics cytogenetic biological dosimetry. We review some biological factors that are especially relevant from a computer vision perspective as they cause variations in image quality. Next, we review historical examples of chromosomal aberration detection methods. We focus on a commercial solution developed by Canadian researchers as its inner workings were explained in details in a series of publications.

Finally, we review some references written by researchers at LRAcc dealing with evaluating another commercial solution in a semi-automated regime. While the computer vision techniques used in this solution were not described in the open scientific literature its evaluation by biological dosimetry specialists provides a benchmark to compare against for the new methods developed in this thesis.

## 2.1 Basics of biological dosimetry

The goal of biological dosimetry is estimating a dose of ionizing radiation using biomarkers instead of measuring it with instruments like a dosimeter [14]. This may be necessary because dosimeter readings are unavailable, in the case of an accidental (or potentially malicious) exposition. While several biomarkers of radiation have been identified, the average number of DC per peripheral blood lymphocytes (also called dicentric yield) is the gold standard according to the IAEA [14].

Those aberrations can be observed during the metaphase step of mitosis, so that an image of a peripheral blood lymphocyte undergoing metaphase is usually abbreviated to "metaphase". Once a blood sample is collected, blood cells are grown and Demelcocine is used to stop mitosis during metaphase. The blood sample is spread on a transparent slide and observed with a bright-field microscope. The slide is backlit, and because chromo-
<figure><img src="image_29.png" /><figcaption>Figure 2.1 – Basic morphology of a metaphase chromosome.</figcaption></figure>

somes are translucent objects, Giemsa staining is used to make them opaque and improve
contrast with the bright background [70], [71] (see Figure 2.2).

DCs have two centromeres, which is the visible constriction of its chromatides (like at the center of a MC). For an example of a DC, see the shape circled in red in Figure 2.2 b). The presence of DCs and fragments in a metaphase (circled in green in Figure 2.2 b) ) are highly correlated. The required cell culture duration for DC scoring is lower than for other methods like micronuclei scoring [70]. Furthermore, the preparation process is simple and relatively cheap. Other imaging modalities, like Fluorescence In Situ Hybridization (FISH) offer significant advantages: probes can be used to highlight other types of aberrations like translocations, or centromeres [72], [73]. However, FISH is significantly more time consuming, and more expensive. MCs and DCs can be discriminated using their morphological differences.MCs are X-shaped (see Figure 2.1), while DCs look like a pair of MCs bound together.

The number of DCs in a metaphase is random. It follows a Poisson distribution, where the parameter λ of this distribution is the average number of DC per cell (also called dicentric yield) for the considered dose. This Poisson model is adequate in the case of a whole body photon exposition. The relationship between DC yield and dose (in Grays) is described by a calibration curve with a linear-quadratic shape. In the case of a partial exposition, only some blood cells are exposed to ionizing radiation. The blood sample is usually modelled as a mixture of an irradiated and a non-irradiated sample [74]. This effectively increases the number of cells without any DC, and partial expositions are often modelled with Zero Inflated Poisson models. In this context, the DC distribution
<figure><img src="image_30.png" /><figcaption>Figure 2.2 – Chronology of a biological dosimetry-based estimation of an ionizing radiation dose. First, a blood sample is taken and cells are grown for 48 hours, Demelcocine is used to stop cells in the metaphase stage of mitosis. Metaphase acquisition is automated with modern microscopy systems, as shown in a), and chromosomal aberrations are counted, as displayed in b). The proportion of chromosomal aberration is linked to the ionizing radiation dose through a linear-quadratic relationship, see c).</figcaption></figure>

showcases overdispersion: the variance of the distribution is greater than the mean.

For a blood sample that was not exposed to ionizing radiation, the dicentric yield is around $10^{-3}$: the vast majority of cell does not contain any aberration. For small doses, a large number of metaphases has to be processed to accurately detect a deviation from the basal aberration rate. This leads to stringent requirements in terms of false positives, to prevent overestimation of low doses. Even at 5Gy, which corresponds to the median lethal dose for a whole body exposure (also called $LD_{50}$), there is only around 1 aberration per cell. This means that in a single metaphase, on average, 45 chromosomes are healthy and a single DC is observed.

To estimate dicentric yield, a sample of metaphase is checked for chromosomal aberrations, a process called DC scoring. In practice, human experts perform binary classification on each chromosome, sorting them between the MC and DC class. Furthermore, multicentrics chromosomes and fragments are identified. To this end, those experts undergo several months of training and are then evaluated on benchmark datasets. This training helps humans achieve an extremely low error rate. However, this error rate comes at a cost: human can process a relatively low number of metaphases per hour (see Figure 2.3). Cognitive load is high, and human experts can only count chromosomal aberrations for a few hours at a time to keep the error rate low. Furthermore, a large number of metaphase is needed to achieve usable confidence intervals for low doses, as illustrated in Figure 2.3.

Within the context of a large scale exposure there may not be enough human experts to accurately triage patients between different therapeutical options. In this case, com-
<figure><img src="image_31.png" /><figcaption>Figure 2.3 – Relative confidence interval size as a function of metaphase number and dose. The lower the dose, the higher number of images is needed to get a narrow confidence interval. Aberration counting takes about one hour for 50 images, and a thousand image may be needed for a narrow enough confidence interval for a small dose.</figcaption></figure>

puterized expertise could provide invaluable help. This is the central motivatoin for the
deep-learning-based ADS proof of concept presented in this thesis.

**2.2 Challenges of automated aberration counting**

In Figure 2.2 a), the DC circled in red is relatively obvious even to the untrained
eye. Unfortunately, chromosomes are flexible objects and the chromatids of a MC may
fold in way that would make a MC indistinguisable from a DC. Human experts can shift
the focal plane of the microscope to disambiguate those cases, as folded chromatids and
<figure><img src="image_32.png" /><figcaption>Figure 2.4 – Four examples of metaphases showcasing variations in acquisition conditions. Variations in chromosome spread, chromosome condensation, chromosome thickness, chromosome texture and illumination are visible. Furthermore, some metaphases contain more debris (cellular nuclei, filaments) than other.</figcaption></figure>

centromeres have different thickness. This can then be used to distinguish MCs from DCs, as DCs have two centromeres. Because the microscope-mounted camera only acquires 2D data, this information is not available in the data we are working with.

Beyond those issues, the process of cell culture and metaphase acquisition leads to variation in image quality. In each cell, mitosis is stopped at a different level of chromosome condensation, which leads to difference in chromosome size or texture. Giemsa staining may bind to the chromosomes in a non-uniform manner. The location of the cell on the
glass slide leads to variations in backlight intensities between metaphases, or even in a
single metaphase, so that a gradient in backlight intensity could be observed even for
a single cell. A lot of metaphases contain cellular debris, which might be circular and
easy to reject but could also be filament-shaped and much harder to discriminate from
a chromosome fragment. In Figure 2.4 a few examples of metaphases showcasing those
variability factors in the image acquisition process are shown.

Any computer vision system aiming to accurately count chromosomal aberrations in
those conditions will need to deal with those sources of variability, either in an explicit or
implicit manner.

## 2.3 An historic example of automated dicentric scoring systems

While biological dosimetry is a relatively niche field of study, it quickly became of
playground for advances in computer vision and image processing. One of the earliest
implementation of an ADS was called the Edinburgh dicentric hunter. It was introduced
in "Radiation dosimetry by automatic image analysis of dicentric chromosomes" by Bayley
et al. [70] in 1991. Because of the lack of sufficiently effective computer vision tools, the
authors choose to significantly restrict the range of DC they would try to detect. For
example, they do not attempt to detect DCs where the centromeres are very close to each
other, or DCs where one of the centromeres is located at the edge of the chromosome
(acrocentric chromosomes). Based on previous references, the authors assume that those
configurations amount to around 30% of all DCs.

The authors also note that one of the significant challenges in evaluating an ADS is
building a ground truth. Indeed, there is usually a high level of inter-observer disagreement
during a manual scoring session of DCs. In [70], DC scoring is divided in 5 steps:

- metaphase finding;

- image acquisition;

- segmentation of chromosomes;

- MC - DC classification;

- human review.

The metaphase finding step was implemented with the Cytoscan 110 metaphase finder,
which was evaluated in [76]. During metaphase acquisition, the images undergo global
thresholding. The software attempts to identify chromosome objects, reject debris and
<figure><img src="image_33.png" /><figcaption>Figure 2.5 – This figure depicts various degrees of SCS present in some Giemsa stained chromosome cell images (a), b) and c)) as well as some lenghty chromosomes characteristic to those prepared at a cytogenetic laboratory (d), e) and f)). (Source: [75])</figcaption></figure>

split chromosome clusters. Metaphases are accepted or rejected on the basis of their object counts. Early and late metaphase cells are rejected by identifying the level of chromosome condensation based on their mean aspect ratio. This prevents the selection of metaphases containing chromosomes with excessive Sister Chromatide Separation (SCS). Sister chromatides are the "branches" of a chromosome. Depending on the preparation of the blood sample, various levels of SCS might be encountered, as shown in Figure 2.5.

Once metaphase finding is completed, chromosome classification is performed. The classifier used by [70] is described in [77]. A first-stage classifier extracts a set of centromere candidates by retrieving points of contour constriction along the chromosome centerline. Because of this, SCS can lead to spurious DC detections as constrictions points might not correspond to centromeres. As the authors point out, chromosome shape variability is usually much higher between cells than inside a single cell. Therefore, the centromere clas-
sifier must be adapted to every cell. This can be achieved by finding a within-chromosome
or within-cell normalization for feature values. A feature set is extracted at candidate
points, as illustrated by Figure 2.6.

<figure><img src="image_34.png" /><figcaption>Figure 2.6 – Chromosomes showing the various structures from which feature values are computed: (a) the chromosome and its computed symmetry axis, with candidate centromere positions marked by horizontal lines, (b) the longitudinal profile of densities, (c) the chromosome boundary from which curvature values are obtained, (d) the "crossing profiles", the density profiles across the chromosome at the centromere candidates. The right-hand chromosome has a false candidate, clearly distinguished by the bimodal crossing profile. (Source: [70])</figcaption></figure>

The set of centromere candidates is re-classified by a second-stage classifier which takes
this feature set as an input. The set of detected DCs is stored, and presented to a human
expert for review. The software retains the position of DCs in every cell, which makes
human review with a microscope (instead of on an image acquired by a camera) possible.
In practice, this possibility was found to increase analysis time, as human experts tended
to prefer lengthy reviews using the microscope over fast review using DCs candidates
extracted by the software.

In terms of performance, The Edinburgh dicentric hunter achieves a false positive rate
of around 0.5 DC per cell according to Bayley *et al*. The system has a recall of around
40%. The sensitivity for low doses is insufficient, which is not surprising given that no
attempt is made to detect around 30% of DCs. The results of Bayley *et al*. [70] were
confirmed by another evaluation of the same software in Finnon *et al*. [78].

To summarize, ADS for Giemsa stained images data back to the mid 1980s. The Ed-
inburgh dicentric hunter is an noteworthy example of early ADS dating back to 1991.
The overall workflow [70] by Bayley *et al*. is still used by the modern method presented
in Section 2.4. Other ADS systems have been designed with FISH imaging, which en-
ables labeling centromeres or specific chromosomes. For example, Roy *et al.* [71] made a comparison of ADS which includes fluorescence-based systems and Huber *et al.* [73] did technical review of a version of Metafer for fluorescence imagery.

## 2.4 ADCI: implementing an ADS framework

This section aims to provide an overview of a commercially available ADS solution developed by Cytognomix called Automated Dicentric Chromosome Identifier (ADCI). The first version of this dose estimation pipeline was described in "Towards large scale automated interpretation of cytogenetic biodosimetry data" by Li *et al.* [79] in 2012. While there are other commercially available ADS, we have more insight in the inner workings of ADCI through a series of paper providing details of its implementation.

ADCI does not use deep-learning-based object detection. DC counting is framed as two successive tasks: chromosome detection, and chromosome classification. A large amount of prior knowledge is used to solve those tasks. Metaphases are selected based on the number of objects, so that images with overlapping objects are rejected, chromosome contours are identified accurately for precise centerline extraction, centromeres are detected for accurate DC classification, etc.

### 2.4.1 Chromosome feature extraction

First, an accurate chromosome extraction algorithm is needed. Noise may lead to incorrect chromosome extraction when global thresholding algorithms (like Otsu thresholding [80]) are used, like in [70]. Accurate segmentation should rely on more than a single pixel for its decision. In [81], the authors found that an active contour method called Gradient Vector Flow (GVF) provided good segmentation results. Active contours ensure that the chromosome segmentation method has a relatively regular boundary.

Once an accurate chromosome contour (and segmentation) has been extracted, the centerline of the chromosome must be estimated. Usually, skeletonization is performed using morphological operators (e.g. thinning).

This step emphasizes the importance of an accurate contour extraction, as the estimated centerline location depends on the width of the estimated contour at every location of the chromosome boundary. Furthermore, biological processes can make centerline extraction difficult, the most likely being Sister Chromatide Separation (SCS) in late
<figure><img src="image_35.png" /><figcaption>Figure 2.7 – Chromosome centerlines proposed by the GVF + DCE algorithm proposed in [81]</figcaption></figure>

metaphase images. A further refinement of centerline extraction was introduced in [82].

Spurious branches may be produced during centerline extraction, especially if the chromosome is bent. In those cases, the authors performed skeleton pruning with a Discrete Curve Evolution (DCE) algorithm which was proposed in [83]. As explained in [81], this centerline extraction algorithm relies on the adjustment of hyperparameters for the snake algorithm. Those were set empirically. Figure 2.7 shows a few examples of extracted centerlines.

This centerline can be used to estimate centromere location and number, which is a common technique in biological dosimetry (see Section 2.3). This is done by building a width profile of the chromosome around the centerline by retrieving the length of scan lines that are perpendicular to the centerline. Identifying constrictions in this width profile provides a way to count centromeres, as shown in Figure 2.8.

### 2.4.2 Chromosome classification

In the first ADCI prototype introduced by Li *et al.* in 2012 [79], monocentric-dicentric (MC-DC) classification relied on counting the minimas of the chromosome width profile, like in early ADS systems like the Edinburgh dicentric hunter introduced by Bayley *et al.* in 1991 [70].

The MC-DC classification was improved in Subasinghe *et al.* and Li *et al.* in 2016 [75], [85]. For centromere detection, a contour is extracted using GVF and DCE, and the
<figure><img src="image_36.png" /><figcaption>Figure 2.8 – Chromosome images processed by ADCI, annotated with key segmentation features. (A) MC and (B) DC. Chromosome contour overlaid in green, long-axis centreline in red. For reference, the minimum bounding box of the contour is also displayed in magenta and green. Yellow and cyan markers on the centerline indicate the top-ranked and 2nd-ranked centromere candidates, respectively, and all other candidates are indicated with a dark blue marker. For each centromere candidate, their corresponding width traceline (crossing through the candidate and running approximately orthogonal to the centerline) are displayed in dark blue. The arc lengths of width tracelines running down the centerline (not all shown) are used to construct a chromosomal width profile. Note that for the MC (A), the top-ranked candidate correctly labels the true centromere location, while the 2nd-ranked candidate labels a minor non-centromeric constriction. Meanwhile, for the DC example (B), both the top and 2nd-ranked candidates label true centromere locations. By comparing features extracted from the top 2 candidates (including width and pixel intensity information), the software determines if the chromosome is a MC or a DC. (Source: [84])</figcaption></figure>

centerline is extracted, as introduced in [82]. Local minima of the width profile of the chromosome along this centerline are considered to be candidate points for centromeres. A set of handcrafted features is built from those locations, like local curvature of the contour, thickness of the chromosome at this point, distance from the candidate to the closest end of the centerline, etc. Those features are used by a SVM classifier to accept or reject candidate points as centromeres. The distance to the separating hyperplane is used as confidence measure, called Candidate Based Centromere Confidence (CBCC). For every chromosome two candidate points are extracted, as shown in Figure 2.8. Additional steps were taken by Liu *et al.* [84] to reduce the number of false positive in MC-DC classification. A set of filters was designed to correct SVM classifications. For example, if the surface area of a candidate DC is too low, it should be rejected. Filters based on chromosome maximum and median width, surface area, oblongness were designed.
<figure><img src="image_37.png" /><figcaption>Figure 2.9 – Outputs of the 2014 version of ADCI at different steps of the pipeline as shown in [86]</figcaption></figure>

### 2.4.3 A pipeline for dose estimation

Several iterations of the ADCI software have been published. All of them follow roughly the same outline, but ADCI was improved over time to deal with its initial shortcomings. Using the centerline extraction algorithm and the centromere detection method proposed by Subasinghe et al. in 2010 [81], [87], Li et al. [79] introduced the first prototype of ADCI and its parallelization over a compute cluster to deal with mass exposition in 2012. Particular care was taken to ensure that the results are as independent as possible of preparative methods in biodosimetry labs, as this is identified as a core factor of ADS performance variation. Using the centerline extraction improvements proposed in [82], the first iteration of the ADCI software beyond the prototype stage was introduced by Rogan et al. in 2014 [86]. A flowchart of this pipeline is visible in Figure 2.9.

A new version of ADCI was introduced by Rogan *et al.* [88] in 2016. This version uses an improved centromere location algorithm and MC-DC classifier introduced by Subasinghe *et al.* and Li *et al.* in 2016 [75], [85]. This version of ADCI is compared to human evaluation in terms of calibration curve fitting. First, the authors note that the
<figure><img src="image_38.png" /><figcaption>Figure 2.10 – Calibration curves of the 2016 version of ADCI introduced in Rogan et al. [88]. The two human expertises are provided by HC, CNL (plotted as AECL). Three calibration curves are plotted for ADCI, depending on the hyperparameter $\sigma$ of the MC-DC classifier.</figcaption></figure>

number of segmented objects is a very reliable feature for metaphase selection: if the number is too low, it may be because its chromosomes are overlapping, therefore reducing the number of detected objects. On the other hand, if the number of objects is too high, the image may contain chromosomes from two cells instead of one. The parameter $\sigma$ of the improved SVM MC-DC classifier (which is the standard deviation of the radial basis function used as a kernel) introduced in [85] can be tuned by the user to balance specificity and sensitivity. If irradiation doses are known for a set of reference blood samples, ADCI has the ability to estimate a calibration curve. In [88], this is done in comparison to manual scoring provided by two labs, Health Canada (HC) and CNL (Canadian Nuclear Laboratories, formerly Atomic Energy of Canada Limited, or AECL). The calibration curve comparing ADCI to those laboratories is displayed in Figure 2.10.

Nevertheless, as shown in Figure 2.10, a large number of DC detections in this version of ADCI are false positives, which tends to skew calibration curves upward at the origin. While the aberration frequency for healthy patients (0 Gy) is around $10^{-3}$, we can see in Figure 2.10 that ADCI estimates a much higher aberration frequency at 0 Gy. While it is possible to reduce that number of false positives by imposing stricter thresholding in the MC-DC classifier, this tends to skew the calibration further down for high doses. The
<figure><img src="image_39.png" /><figcaption>Figure 2.11 – Figure displaying the improvements brought by applying false positive filters implemented in the 2017 ADCI version introduced in Liu *et al.* in [84]. A) shows the calibration curve on a sample prepared by Health Canada (HC), while B) shows the same for a sample prepared by Canadian Nuclear Laboratories (CNL). Green curves show ADCI performance without false positive filters, while cyan curves show the performance after DC reclassification using false positive filters. (Source: [84])</figcaption></figure>

following paragraph taken from [88] explains the relationship between false positives and
dose estimation very well:

" False positives have a proportionately larger impact at low doses, causing a stronger response in ADCI curves. At higher levels of radiation exposure, false positives tend to be balanced out by an increase in false negatives, i.e. DCs missed by the SVM. Beyond the level at which the number of false negatives exceed the number of false positives, the ADCI calibration curves exhibit a weaker response than manually derived curves. "

Following versions of ADCI implemented changes that improved this false positive rate.
The 2017 version of ADCI introduced in Liu *et al.* [84] introduced the MD-DC classifier
improvements described in Section 2.4.2. Furthermore, metaphase quality has a large
impact on MC-DC classification performance. Therefore, this version of ADCI provides
improvements in metaphase selection. To address this issue, a set of filters were designed
to rank metaphases based on morphological features like average length-width ratio of
chromosomes, centromere density (images containing closely clustered chromosomes are
<figure><img src="image_40.png" /><figcaption>Figure 2.12 – Figure displaying the improvements brought by applying false positive filters and manually curating metaphases on the Health Canada (HC) sample. Although [84] implements metaphase selection algorithms, this figure displays the impact of manual metaphase selection. Green curve is evaluated on a un-curated dataset, without false positive filters. Red curve is curated, but estimated without false positive filters. Cyan curves uses false positive filters but is evaluated on an un-curated dataset. Finally, blue curve is evaluated on a curated dataset with false positive filters.</figcaption></figure>

ranked lower), chromosome contour smoothness, total object count (metaphases with less
than 35 objects or more than 50 are ranked lower), ratio of identified to unidentified
objects, etc. A global image score was built as a linear combination of the previously
mentioned image features. The weights of this linear combination were optimized with a
grid search by minimizing the error to a known calibration curve. Finally, this weighted
sum was used to rank images to be processed by ADCI. The improvements caused by
re-classification of false positives can be seen in Figure 2.11. This figure also shows the
relationship between sample provenance and calibration curve shape, as the radiation
response of the Health Canada (HC) sample is almost linear, while the quadratic shape of
the calibration curve is much more visible for the Canadian Nuclear Laboratories (CNL) sample.

False positives reclassification and image selection lead to a very low level of false positives, and prevents the usual over-estimation of low doses in ADS systems. False positives can be further reduced by curating metaphases. The impact of this curation process can be seen in Figure 2.12.

the image selection models introduced in [84] rely on user-specified thresholds. In the 2019 version of ADCI introduced in [89], Li et al. demonstrate a procedure to find the parameters of those image selection models. A laboratory builds a calibration curve dataset comprised of a set of slides irradiated at different doses. As the ground-truth dose is known for each sample, it is possible to evaluate the dose estimation error. The parameters of the image selection models described in [84] are enumerated exhaustively, and the best model is chosen with respect to the calibration curve fitting error. Dose estimation is therefore completely automated, as long as a reference sample from a laboratory can be made available. On a laptop equipped with a GTX 960M, this implementation of ADCI can process around 100 metaphases per minute.

## 2.5 Deep learning for biological dosimetry

While all the solutions described in this section rely on conventional computer vision techniques, recent papers have used deep learning to solve the DC detection task at the core of biological dosimetry.

In [90], Jang et al. evaluate Faster R-CNN as a model for DC detection. First, metaphases are selected with a ResNet classifier, to reject images containing two metaphases, under-spread metaphases or metaphases with exposition issues. Once this image selection step is finished two Faster R-CNN-based models are trained: the Counting Network (CN) and Identification Network (IN). Dicentric yield is estimated using the number of chromosomes counted by the CN and the number of DC counted by the IN. Figure 2.13 shows the objects detected by the two models.

Both the Precision and Recall for the IN are around 90%, which means that around 10% of DC candidates are false positives, which is a relatively high false positive rate for biological dosimetry. Therefore, while the pipeline is accurate for high doses, the performance is unsatisfactory at low doses (below 1 Gy).

This pipeline is used to estimate a calibration curve on an independent sample (see
<figure><img src="image_41.png" /><figcaption>Figure 2.13 – Figure displaying the input image (left) detections of the Counting Network (CN) in the middle and Identification Network (IN), on the right.</figcaption></figure>

Figure 2.14). Note that in this plot, the "dicentric rate" (and not the usual dicentric yield) is plotted as a function of dose. The dicentric rate corresponds to the probability for a single chromosome to be dicentric. The dicentric yield, on the other hand, is the average number of DC per metaphase. Assuming that every metaphase contains 46 chromosomes, the dicentric rate should be multiplied by 46 to get the usual dicentric yield. At 0 Gy, the model presented in [90] achieves a dicentric yield of around 0.01, which is approximately 10 times higher than the manual dicentric yield (depending on laboratory counting protocol). Therefore, low doses are overestimated.

The approach proposed by Jang *et al.* partially circumvents the sophisticated metaphase selection algorithms presented in Section 2.4. Indeed, Jang *et al.* treat wether or not any chromosome is a DC as a Bernoulli-distributed random variable. Therefore, the radiation dose is linked to the mean of this distribution, instead of the average number of DC cell. Conventional calibration curves that link dicentric yield and dose rely on the implicit assumption that every metaphase contains 46 chromosomes which makes metaphase selection crucial. Incomplete metaphases quickly lead to an erroneous dicentric yield estimation: if an incomplete metaphase contains 39 chromosomes and 1 DC, a predicted DC count of 1 is likely to be biased downwards as the true (unknown) DC count can only match or exceed 1. In this scenario, Jang *et al.* effectively weight the DC count of 1 by a factor of $39/46$. As long as one can count chromosomes accurately and detect overlapped chromosomes estimating the dicentric rate is always possible even if some metaphases are incomplete.

Note that however, the training dataset of both neural networks relies on metaphases annotated by humans. As those metaphases were selected with human quality standards (good spread, object counts consistent with a single cell), some metaphase selection is still
<figure><img src="image_42.png" /><figcaption>Figure 2.14 – Figure displaying the calibration curve estimated by the pipeline proposed by Jang *et al.* in 2021 in [90]. This curve displays dicentric rate instead of dicentric yield. Multiplying dicentric rate by 46 (number of chromosomes in a normal metaphase) to get dicentric yield.</figcaption></figure>

needed to ensure that the training and evaluation distribution match closely. However, this step is considerably simplified as only a ResNet-based binary classifier is needed in the proposed pipeline. Even if an incomplete metaphase is accepted the counting network weights the DC count accurately instead of completely rejecting the metaphase. This prevents the rejection of some DCs and improves recall which is especially relevant for low doses as the blood sample might only contain a single-digit number of DCs. Overall, the importance of metaphase selection is reduced but the step is not eliminated completely although it certainly leads one to wonder whether or not it would be feasible.

Another approach was proposed by Shen *et al.* in [91]. Instead of training a model in an end-to-end fashion to detect bounding boxes, chromosomes are detected with conventional computer vision techniques and then classified with a ConvNet. First, the metaphase image is segmented with K-Means. Chromosomes that are not adequately separated are identified based on their morphology, and individual chromosomes are separated with a watershed algorithm. Next, a binary classifier is trained with stochastic gradient descent to solve the MC-DC classification problem. The authors of [91] do not estimate a calibration
curve with their DC detection method which makes objective comparison with other
approaches in this Chapter difficult.

## 2.6 DCScore: evaluating an ADS in a semi-automatic regime

2.6.1 Evaluating DCScore in realistic scenarios

ADCI aims to provide a fully-automated ADS. This is challenging as the number of
false positives need to be controlled very tightly to prevent an overestimation of the dose.
In a series of papers, researcher at LRAcc evaluated DCScore, a commercial ADS system
provided by MetaSystems in a semi-automatic regime. An automated metaphase finder
is used, and metaphases are selected by hand based on their quality. DCs are detected
automatically and undergo a manual review by an expert. False positives are rejected
by hand during this review. While this approach still need human experts, expertise
time is reduced. It provides a good compromise between manual scoring (MS) and fully
automated approaches like ADCI.

Within the context of a mass exposure where physical dosimetry information is not
available, the IAEA established guidelines to ensure that victims receive adequate care
based on the received dose [14]. Once a blood sample is taken and the cell culture step is
done (as shown in Figure 2.2), an initial manual scoring of 50 metaphases (called the triage
step) allows one to classify victims in range of doses. This first step is mostly useful to
identify patients who received large doses. Once this is done, 500 images should be scored
to refine the dose estimation. As manual scoring is very time-consuming, the triage step
ensure that therapeutically relevant data is produced as soon as possible after exposure
to radiation.

In real-world cases, the triage step has been shown to almost always underestimate
the real dose as shown by Voisin et al. in 2001 [92]. ADS systems provide an alterna-
tive approach to this two step process, as the scoring time can be drastically reduced. In
2009, Vaurijoux et al. [93] compared the triage step to a different strategy where several
hundreds of metaphases were scored using DCScore. The number of metaphases scored
automatically was chosen so that the analysis time is roughly the same for the two strate-
gies. Both metaphases and DCs were selected semi-automatically. In both cases, DCScore
made suggestions, and a human operator performed a review. Vaurijoux et al. use real
data from 46 victims of the Dakar 2006 accident, where an iridium-192 source from a gammagraphy instrument was not properly returned to its shielded container. The comparison between the DCScore and manual scoring calibration curve is shown in Figure 2.15.

<figure><img src="image_43.png" /><figcaption>Figure 2.15 – Calibration curves for manual scoring (bold black dashed line) and semi-automated scoring (bold black line) with manual review to suppress false positives. (Source: [93])</figcaption></figure>

Note that this calibration curve showcases significant differences from calibration curves produced by other ADS systems like ADCI, as seen in Figure 2.10. This is because of the human review, which suppresses false positives very effectively. However, this leads to a systematic underestimation of doses. This is because the recall is low which leads to an underestimation of the true number of DCs.

Furthermore, the authors compare DCScore and manual scoring in their ability to classify patients in range of doses. First, they note that misclassification during triage is common with the 2006 Dakar data: 50% of patients were classified in a different dose range depending on whether scoring was performed on 50 metaphases of 250. While this may seem very high, this needs to be recontextualized as the dose range for the Dakar accident is low, and it is likely that triage would perform better in higher dose scenarios.

On the contrary, dose range classification of DCScore and manual scoring are in agree-
ment in 96% of cases. Like with ADCI, slide-to-slide quality variations are identified as a key factor of performance variation for DCScore. Again, shown in Figure 2.15 the low-dose context of the 2006 Dakar accident is where DCScore shows the best performance, i.e., the lowest difference to the manual scoring calibration curve.

In [93], expositions are assumed to be uniform. However, accidental expositions may also be partial. In the uniform case, the number of DCs follows a Poisson distribution, as explained in Section 2.1. However, in the partial exposition case, the blood sample is effectively a mixture of an unexposed sample and an exposed sample with an unknown mixing factor. The distribution of aberration count follows a Zero-Inflated (ZIF) distribution, where the probability of a cell containing zero aberration is much higher than in the Poisson case. Wether or not DCScore can provide accurate triage in the case of partial expositions was studied by Vaurijoux et al. [94] in 2012 using in vitro experiments. Irradiated and un-irradiated samples are mixed in different proportions and the ability of DCScore to provide an accurate dose estimation in this challenging context is evaluated. The authors again observe high agreement between manual scoring results on 500 metaphases and ADS. Furthermore, as DCScore lets the authors process a much larger number of cells than MS, they found that "ADS is able to detect a dose significantly different from zero in samples having a lower fraction of irradiated cells compared to the MS method".

The results of [94] were refined in 2013 by Gruel et al. in [95] where a large scale emergency is simulated in vitro. The authors provide finer bounds on the number of cells needed for accurate triage with DCScore compared to manual scoring. As the suspected dose gets lower, the number of metaphases to score in order to accurately estimate the dose should be larger. For doses above 1 Gy, 1500 images scored with DCScore are sufficient, but lower doses especially in the context of a partial exposition may require 3000 images to be accurately estimated. Results from several european laboratories using this semi-automated approach for large scale expositions were compared in 2013 by Romm et al. [96]. The overall biological dosimetry pipeline (from blood sampling to dose estimation) was found to be sufficiently normalized, so that the difference in calibration curve estimation between laboratories was minimal.

## 2.6.2 Exploring DCScore shortcomings

Metafer has several tunable parameters for metaphase discovery, like microscope scanning speed and detector sensitivity. In this section, we discuss two recent papers (2019 and
<figure><img src="image_44.png" /><figcaption>Figure 2.16 – Calibration curves for automated (blue) and semi-automated (orange) DC-Score. Left panel shows calibration curve for X-rays, right panel for γ-rays. Note that the authors did not fit a linear-quadratic model to the dicentric yield point cloud, which explains the unusual look of the calibration curve compared to other examples in this thesis. (Source: [97])</figcaption></figure>

2020) exploring the relationship between metaphase selection performance and accurate dose estimation. In [97], Ryan *et al.* evaluated the relationship between those parameters and dose estimation accuracy. Furthermore, they compared automated DC detection with semi-automated detection. They found that for high doses (between 1 and 3 Gy), DCScore can accurately estimate doses in an automated fashion, and is good enough for preliminary triage of patients. Furthermore, manual review of metaphase images does not improve performance. Like previous papers [93]–[95], they find that the false positive rate is too high for accurate low dose estimation, and that high doses tend to be underestimated because of the insufficient recall. This pattern is very clear in Figure 2.16.

The performance of DCScore was found to depend strongly on the accurate rejection of unusable images (images containing several cells, underspread cells or poorly exposed cells). This is because dicentric yield estimation relies on the assumption that metaphases are complete, i.e they contain all 23 pairs of chromosomes. However, automatically detected metaphases may be incomplete because of failures in the metaphase detection software. A lower number of chromosomes reduces the actual probability of observing a DC in a specific cell, which means that a larger number of cells does not contain any dicentric. This can be wrongly interpreted as a partial exposition. In [98], Endesfelder *et al.* use the chromosome count provided by DCScore to adjust the uncertainty around the dose estimated by DCScore. This strategy is analogous to the one used by [90] and de-
<figure><img src="image_45.png" /><figcaption>Figure 2.17 – Figure showing the impact of chromosome count on calibration curves. [98] uses three different blood samples, shown respectively with blue, red and green calibration curves. Black calibration curve corresponds to pooled dataset. Chromosome counts are shown in lower right corner of all quadrants. (Source: [98])</figcaption></figure>

scribed in Section 2.5. The authors confirm that overdispersion is present (mean is greater
than variance) at high doses, suggesting that correcting dicentric yield estimation with
chromosome count is necessary. Excluding metaphases with low chromosome counts tends
to reduce overdispersion. This effect is shown in Figure 2.17.
## 2.7 Conclusion

In biological dosimetry, a radiation dose is inferred by estimating a dicentric yield. This yield corresponds to a proportion of DCs per cell. Accurately identifying DCs is a challenging computer vision problem, which has received significant attention over the years, even going back to the late 1980s.

Those ADS systems relied on a pipeline of object detection and object classification. DCs are identified based on their number of centromeres, which are identified by counting constrictions in chromosomes width profile. The Edinburgh dicentric hunter is an example of those early ADS systems [70]. Over the years improvements in computer vision and pattern recognition led to better performing ADS systems. The ADCI is the leading example for fully automated ADS systems, and received continuous upgrades over 7 years, culminating in Li et al. [89], which was published in 2019. Another example is DCScore, which was not described as precisely as ADCI. However, its potential as an ADS system in emergency situations was assessed by several different labs, as described in Section 2.6.

All of the ADS examples presented in this section are pipelines of processing step. Object detection and classification are treated in consecutive, separate step. On one hand, this offers unparalleled granular control over every step. This makes each step interpretable, and it is usually very easy to know why a specific object is misclassified from its morphological features. On the other hand, a considerable amount of design work was needed. Furthermore, errors in pipeline tend to compound: a segmentation error will often lead to a classification error, which could lead to a dose estimation error. Therefore, each step must be highly reliable, as the performance of every algorithm is strongly correlated.

While DCScore does not incorporate sophisticated metaphase selection algorithms for fully automated operations, it has been shown to be very effective assistants to human experts. In a series of papers published by LRAcc researchers [93]–[95] it was shown that human post-processing of DCScore detections could provide usable dose estimation faster than completely manual scoring. However manual scoring is more accurate, especially for low doses. This remains true for challenging contexts like partial exposition detection.

More recently, several papers investigated ConvNets for DC detection. In a 2023 paper [91], Shen et al. propose a conventional pipeline where chromosomes are first detected and then classified. The binary MC-DC classifier is a ConvNet trained with SGD. While the authors demonstrate high levels of precision and recall they do not present any calibration curve estimate. In a 2021 paper [90] Jang et al. trained Faster R-CNN to detect DCs.
While annotating metaphases with bounding boxes is costly, this removes the need for a
sophisticated chromosome detection step.
<figure><img src="image_46.png" /></figure>
TWO-STAGE CHROMOSOMAL
ABERRATION DETECTION WITH PATCH
CLASSIFICATION

In this chapter, we work with a dataset of chromosome patches as a preliminary study.
Building such a patch classification dataset is much easier than the spatially labelled
datasets used later in this thesis. However, as we will see below, building components of
an ADS system using only this dataset is very challenging. We focus on several different
potential uses for this dataset within the context of an ADS system.

First, we investigated the performance of deep classifiers, and compare it to conven-
tional methods. While reaching a high level of performance (even in a semi-supervised
regime) on this classification dataset is relatively easy, it does not transfer to real metaphase
images. We provide potential explanations for this fact.

If we can learn to synthesize chromosomes in an unsupervised fashion, we would be
able to produce synthetic, annotated metaphases. This could be used to pre-train object
detection models. We investigate two chromosome simulation strategies, VAEs and a
specific GAN architecture called pix2pix (see Section 1.2.1). Both of those strategies
struggle with synthesizing diverse and realistic chromosomes without artifacts.

Finally, given our lack of success in synthesizing satisfying chromosomes, we build a
synthetic annotated metaphase dataset from chromosome "skeletons". Using this dataset,
we demonstrate that model performance depends on the number of images in the training
set, and on the dose. Models trained on one given dose usually do not reach the same
level of performance when evaluated on a different dose.
## 3.1 Datasets

### 3.1.1 Patch dataset

Our dataset is comprised of 17061 grayscale chromosome patches of varying sizes. This dataset contains two classes, MCs and DCs. 6349 of those chromosomes are DCs, while 10712 are MCs. In the remainder of this chapter, we denote this dataset as $\mathcal{D} = \{(x_1, y_1), \cdots, (x_N, y_N)\}$, where $N$ is the number of images. $\mathbf{X} = \{x_1, \cdots, x_N\}$ is the set of input image $\mathbf{y} = \{y_0, \cdots, y_N\}$ is the set of labels, with $y_k \in \{0, 1\} \forall 1 \le k \le N$. $y_k = 0$ if the $k$-th image is a MC, and 1 if it is a DC

The chromosome patches underwent a significant amount of pre-processing. First, the images where segmented with K-Means to separate background and foreground. K-Means was chosen instead of a global thresholding method because it dealt with the intensity differences of chromosomes better. For example, the region between sister chromatides usually exhibit higher intensity, and Otsu thresholding produces noisy segmentations in this case. Once this step was completed, secondary object are rejected with connected components, on the basis of their size. We keep the largest object in the patch, assuming that it was the chromosome object. Finally, the image intensity values are normalized between 0 and 1 in the following fashion:

$$ \tilde{x} = \frac{x - \min(x)}{\max(x) - \min(x)}, \qquad (3.1) $$

where $\tilde{x}$ is the normalized image.

### 3.1.2 Skeleton datasets

Using the aforementioned patch dataset, we built a "skeleton" dataset. First, for every chromosome, we extracted the centerline, and computed the distance between the boundary and the centerline. Centromeres were identified by locating minimas of this distance. A stylized chromosome representation was built using this information. A grid of chromosomes using this skeleton representation is shown in Figure 3.2. To ensure that every image has the same size, they are padded to a size of 128 by 128 pixels, i.e., $x \in [0, 1]^{H \times W}$, with $H = W = 128$. A grid of padded chromosome images is displayed in Figure 3.1. We also built a labelled metaphase datasets using the chromosomes skeletons introduced above. First, we extracted the locations of chromosomes in real metaphase images. Then,
<figure><img src="image_47.png" /><figcaption>Figure 3.1 – A grid of padded chromosome images.</figcaption></figure>

we sampled from our skeleton dataset at random and we overlaid those skeletons on a
white background, at the locations extracted previously. Chromosomes skeletons where
randomly rotated and resized, to increase metaphase diversity. Finally, we sampled a re-
alistic number of DCs using the manual calibration curve established by the LRAcc (see
Chapter 2). We built 3 different datasets at dose levels of 1 Gy, 2 Gy and 5 Gy. For every
dose level, we retrieved the corresponding average number of DC per metaphase using
<figure><img src="image_48.png" /><figcaption>Figure 3.2 – A grid of samples of chromosome skeletons.</figcaption></figure>

the manual calibration curve, and sampled a Poisson distribution with this mean to get a DC count per metaphase. In our first experiments with Faster R-CNN, we quickly found that if we labelled every MC, the training of the model falls into a local minima where minority classes are never predicted. Therefore, we only labelled chromosome aberrations and MCs are effectively treated as background. Figure 3.3 gives an example of cropped training image. For each dose, we built a test dataset of 5000 images.
<figure><img src="image_49.png" /><figcaption>Figure 3.3 – Cropped image of the synthetic dataset. The two DCs in this image are indicated by red bounding boxes.</figcaption></figure>

**3.2 Chromosome classification**

**3.2.1 Resnet-based classifier**

To establish a DNN-based, fully supervised chromosome classification baseline, we used
a ResNet model [8]. Our goal was to demonstrate that conventional DNN-based classifiers
could reach a high level of performance. This model is trained with Adam [12] (a variant
of SGD), with a learning rate of 0.0001, on the dataset D, as described in Section 3.1.1.
We did not use any data augmentation or learning rate scheduling. Furthermore, we did
not perform any architecture search or hyperparameter tuning. Therefore, we do not use
<figure><img src="image_50.png" /><figcaption>Figure 3.4 – Reconstruction grid of unseen chromosome patches. The first 2 rows are unseen input images, the last 2 rows show their reconstruction by the VAE. Notice that the reconstructions are smoother than the inputs. While the $L_2$ criterion is a common reason for this effect, the output of the same model trained with an $L_1$ criterion are very similar.</figcaption></figure>

any validation set, and the performance was evaluated on a test set comprised of 20% of the training data. The training dataset is comprised of the remaining 80 % of $\mathcal{D}$.
### 3.2.2 Chromosome patches autoencoder

We trained a convolutional VAE as a simple unsupervised feature learning technique. This model is trained on the dataset X with an $L_2$ reconstruction criterion and a KL weighting parameter of 1 (see Chapter 1, Section 1.2.1 for an overview of VAEs). We chose a latent space dimension of 256. In Figure 3.4, we show some examples of patch reconstruction on unseen data from this VAE.

### 3.2.3 Latent space classifier

The trained VAE model mentioned in Section 3.2.2 provides a way to map image samples into a lower dimensional latent space of dimension $d$. We set $d = 256$ in our experiments. Knowing the image labels associated with each low-dimensional feature vector, we can train simple binary classifiers, like a logistic regression model. Again, the goal here was not to demonstrate the highest level of performance possible. Other methods could have been studied, like gradient boosting, random forests or SVMs. Because of the high number of latent vectors, this logistic classifier was trained with a batch stochastic gradient method called SAGA [99]. We used the scikit-learn [100] implementation of this algorithm in our experiments.

### 3.2.4 Performance results and commentary

First, we computed the performance of our ResNet classifier baseline in a conventional fashion, with a test set comprised of 20% of overall images. Then we explored the potential of classifiers trained on the embeddings of a variational autoencoder. Our goal was to evaluate the potential of a logistic regression classifier trained on a small number of embeddings from a VAE. If this strategy succeeded it would be possible to train an effective MC-DC classifier using a large, unlabelled dataset of patches (to train a VAE) and a few labelled patches. Therefore, we intentionally made the training set of our logistic classifier small (at most 10% of the complete dataset). This classifier takes a VAE embedding as input and outputs a probability distribution over both classes. Furthermore, we also wanted to evaluate whether our logistic classifier was sensitive to shifts in the distribution of its training data.

Therefore, we built $k$ non-overlapping splits of the dataset $\mathcal{D}$. Each one of those splits $S_k$ is composed of image-label pairs, i.e $S_k = \{\{\mathbf{X}_k^T, \mathbf{y}_k^T\}, \{\mathbf{X}_k^E, \mathbf{y}_k^E\}\}$, where $\mathbf{X}_k^T \in \mathbb{R}^{n_k^T \times d}$ where $d$ is the latent dimension of the autoencoder, so that $n_k^T + n_k^E = N$. Here, $\{\mathbf{X}_k^T, \mathbf{y}_k^T\}$
<figure><img src="image_51.png" /><figcaption>Figure 3.5 – Description of our weakly-supervised experimental setup. The fully supervised ResNet baseline is trained on a simple train-test split. The VAE is trained in an unsupervised fashion (it does not use $\mathbf{y}$) and its performance is evaluated on a test split. Once the VAE is trained its encoder is used to compute the latent representation for all images in $\mathbf{X}$. We build $K$ different splits of this "embedded" dataset. Here, the training dataset is always less than 10% of the complete dataset. Performance is evaluated on the test set. As we have $K$ different splits, the mean and variance of performance can be computed, to evaluate the sensitivity of the logistic classifier to data sampling. Results are visible in Figure 3.6.</figcaption></figure>
is a training set while {$\mathbf{X}_k^E, \mathbf{y}_k^E$} is a test set. Finally, we have $\mathbf{X}_k^T \cap \mathbf{X}_k^E = \emptyset \forall k$. Figure 3.6 shows the results of this evaluation.

The highest performance model is the fully-supervised ResNet baseline, reaching an AUC of 0.98. However, it is possible to reach a very high level of performance with much less labelled data and a strong feature learner, as the solid blue, purple and green lines show. We trained a VAE on $\mathbf{X}$, and computed the latent representation of every image in $\mathbf{X}$. We then used the embeddings and the corresponding labels to train classifiers. Those classifiers are trained on respectively 10, 5 and 1% of the training data, and evaluated on the rest. One can get over 90% of the performance of the supervised baseline with only 10% of the labelled data ! The non-linear embedding provided by the VAE leads to higher downstream performance than the linear embedding of PCA, as the dashed curves show. Furthermore, there is more performance variation in the downstream models that use PCA embeddings. While this level of performance might lead us to think that "two-stage" models which successively detect and classify chromosomes might be very simple to build, in practice we did not succeed in reproducing this high level of performance with patches extracted from metaphases.

The first reason we identified is the lack of representativity. In our dataset, all mono-centric chromosomes correspond to chromosomes detected by DCScore (See section 2.6 for a review of papers examinating DCScore performance). We used a first patch dataset where the DCs were also extracted by DCScore. However, DCScore is not able to recover all DCs. For example, in some intensity conditions, the detection performance of DCScore worsens, which means that especially bright or dark DCs are not present in the dataset. This led to a phenomenon analogous to label leakage; it became possible to separate MCs and DCs with simple thresholds on the mean and variance of the patch intensity. We corrected this problem by enlarging the dataset with manually detected DCs, to build a more representative dataset. However, we still observed difference in validation performance between the two subset of our dataset, depending on wether the DCs were detected by DCScore or a human expert.

We identified another explanation for this low level of DC detection performance. The performance of the detection model and the classifier are strongly linked. For example, if the second centromere of a DC is very close to the end of this chromosome, slight segmentation errors might lead to a bounding box that does not include this second centromere. In this case, there is no chance of accurately classifying this chromosome as a DC. In the same way, rejecting overlapping chromosomes is essential to ensure a high level
<figure><img src="image_52.png" /><figcaption>Figure 3.6 – Figure comparing the performance of several models. The dashed redline indicates the performance of the ResNet supervised baseline. The solid blue line indicates the mean performance of a logistic classifier trained on VAE latent features, if 10% of the data is labelled. The purple and green lines depict the performance of the same model, for supervision rates of 5% and 1% respectively. The shaded area shows performance variation across all splits (plus / minus one standard deviation). Finally, the dashed lines depicts the same performance for PCA embeddings, instead of VAE.</figcaption></figure>