main.tex

\documentclass[
  digital,     %% The `digital` option enables the default options for the
               %% digital version of a document. Replace with `printed`
               %% to enable the default options for the printed version
               %% of a document.
%%  color,       %% Uncomment these lines (by removing the %% at the
%%               %% beginning) to use color in the printed version of your
%%               %% document
  oneside,     %% The `oneside` option enables one-sided typesetting,
               %% which is preferred if you are only going to submit a
               %% digital version of your thesis. Replace with `twoside`
               %% for double-sided typesetting if you are planning to
               %% also print your thesis. For double-sided typesetting,
               %% use at least 120 g/m² paper to prevent show-through.
  nosansbold,  %% The `nosansbold` option prevents the use of the
               %% sans-serif type face for bold text. Replace with
               %% `sansbold` to use sans-serif type face for bold text.
  nocolorbold, %% The `nocolorbold` option disables the usage of the
               %% blue color for bold text, instead using black. Replace
               %% with `colorbold` to use blue for bold text.
  lof,         %% The `lof` option prints the List of Figures. Replace
               %% with `nolof` to hide the List of Figures.
  lot,         %% The `lot` option prints the List of Tables. Replace
               %% with `nolot` to hide the List of Tables.
]{fithesis4}
%% The following section sets up the locales used in the thesis.
\usepackage[resetfonts]{cmap} %% We need to load the T2A font encoding
\usepackage[T1,T2A]{fontenc}  %% to use the Cyrillic fonts with Russian texts.
\usepackage[
  main=english, %% By using `czech` or `slovak` as the main locale
                %% instead of `english`, you can typeset the thesis
                %% in either Czech or Slovak, respectively.
  english, czech, slovak %% The additional keys allow
]{babel}        %% foreign texts to be typeset as follows:
%%
%%   \begin{otherlanguage}{german}  ... \end{otherlanguage}
%%   \begin{otherlanguage}{russian} ... \end{otherlanguage}
%%   \begin{otherlanguage}{czech}   ... \end{otherlanguage}
%%   \begin{otherlanguage}{slovak}  ... \end{otherlanguage}
%%
%% For non-Latin scripts, it may be necessary to load additional
%% fonts:
\usepackage{paratype}
\def\textrussian#1{{\usefont{T2A}{PTSerif-TLF}{m}{rm}#1}}
%%
%% The following section sets up the metadata of the thesis.
\thesissetup{
    date        = \the\year/\the\month/\the\day,
    university  = mu,
    faculty     = fi,
    type        = bc,
    department  = Department of Visual Computing,
    author      = Bruno Petrus,
    gender      = m,
    advisor     = {Assoc. Prof. RNDr. Martin Maška, Ph.D.},
    title       = {Segmentation of Membrane-Stained Cells in Image Data of Organoids},
    TeXtitle    = {Segmentation of Membrane-Stained Cells in Image Data of Organoids},
    keywords    = {segmentation, image processing, scikit-image, watershed,
                   thresholding, fluorescence microscopy, zero-crossing},
    TeXkeywords = {segmentation, image processing, scikit-image, watershed,
                   thresholding, fluorescence microscopy, zero-crossing},
    abstract    = {%
    In this thesis, we present a two-stage image analysis pipeline for segmenting membrane-stained and nuclei-stained image data of organoids captured using a fluorescence microscope. In the first stage, a coarse segmentation of individual cells in the membrane channel is obtained by analysing the Laplacian operator's response and watershed-driven tesselation. In the second stage, the coarse segmentation is projected into the nuclear channel, and individual nuclei regions are refined via local gradient thresholding. The performance of the developed pipeline is evaluated using a representative, manually annotated dataset, with a focus on detection and segmentation accuracy and execution time. The developed pipeline is implemented in Python and complemented with a detailed documentation of its usage. 
    },
    thanks      = {%
       I am extremely grateful to my advisor, Assoc. Prof. RNDr. Martin Maška,
        Ph.D., for his continuous guidance, helpful advice, and \mbox{feedback}. I am
       also thankful to my family and close friends for their support. \mbox{Finally},
       I am grateful to Dr. Zuzana Sumbalová Koledová and \mbox{Dr. Matea} Brezak for providing
       me with the annotated image data \mbox{analysed} in this thesis.
    },
    bib         = bibliography.bib,
    %% Remove the following line to use the JVS 2018 faculty logo.
    facultyLogo = fithesis-fi,
}
\usepackage{makeidx}      %% The `makeidx` package contains
\makeindex                %% helper commands for index typesetting.
%% These additional packages are used within the document:
\usepackage{paralist} %% Compact list environments
\usepackage{amsmath}  %% Mathematics
\usepackage{amsthm}
\usepackage{amsfonts}
\usepackage{url}      %% Hyperlinks
\usepackage{markdown} %% Lightweight markup
\usepackage{listings} %% Source code highlighting
\usepackage{graphicx}
\lstset{
  basicstyle      = \ttfamily,
  identifierstyle = \color{black},
  keywordstyle    = \color{blue},
  keywordstyle    = {[2]\color{cyan}},
  keywordstyle    = {[3]\color{olive}},
  stringstyle     = \color{teal},
  commentstyle    = \itshape\color{magenta},
  breaklines      = true,
}
\usepackage{floatrow} %% Putting captions above tables
\floatsetup[table]{capposition=top}
\usepackage[babel]{csquotes} %% Context-sensitive quotation marks
\usepackage{svg}
\usepackage{caption}
\usepackage{subcaption}

%% Specify new commands
\newcommand*{\R}{\ensuremath{\mathbb{R}}}
\newcommand*{\Z}{\ensuremath{\mathbb{Z}}}

\begin{document}
%% The \chapter* command can be used to produce unnumbered chapters:
\chapter*{Introduction}
%% Unlike \chapter, \chapter* does not update the headings and does not
%% enter the chapter to the table of contents. I we want correct
%% headings and a table of contents entry, we must add them manually:
\markright{\textsc{Introduction}}
\addcontentsline{toc}{chapter}{Introduction}

In the field of biology, the development of confocal laser scanning and
fluorescence microscopy has allowed scientists to create high-resolution multidimensional images of various plant and animal tissues \cite{stegmaier2016}. One of the areas where these methods are actively utilised is in skin cancer research \cite{gupta2020} or in the analysis of ovarian tissues from oncological patients \cite{fabbri2014}. Since cancer is the second-leading cause of death in the world \cite{mayo-clinic-cancer}, further research into this area and the creation of better tools for image processing are essential to developing cures for the disease.

One of the critical processes in researching this topic is the segmentation of cells to study their behaviour and shape; however, this often has to be done manually, taking valuable time from researchers to focus on more significant parts of their research. It is, therefore, necessary to develop robust automatic segmentation techniques that scale with the data's growing size.

This thesis focuses on developing an image analysis pipeline for segmenting membrane-stained and nuclei-stained cells in fluorescence microscopy image data of organoids without reference annotations. This pipeline is delivered in an easy-to-use form, and a quantitative analysis of the segmentation quality is performed. The developed method consists of two steps. In the first step, the program analyses the second partial derivate of the image data and uses the watershed algorithm to create a coarse segmentation of individual cells; then, in the second step, this is utilised to create a better-fitting segmentation by locally thresholding the nuclei-stained channel.

The first chapter defines the basic terms required for understanding the rest of the thesis. It also formalises the problem of edge detection, thresholding, and segmentation. The next chapter describes the available image data and the common issues that must be dealt with before analysing the data. In the third chapter, a high-level overview of all the steps taken in the segmentation process is detailed and visualised. The fourth chapter details the implementation of the non-trivial steps in the pipeline and provides the motivation behind choosing Python and the required libraries. The fifth chapter discusses how the quantitative analysis was done and its results, alongside the program's performance on the analysed dataset. Lastly, the appendix serves as a documentation for the numerous tweakable parameters.

\chapter{Basic terms}

Image processing involves mathematical techniques and algorithms
to extract useful information or somehow alter the image. The purpose of this
section is to provide a comprehensive and shallow overview of the fundamental
concepts used in the developed image processing pipeline.

\section{Digital Image}
Our eyesight is undoubtedly one of our most helpful senses, allowing us to
observe, study and comprehend the world around us. However, the way a computer
processes visual information is quite different. To enable us to study the world
using computers, numerous tools and algorithms have been developed. In this
section, I will begin by defining the concept of images in a formalized manner,
so that mathematics and computer science knowledge can be utilized to analyze
images.

Essentially we can think of an image as an n-dimensional function $f(x_1, x_2,
..., x_n)$, where $x_1, x_2, ..., x_n$ are coordinates inside a spatial domain,
while the function value at those coordinates specifies the magnitude or
intensity of the signal at that point. We can think about images as functions of
type $f:\R^n \rightarrow \R^m$ where $n$ indicates the number of spatial
dimensions and $m$ specifies the number of channels. For example, when we are
talking about two-dimensional grayscale images, $m$ is equal to 1 and $n$ is
equal to 2, while a typical coloured photograph has 3 colour channels. Another
important image type is the binary image, where the function's range is just two
values.

Once we start using computers to work with images, we no longer can work with
real-valued spatial dimensions and intensities, because we only have a finite
number of bits. The process of discretizing images in a finite grid and assigning
intensities from a finite range is called sampling and quantization. In essence,
we create a digitalized version of the original signal, which can be represented
as a one- or \mbox{multi-dimensional} array of values.

\section{Edge detection}

In digital image processing, edges typically separate different regions of
interest from one another. They are essential in digital image processing
as they can be used in many workflows --- such as image segmentation and
classification --- to extract various interesting features and regions
from the image.

There is no authoritative definition of an edge, but it can be thought of as a
set of connected pixels that lie on the boundary between two regions
\cite{gonzalez2002}. Usually, there is a gradual difference in intensities
between two regions, as is shown in Figure \ref{fig:edge_intensities}. It is
hard to define precisely where an edge begins and ends, so the first and the
second derivative of the image is usually calculated and analysed.

Looking at the second derivate of the image, we can see two distinct spikes,
which are called zero-crossings, and they will be utilised by the program for finding edges.

\begin{figure}
    \begin{center}
        \includegraphics[width=6cm]{resources/inkscape/gradient.png}
    \end{center}
    \caption{Intensity of the signal, and its first and second derivative.}
    \label{fig:edge_intensities}
\end{figure}

\subsection{Gradient operator}

The first common way of working out where the edge is located is by looking at
the gradient of the image. The gradient of an image $\nabla f(x,y)$ is defined
as a vector \cite{gonzalez2002}:
$$\nabla f =
\begin{bmatrix}
    \dfrac{\partial f}{\partial x}\\[2ex]
    \dfrac{\partial f}{\partial y}
\end{bmatrix}.$$
This vector represents the rate of change in the image intensities, and it points in
the direction of the maximum rate of change of intensity. Often just the bare
vector is not of interest; instead, the magnitude of the gradient, denoted as
$|\nabla f|$, is calculated too. Confusingly enough, sometimes this is also called
the gradient in the literature, but here I will differentiate it by putting
absolute value around it. Formally it is defined as:
$$|\nabla f| = \sqrt{\left(\frac{\partial f}{\partial x}\right)^2 +
\left(\frac{\partial f}{\partial y}\right)^2}.$$

Now the question is how to calculate the partial derivatives of some digital image.
It is impossible to calculate the partial derivatives using the standard
definition from mathematical analysis because digital images are discrete; thus,
we must use an approximation. To calculate the $\frac{\partial f}{\partial x}$
we can use the following approximation:
$$\frac{\partial f}{\partial x}(i, j) \approx \frac{f(i+1,j)-f(i-1,j)}{2}$$
and similarly for the $\frac{\partial f}{\partial y}$ we can use:
$$\frac{\partial f}{\partial y}(i, j) \approx \frac{f(i, j+1)-f(i,j-1)}{2}.$$

Using them, we can reasonably approximate the magnitude of the gradient as:
$$|\nabla f(i, j)| \approx \frac{1}{2}\sqrt{(f(i+1, j) - f(i-1, j))^2 + (f(i, j+1) - f(i, j-1))^2}.$$

An example of the magnitude of the gradient is visible in Figure
\ref{fig:grad-thresh-grad}.
We can then extract the edges with further processing but since the program
does not rely on this method, I will skip the details; however, the program uses
the magnitude of the gradient for calculating the gradient threshold, described later
in this chapter.

\subsection{Laplacian operator}

The second common way --- and the one we will use for finding edges --- is using the second derivative, or more
precisely, using the Laplacian operator. For a 2D function $f(x, y)$ the Laplacian operator is
defined as \cite{gonzalez2002}:

$$\nabla^2 f = \frac{\partial^2 f}{\partial x^2} + \frac{\partial^2 f}{\partial y^2}.$$

However, in the context of digital image processing, our signal is discrete;
therefore, a discrete approximation is needed again. A fairly standard way is to use
the following approach to calculate the partial second derivative in the $x$
axis:
$$\frac{\partial^2 f}{\partial x^2}(i, j) \approx f(i + 1, j) + f(i - 1, j) - 2f(i, j)$$
and similar in the y axis:
$$\frac{\partial^2 f}{\partial y^2}(i, j) \approx f(i, j + 1) + f(i, j - 1) - 2f(i, j).$$
Substituting into the prior definition, we get the following approximation:
$$\nabla^2 f \approx f(i+1, j) + f(i-1, j) + f(i, j+1) + f(i, j-1) - 4f(i,j).$$
To calculate the Laplacian of an image, this formula is applied to every pixel.
The outer pixels can be dealt with in many ways, but the developed program simply
ignores them.


One of the disadvantages of this approach is the sensitivity to noise. For
demonstration purposes, the effects are shown in 1D in Figure
\ref{fig:filtering-noise}. While the original
signal does not look too degraded, it is impossible to find the correct
zero-crossings. Therefore we have to help the operator by first smoothing the
image. A common approach is first to use Gaussian smoothing; this is sometimes
called a Laplacian of Gaussian operator or a Mexican hat filer. The idea is that
by averaging, we reduce some of the effects of the noise
\cite{hipr-mexican-hat}. As is shown in the figure, without the additional Gaussian
smoothing, we cannot deduce much from the second derivative, but after the smoothing
is applied, the zero-crossing becomes visible.

\begin{figure}
    \begin{subfigure}[t]{0.48\textwidth}
        \frame{\includegraphics[width=\textwidth]{./resources/small_noise.png}}
    \end{subfigure}
    \begin{subfigure}[t]{0.48\textwidth}
        \frame{\includegraphics[width=\textwidth]{./resources/filtered_noise.png}}
    \end{subfigure}
    \caption{Demonstration of the effect of noise and filtering on a simple
    linear increase in signal intensity. In the second row, the first derivate of the
    signal is displayed, and in the third row is the second derivative. On the
    left, an artificial noise has been added to the signal. On the right, a
    Gaussian smoothing has been applied on top of the noisy signal.}
    \label{fig:filtering-noise}
\end{figure}

After applying the Laplacian operator, we can find the zero-crossings and the edges. 
A simple algorithm is to go through the whole image and checks
if a particular pixel at a given position is positive while its neighbours are
negative. If this condition is true, we set this position to white in the
resulting image and black otherwise. This creates the following image in
Subfigure \ref{fig:zc}. Here is another valuable property of the zero-crossing
algorithm. The produced edges are thin and form closed contours. This algorithm is
described in more detail in Section \ref{section:zc}.

\section{Thresholding}

Another fundamental area of image processing is thresholding. During thresholding
we try to find a suitable threshold value which separates logical parts of the image,
such as the background and foreground.

\subsection{Otsu method}

The Otsu thresholding method is a popular technique to determine a threshold
value between two classes (usually the foreground and the background). As such,
it works well if the image's histogram has a clear bimodal distribution. The
main idea of this approach is to see the problem of finding the threshold as an
optimization problem, where the inter-class variance of the two classes is to be
minimized. More information is available in \cite{otsu1979}. A demonstration of the
Otsu thresholding is shown in Figure \ref{fig:otsu-demonstration}.

\begin{figure}
    \begin{subfigure}[t]{0.45\textwidth}
        \centering
        \frame{\includegraphics[width=\textwidth]{resources/otsu-orig.png}}
        \caption{}
        \label{fig:otsu-orig}
        %\caption{}
    \end{subfigure}
    \begin{subfigure}[t]{0.45\textwidth}
        \centering
        \frame{\includegraphics[width=\textwidth]{resources/otsu-thresholded.png}}
        \caption{}
        \label{fig:otsu-thresholded}
        %\caption{}
    \end{subfigure}
    \begin{subfigure}[t]{0.6\textwidth}
        \centering
        \includegraphics[width=\textwidth]{resources/otsu-histogram.png}
        \caption{}
        \label{fig:otsu-histogram}
        %\caption{}
    \end{subfigure}
    \caption{Demonstration of the Otsu thresholding method. In (\subref{fig:otsu-orig}) an image from the membrane-stained channel is shown. In (\subref{fig:otsu-thresholded}) the Otsu thresholding is applied. In (\subref{fig:otsu-histogram}) the histogram of the image data is shown with a line indicating the threshold.}
    \label{fig:otsu-demonstration}
\end{figure}

%todo: Finish

\subsection{Gradient thresholding}

Another useful thresholding method is the gradient thresholding algorithm. Unlike
many other thresholding algorithms, it does not purely rely on the analysis of
the image's histogram; instead, the threshold is calculated as a weighted sum of
the image's intensities, where the weights are dictated by the magnitude of the
gradient \cite{pb130}.

The idea is that the background tends to be rather uniform compared to the
object in the foreground; hence, the magnitude of the gradient is the greatest
in the foreground, and the edge separating it from the background. That can be
utilized by using the magnitude as the weight in the summing of intensities.
This way the background's intensities play lesser role.

Formally the threshold $a$ is defined as:
$$a = \sum_{u, v} I(u, v) \frac{|\nabla I(u, v)|}{\sum_{i, j} |\nabla I(i,j)|}$$
where $I(i, j)$ stands for the image intensity at the $(i, j)$ position.

\begin{figure}
    \begin{subfigure}{0.3\linewidth}
        \centering
        \includegraphics[width=\textwidth]{./resources/grad-thresh-orig.png}
        \caption{}
        \label{fig:grad-thresh-orig}
    \end{subfigure}
    \begin{subfigure}{0.3\linewidth}
        \centering
        \includegraphics[width=\textwidth]{./resources/grad-thresh-grad.png}
        \caption{}
        \label{fig:grad-thresh-grad}
    \end{subfigure}
    \begin{subfigure}{0.3\linewidth}
        \centering
        \includegraphics[width=\textwidth]{./resources/grad-thresh-final.png}
        \caption{}
        \label{fig:grad-thresh-final}
    \end{subfigure}
    \caption{Demonstration of the gradient thresholding. Image data from the
    membrane-stained channel (\subref{fig:grad-thresh-orig}) was thresholded
    (\subref{fig:grad-thresh-final}) using the gradient thresholding. In
    (\subref{fig:grad-thresh-grad}) the magnitude of the gradient is shown.}
    \label{fig:demo_grad_thresh}
\end{figure}

In Figure \ref{fig:demo_grad_thresh} an example image is shown. There, the
gradient thresholding algorithm was used to segment the cell from the
background. As can be seen, the magnitude of the gradient is biggest in the
cell; hence, it contributes more to the value of the threshold.

\section{Mathematical morphology}

Mathematical morphology is an area of research studying how to extract valuable
information and how to filter images. A few examples of what mathematical
morphology is used for are edge-detecting algorithms, finding convex hulls, and image
filtering. In the simplest terms, mathematical morphology looks at images as
sets of numbers and modifies them using structural elements. Those can take many
forms, but their shape varies based on the task we are trying to achieve. The
main use of mathematical morphology in this thesis is to pre- and post-process
image data generated by other image processing methods. First, some basic
operations on binary images are described. Note that this area is vast and still
being actively researched, and I will cover only parts which are helpful later
in my prototype.

In this branch of image processing, we look at the images through set theory
instead of looking at images as discrete-valued functions. In this context, an
image is a subset of 2-D integer space $\Z^2$, which can easily be generalized
to more dimensions. Each element of this subset $(x, y)$ represents a present
pixel at those coordinates in the image. In other words, if we look at an image
as an array of rectangular shape, then each foreground pixel's coordinates will
be included in the set representation.

\subsection{Preliminary}
Before we look at the two main operations --- erosion and dilation --- we have to
define two not-so-common operations on sets which are heavily used in
mathematical morphology.

During the following sections, we need to be able to talk about how to create a
reflection of the image around its origin, as shown in Figure
\ref{fig:morp_refl}. Reflection  of a set $B$ is formally defined as \cite{gonzalez2002}:
$$\hat{B} = \{w | w=-b \text{ , for } b \in B\}$$
\begin{figure}
    \begin{center}
        \includegraphics[width=0.5\textwidth]{resources/inkscape/morphology-reflection.png}
    \end{center}
    \caption{Morphological reflection of image $A$, where $C=\hat{A}$} % todo: replace
    \label{fig:morp_refl}
\end{figure}

While manipulating images we commonly have to move parts of the image around.
In principle, a similar operation, called translation, is defined on sets. The
formal definition of moving a set $A$ by some translation vector $z$ is
\cite{gonzalez2002}:
$$(A)_z = \{c | c = a + z\text{ , for } a \in A\}$$
The operation is illustrated in Figure \ref{fig:morph_translation} where the set
$A$ is moved downwards by a vector $z$.
\begin{figure}
    \begin{center}
        \includegraphics[width=0.5\textwidth]{resources/inkscape/morphology-translation.png}
    \end{center}
    \caption{Morphological translation}
    \label{fig:morph_translation}
\end{figure}

\subsection{Dilation and Erosion}
Now it is possible to properly define two fundamental operations in the field of
mathematical morphology, upon which many other operations are built. Both of
these operations include a structural element. A structural element is again a
set, but it can be used to describe how far and in which way a particular
operation will influence the result.


\subsubsection{Dilation}
Intuitively we can think about the dilation operator as a way to dilate spots
inside an image into a bigger space. Formally it is defined as
\cite{gonzalez2002}
$$A \oplus B = \{z | (\hat{B}_z) \cap A \neq \emptyset\}$$
where $A$ and $B$ are both sets in $\Z^2$.

An example of dilation can be seen in Figure \ref{fig:morph_dilation}. In the
first row, we have a square image with a side length of $d$ and a structuring
element whose side length is a fourth of the original. If we apply the dilation
operation on these two images, we get the one on the right. You can imagine
taking the smaller square and gliding it over the bigger one. If there is an
overlap, we can add the centre of the structural element to the result.

Closer to the mathematical definition, all possible translations are applied to
the structural element and if there is an overlap, we add the centre of the
structural element to the resulting set.

For example, this operation can be used to join disconnected segments in an
image by choosing a large enough structural element. Its size will depend on the
particular example and does not necessarily have to be of regular size \cite{gonzalez2002}.

\begin{figure}
    \begin{center}
        \includegraphics[width=6.3cm]{resources/morph_dilation.jpg}
    \end{center}
    \caption{Morphological dilation (Gonzalez 2002)} % todo: Replace
    \label{fig:morph_dilation}
\end{figure}

\subsubsection{Erosion}

The next essential operation is erosion. As the name suggests, we use it to
erode parts of the image which we consider unnecessary. An erosion $A \ominus B$
of two sets $A$ and $B$, both being part of $\Z^2$ space, is defined as
\cite{gonzalez2002}: $$A \ominus B = \{z | (B)_z \subseteq A\}$$

Again $A$ can be thought of as an image while $B$ can be looked at as a
structural element. In contrast to the dilation operation, we are now looking if
the whole structural element fits inside $A$ as can be seen in Figure
\ref{fig:morph_erosion}. There a smaller structural element is used to cut away
parts of the bigger square.

As opposed to the dilation, we can use erosion to divide lightly connected
segments, and to suppress noise and other unwanted disturbances in an image.\cite{gonzalez2002}

\begin{figure}
    \begin{center}
        \includegraphics[width=6.3cm]{resources/morph_erosion.jpg}
    \end{center}
    \caption{Morphological erosion (Gonzalez 2002)} % todo: Replace
    \label{fig:morph_erosion}
\end{figure}

\subsection{Opening and closing}

Now that the basic two operations have been defined, we can move on to two
very useful operators, morphological opening and morphological closing.

\subsubsection{Opening}

The morphological opening is defined as erosion followed by a dilation with a
reflected structural element, formally the opening $\circ$ of a set $A$ by a structuring
element $B$ is defined as\cite{soile2004}:

$$A \circ B = (A \ominus B) \oplus \hat{B} $$

This process in visualised in Figure \ref{fig:opening}, where the square
structural element is applied on the binary image. Morphological opening
preserves larger structural components, and removes small objects from the
image. This means, it is a good candidate for removing noise, small objects or
thin structures, which may interfere with later processing. \cite{soile2004}

\begin{figure}
    \begin{subfigure}[t]{0.4\textwidth}
        \centering
        \includegraphics[width=\textwidth]{resources/inkscape/opening_orig.png}
        \caption{}
        \label{fig:opening_orig}
    \end{subfigure}
    \begin{subfigure}[t]{0.4\textwidth}
        \centering
        \includegraphics[width=\textwidth]{resources/inkscape/opening_erosion.png}
        \caption{}
        \label{fig:opening_erosion}
    \end{subfigure}
    \begin{subfigure}[t]{0.2\textwidth}
        \centering
        \includegraphics[width=\textwidth]{resources/inkscape/opening_erosion_se.png}
        \caption{}
        \label{fig:opening_erosion_se}
    \end{subfigure}
    \begin{subfigure}[t]{0.4\textwidth}
        \centering
        \includegraphics[width=\textwidth]{resources/inkscape/opening_dilation.png}
        \caption{}
        \label{fig:opening_dilation}
    \end{subfigure}
    \begin{subfigure}[t]{0.2\textwidth}
        \centering
        \includegraphics[width=\textwidth]{resources/inkscape/opening_dilation_se.png}
        \caption{}
        \label{fig:opening_dilation_se}
    \end{subfigure}
    \caption{Morphological opening (\subref{fig:opening_dilation}) of the binary image (\subref{fig:opening_orig}) . Centres of structural elements are indicated by the circles. The gray pixels indicate which pixels
    were lost (erosion) or added (dilation) after the transformation. In (\subref{fig:opening_erosion}) the erosion of the binary image with a structuring element (\subref{fig:opening_erosion_se}) is shown, and in (\subref{fig:opening_dilation}) the dilation of the previous image by the structuring element (\subref{fig:opening_dilation_se}) is shown.}
    \label{fig:opening}
\end{figure}

\subsubsection{Closing}

The morphological closing is defined similarly to the opening, but the order of
operations is reversed. The closing $\bullet$ of a set $A$ by a structuring
element $B$ is defined as\cite{soile2004}:

$$A \bullet B = (A \oplus B) \ominus \hat{B} $$

\begin{figure}
    \begin{subfigure}[t]{0.4\textwidth}
        \centering
        \includegraphics[width=\textwidth]{resources/inkscape/closing_orig.png}
        \caption{}
        \label{fig:closing_orig}
    \end{subfigure}
    \begin{subfigure}[t]{0.4\textwidth}
        \centering
        \includegraphics[width=\textwidth]{resources/inkscape/closing_dilation.png}
        \caption{}
        \label{fig:closing_dilation}
    \end{subfigure}
    \begin{subfigure}[t]{0.2\textwidth}
        \centering
        \includegraphics[width=\textwidth]{resources/inkscape/opening_erosion_se.png}
        \caption{}
        \label{fig:closing_erosion_se}
    \end{subfigure}
    \begin{subfigure}[t]{0.4\textwidth}
        \centering
        \includegraphics[width=\textwidth]{resources/inkscape/closing_erosion.png}
        \caption{}
        \label{fig:closing_erosion}
    \end{subfigure}
    \begin{subfigure}[t]{0.2\textwidth}
        \centering
        \includegraphics[width=\textwidth]{resources/inkscape/opening_dilation_se.png}
        \caption{}
        \label{fig:closing_dilation_se}
    \end{subfigure}
    \caption{Morphological closing (\subref{fig:closing_erosion}) of the binary image (\subref{fig:closing_orig}). The circles indicate centres of structural elements. The grey pixels indicate which pixels were lost (erosion) or added (dilation) after the transformation. In (\subref{fig:closing_dilation}) the dilation with structuring element (\subref{fig:closing_erosion_se}) is applied on (\subref{fig:closing_orig}). In (\subref{fig:closing_erosion}) the erosion of (\subref{fig:closing_dilation}) with structuring element (\subref{fig:closing_dilation_se}) is shown.}
    \label{fig:closing}
\end{figure}

As for the previous operator, this is visualised in Figure \ref{fig:closing},
where the square structural element is applied on the binary image.
Morphological closing, like the opening, keeps the larger structural elements
but also fills in holes, which might have been created, for example, by noise.
Also, it does not remove small objects.

While mathematical morphology is a massive topic of interest to many
researchers, here I will only use the aforementioned operations. They can be
further used to define more operations such as performing top hat transforms,
skeletonizing images, finding perimeters etc \cite{mathworks_morp_oper}.

\section{Segmentation}

Image segmentation is a fundamental step in the image processing. Unlike in the
previous sections, we are no longer interested in manipulating and changing the
image. The point of image segmentation is to split the image into regions with
some specific meaning, e.g. in our case, to split the data into regions of
individual cells \cite{gonzalez2002}.

Numerous techniques can be applied, but in this thesis, we are interested in
one particular technique --- the watershed algorithm --- and segmentation based
on the Laplacian operator in the zero-crossing algorithm.

\subsection{Watershed algorithm}

The watershed algorithm is a widespread technique for image segmentation. The
fundamental principle of the watershed algorithm is to view the image as a
topographic map, where the intensities indicate the height of the landscape, in
other words, peaks in the image represent mountains or hills, while low values
represent valleys.

The algorithm then creates water springs in these valleys, from which the water
flows at a uniform rate. As the waters from different springs rise, they might
start to mix. At these points, we build dams to separate the various flooded
areas. The water continues to flow until just the dams are visible above the
water \cite{gonzalez2002}. At this point, the image is segmented using the dams,
and every catchment basin is uniquely labelled.

\begin{figure}
    \begin{subfigure}[t]{0.48\textwidth}
        \frame{\includegraphics[width=\textwidth]{./resources/watershed-basic.png}}
        \caption{}
        \label{fig:demo-watershed-orig}
    \end{subfigure}
    \begin{subfigure}[t]{0.48\textwidth}
        \frame{\includegraphics[width=\textwidth]{./resources/watershed-complex.png}}
        \caption{}
        \label{fig:demo-watershed-run}
    \end{subfigure}
    \caption{Demonstration of the watershed algorithm. In (\subref{fig:demo-watershed-orig}) a binary image with two overlapping circles is shown. Running the algorithm produces (\subref{fig:demo-watershed-run}).}
    \label{fig:demo-watershed}
\end{figure}

In Figure \ref{fig:demo-watershed}, the result of running the watershed algorithm 
on a simple image is shown.

\subsection{Distance transform}

One of the last required terms to define is the distance transform, which is
used in the implementation to calculate the watershed algorithm.

The distance transform is defined on binary images. After the distance transform
is applied on an image, each intensity is replaced by the closest distance to
the background. In principle any metric can be used, but the Euclidean distance
is used in our case. A simple example is shown in Figure
\ref{fig:demonstration-distance-map}.

\begin{figure}
    \begin{subfigure}[t]{0.48\textwidth}
        \frame{\includegraphics[width=\textwidth]{"./resources/distance-map-cirle.png"}}
    \end{subfigure}
    \begin{subfigure}[t]{0.48\textwidth}
        \frame{\includegraphics[width=\textwidth]{"./resources/distance-map-transform"}}
    \end{subfigure}
    \caption{Demonstration of the distance transform. The farther the pixel is
    from the background, the higher the pixel's intensity.}
    \label{fig:demonstration-distance-map}
\end{figure}

\clearpage

\chapter{Image data}

The image data --- of mouse mammary organoids --- was acquired using a confocal
laser scanning microscope. The original data was captured in 2021 by a group of
researchers at the Faculty of Medicine, Masaryk University. The captured field of view is around 500
microns by 500 microns by 80 microns, and a single voxel's size is $0.233\mu m$ by
$0.233\mu m$ by $0.586\mu m$.

In confocal laser scanning microscopy, researchers can stain specific structures in
their specimen so that the studied structures emit light when the proper wavelength
is shined upon them. This allows them to capture the expression of studied proteins
inside the specimen. In this particular case, the researchers captured the
expression of 3 different proteins, of which one is part of membranes, while the
other two are primarily located inside the cells' nuclei.

The data, totalling around 576 GB, consists of 12 different time steps, where
each step has a size of 1920x1920x388 voxels. For every protein, the researchers
captured the image using two different optical setups; thus, there are six
channels. The data is stored in numerous TIFF files, and every channel contains
grayscale data with 16 bits per voxel.
\begin{figure}
    \begin{subfigure}[t]{0.4\textwidth}
        \centering
        \includegraphics[width=\textwidth]{resources/C3-t006-200-scaled.jpg}
        \caption{}
        \label{fig:data_example_membraine}
    \end{subfigure}
    \begin{subfigure}[t]{0.4\textwidth}
        \centering
        \includegraphics[width=\textwidth]{resources/C2-t006-200-scaled.jpg}
        \caption{}
        \label{fig:data_example_nuclei}
    \end{subfigure}
    \caption{On the left panel a sample from the membrane-stained channel is shown, while on the right panel a corresponding slice from the nuclei-stained channel is shown.}
    \label{fig:data_example}
\end{figure}
In Figure \ref{fig:data_example}, a single slice of the image data was taken from
two channels. They are from the same time point and the same place.

\section{Typical challenges to deal with}

The size of the data and the inherent phenomenons of fluorescence microscopy
provide the main challenges in analysing the data. It is very time-consuming to
analyse the whole data set; therefore, all work was done on smaller data slices.

\begin{figure}
    \centering
    \begin{subfigure}[t]{0.5\textwidth}
        \centering
        \frame{\includegraphics[width=\textwidth]{./resources/bad-edges.png}}
        \caption{}
        \label{fig:bad_edges}
    \end{subfigure}
    \hfill
    \begin{subfigure}[t]{0.3\textwidth}
        \centering
        \frame{\includegraphics[width=\textwidth]{./resources/bleading-edge.png}}
        \caption{}
        \label{fig:bleeding_edges}
    \end{subfigure}
    \hfill
    \begin{subfigure}[t]{0.5\textwidth}
        \centering
        \frame{\includegraphics[width=\textwidth]{./resources/light-scattering.png}}
        \caption{}
        \label{fig:scattering}
    \end{subfigure}
    \caption{Typical issues with the data. In (\subref{fig:bad_edges}) many of the membranes are missing, in (\subref{fig:bleeding_edges}) there is not a clear boundary between the specimen and the background, and in (\subref{fig:scattering}) the scattering of light inside the organoid, causes large problems with the segmentation.}
\end{figure}

\subsection{Noise}

The first problem with the data set is the quality. In most places, the noise is
too strong to try to use in any analysis.

In the membrane-stained channel, two main problems can be seen. To begin, most
membranes do not form connected segments; instead, a part of the cell's membrane is
usually missing, or it is joined with the neighbouring cell's membrane, as can be
seen in Figure \ref{fig:bad_edges}. This complicates the analysis, as different
regions may require different filters to be properly preprocessed. Moreover,
this makes edge detection more difficult because it relies on the magnitude of
the gradient of the image. Moreover, sometimes the intensity change between the cell and the
background is too smooth for the program to find the divisive line as in Figure
\ref{fig:bleeding_edges}.

Another issue is that the light gets scattered inside the organoid, as shown
in Figure \ref{fig:scattering}. That means that it is possible to reliably analyse only the
outermost layers of cells of the organoid.

The nuclei-stained channel comes with its own set of problems, some of them different.
It is still true that the sharpest data comes, from the outermost layer
of the organoid, severely limiting the amount of usable data. Another issue is that the
individual cells are sometimes quite fuzzy, which makes it difficult to decide where
exactly the cell ends or if there are multiple cells instead of one.
% todo: reword

\subsection{Size}

The size of the dataset poses the second problem, as it would be computationally
expensive to try to process the whole dataset at once; however, a lot of the
data is either unusable or of no interest. In most instances, the image data is
predominantly composed of the background, and at several moments in time, no
specimen can be seen; therefore, to reduce the computational demands further,
only specific regions of interest were manually selected to create the final
dataset. Ultimately, the reduced dataset is only a few dozen megabytes large.

\section{Dataset}

For reasons specified in the previous sections, it would not be reasonable
to run the evaluation algorithm on the whole image data. Moreover, it is unclear
how prediction on the whole image data would be scored, not to mention the large
computing and storage requirements; thus, we decided to create a set of images
on which the algorithm would be evaluated. The set contains 11 images, which I
hand-picked from the raw image data. I tried to include sections from various
channels, times, and of various quality. This was done to test the
program's performance under diverse circumstances.

All of the images are included as attachments in the thesis. The general
name scheme is \texttt{txxx-roi-y-z}, where \texttt{xxx} stands for the time,
\texttt{roi-y} stands for the region of interest with index \texttt{y}, and the
\texttt{z} stands for the central frame. As the program can use
multiple frames from the membrane-stained channel, it is possible to provide
several slices in a single file. The program then runs the maximum projection
in 3D, but the slice in the middle will be picked for further processing. On the
other hand, the nuclei-stained data is always two-dimensional.

This dataset was then manually annotated as follows: Parts that were considered to be background were labelled using $0$.
Sections which were too complicated or unclear were labelled using $65535$, and
the rest of the cells were uniquely labelled with numbers starting \mbox{at $1$}. In
Figure \ref{fig:label_ground_truth} the manual segmentation of the corresponding raw image data (Fig \ref{fig:label_membrane_channel}) created by the expert is shown.

\begin{figure}
    \begin{subfigure}[t]{0.4\textwidth}
        \centering
        \frame{\includegraphics[width=\textwidth]{./resources/membrane-channel.png}}
        \caption{}
        \label{fig:label_membrane_channel}
    \end{subfigure}
    \begin{subfigure}[t]{0.4\textwidth}
        \centering
        \frame{\includegraphics[width=\textwidth]{./resources/ground-truth.png}}
        \caption{}
        \label{fig:label_ground_truth}
    \end{subfigure}
    \caption{In (\subref{fig:label_ground_truth}) the labelling created by an expert of the (\subref{fig:label_membrane_channel}) is shown. White means a complicated area, black is the background, and the rest are the cells.}
\end{figure}

\clearpage

\chapter{The developed segmentation method}
The goal of the algorithm is to segment individual cells and uniquely label each
cell; however, as previously discussed, they are often not clearly separated
from each other. A naive Otsu thresholding leads to unsatisfactory results, as is
shown in Figure \ref{fig:otsu-nuclei-demo}. Thus, a more sophisticated algorithm
had to be developed to segment the cells' nuclei finer.
\begin{figure}
    \begin{subfigure}[t]{0.45\textwidth}
        \centering
        \frame{\includegraphics[width=\textwidth]{resources/otsu-nuclei.png}}
        \caption{}
        \label{fig:otsu-nuclei}
    \end{subfigure}
    \hfill
    \begin{subfigure}[t]{0.45\textwidth}
        \centering
        \frame{\includegraphics[width=\textwidth]{resources/otsu-naive.png}}
        \caption{}
        \label{fig:otsu-naive}
    \end{subfigure}
    \caption{Otsu thresholding (\subref{fig:otsu-naive}) of the image (\subref{fig:otsu-nuclei}) does not lead to satisfactory segmentation.}
    \label{fig:otsu-nuclei-demo}
\end{figure}
We settled on a two-step algorithm. In the first step, the membrane-stained
channel is used, while in the second step, the nuclei-stained channel is used.
Every cell in our specimen has a membrane and nuclei; hence, we should be able
to use the information provided by the first channel to improve the thresholding
in the second channel.

In the first step, a coarse segmentation step is employed on the raw image data (Figure \ref{fig:max-proj-before}) to roughly find the membranes separating the individual cells (Figure \ref{fig:watershed}). The coarse segmentation is then utilized in the second step; the algorithm overlays the previous
segmentation on top of the nuclei-stained channel and looks at every label.
Each one represents a single cell and some of its surroundings. Finally, a
local version of the gradient thresholding is performed on each region (Figure \ref{fig:grad-thresh-post}).
\begin{figure}
    \begin{center}
        \includegraphics[width=0.8\linewidth]{./resources/inkscape/segmentation-steps.png}
    \end{center}
    \caption{All steps in the segmentation algorithm.}
    \label{fig:segmentation_steps}
\end{figure}

\section{Coarse cell segmentation}
The goal is to provide a coarse segmentation of the cells using the
membrane-stained channel.

\subsection{Preprocessing}
First, the input image data is pre-processed to reduce the amount of noise in it.
As shown in Figure \ref{fig:max-proj}, we need to deal with the
fact that the intensities inside every cell are quite uneven. To even out the
signal inside cells, a maximum projection (Section \ref{sec:max-proj}) is
performed. The idea behind this step is that cell membranes have higher intensities than
the cells' interiors; therefore, they should not be very influenced by this step,
while the uneven illumination inside cells would be smoothed out. That happens
because the highest intensity areas will propagate out to the surrounding area.

\begin{figure}
    \begin{subfigure}[t]{0.45\textwidth}
        \centering
        \frame{\includegraphics[width=\textwidth]{./resources/max_proj_before.png}}
        \caption{}
        \label{fig:max-proj-before}
    \end{subfigure}
    \begin{subfigure}[t]{0.45\linewidth}
        \centering
        \frame{\includegraphics[width=\textwidth]{./resources/max_proj_after.png}}
        \caption{}
        \label{fig:max-proj-after}
    \end{subfigure}
    \caption{The uneven illumination inside individual cells (\subref{fig:max-proj-before}) is mitigated using the
    maximum projection technique (\subref{fig:max-proj-after}).}
    \label{fig:max-proj}
\end{figure}

After dealing with the uneven lighting, additional smoothing is applied
to reduce the effect of residual noise further. Since the voxel's dimensions are
anisotropic, an anisotropic Gaussian smoothing is applied.

\subsection{Zero-crossing algorithm}
\label{section:zc}
Applying the zero-crossing algorithm at this stage would result in many
artificial cells being found in the background due to the presence of noise. To
remove these visual artefacts, we need to create a mask to differentiate the
background from the foreground. For this task, the Otsu thresholding method was
employed due to the bimodal distribution of the image intensities. Sometimes the interface between
the organoid's contour and the background is unclear, as depicted in Figure
\ref{fig:max-proj-before}. However, it is a considerably smaller issue if part of the
the background is categorised as the foreground; hence, this issue can be fixed by
running morphological erosion with a circular structural element with a diameter similar to the average cell diameter. This final mask then approximately distinguishes between the 
background and the foreground. This mask is then utilised in the next steps, which
ignore all background image elements during their calculations.

\begin{figure}
    \begin{subfigure}[t]{0.45\textwidth}
        \centering
        \frame{\includegraphics[width=\textwidth]{./resources/threshold-before.png}}
        \caption{}
        \label{fig:cell-mask-bef}
    \end{subfigure}
    \begin{subfigure}[t]{0.45\linewidth}
        \centering
        \frame{\includegraphics[width=\textwidth]{./resources/threshold-expanded.png}}
        \caption{}
        \label{fig:cell-mask-aft}
    \end{subfigure}
    \caption{The Otsu method might sometimes classify parts of the specimen as
    background (\subref{fig:cell-mask-bef}); thus, the mask must be expanded (\subref{fig:cell-mask-aft}) using morphological erosion.}
    \label{fig:cell-mask}
\end{figure}

The zero-crossing algorithm provides a reasonable outline of the cells, as is
shown in Figure \ref{fig:zc}. The user sometimes has to play around with
the strength of the Gaussian smoothing (parameter \texttt{zc\_sigma} described in later chapters) to improve the quality of the edge detection. When the strength is too low, the program tends to find many small artificial regions; conversely, strong smoothing causes the algorithm to join individual cells together. 
This step is further complicated by the fact that some of the cells are missing parts of
their membranes; hence, the performance of this step is one of the limiting factors of
the whole segmentation pipeline.

\begin{figure}
    \begin{subfigure}[t]{0.45\textwidth}
        \frame{\includegraphics[width=\textwidth]{./resources/zero-crossing.png}}
        \caption{}
        \label{fig:zc}
    \end{subfigure}
    \begin{subfigure}[t]{0.45\textwidth}
        \frame{\includegraphics[width=\textwidth]{./resources/mask.png}}
        \caption{}
        \label{fig:cells}
    \end{subfigure}
    \begin{subfigure}[t]{0.45\textwidth}
        \frame{\includegraphics[width=\textwidth]{./resources/mask-postprocessed.png}}
        \caption{}
        \label{fig:cells-post-processed}
    \end{subfigure}
    \begin{subfigure}[t]{0.45\textwidth}
        \frame{\includegraphics[width=\textwidth]{./resources/watershed.png}}
        \caption{}
        \label{fig:watershed}
    \end{subfigure}
    \caption{Steps taken from zero-crossing algorithm towards watershed. First, in (\subref{fig:zc}) the edges found by the zero-crossing algorithm are shown. Then in (\subref{fig:cells}) the mask of cells is displayed. Then the mask is post-processed (\subref{fig:cells-post-processed}) and finally the result after watershed is shown (\subref{fig:watershed}).}
    \label{fig:zc-watershed}
\end{figure}

\subsection{Watershed}
The next step is to utilize the outlines to segment individual cells. First, we need
to take the previous result and fill in the inside of the cells, as can be seen
in Figure \ref{fig:cells}. How this is achieved is explained in Section
\ref{sec-markers}. After this step is done, the image is further post-processed,
starting by using erosion to remove the leftover edges found by the zero-crossing, and then a morphological opening with
a disk-like structuring element is performed to remove small artificial cells caused
by noise. Finally, dilation is applied to restore some of the unwanted losses
caused by the morphological opening, leading to the result in Figure 
\ref{fig:cells-post-processed}.

Finally, the watershed algorithm is used to label the individual cells uniquely.
Since the way it is implemented, even if two or more cells are joined in the
previous step, the watershed algorithm may be able to divide them. The result
is shown in Figure \ref{fig:watershed}. A more in-depth discussion is in Section
\ref{sec-markers}.

\section{Fine nucleus segmentation}
Unlike in the membrane-stained channel, in the nuclei-stained channel the
pipeline is lot less complicated. The whole pre-processing consists of a single
Gaussian smoothing, that is done because the data is quite grainy, and it helps
the gradient thresholding to produce more reliable results.

Applying the Otsu thresholding method in the nuclei-stained channel leads to
poor results, as shown in Figure \ref{fig:otsu-naive}. The problem is that the
lighting is quite uneven. Some cells are significantly brighter than others;
thus, global thresholding is not a viable option.

However, we can use the previously obtained segmentation from membrane-stained
channel to improve the quality of the thresholding. By looking at just specific
areas denoted by the labels, the thresholding yields more suitable results. In
essence, a local version of some thresholding algorithm can be used. Applying it
this way, the resulting thresholding is shown in Figure \ref{fig:grad-thresh}.

In our testing, the Otsu thresholding method was not suitable as in many cases, the
proximity of cells shrinks the background area, causing the algorithm to
categorise some of the foreground as the background. After some experimentation,
we concluded that the gradient thresholding algorithm provides better results.
Also, others easily accessed thresholding methods in the \texttt{scikit-image}
library were experimented with, but the gradient thresholding method provided
reasonable results in many cases.

Even though the data was smoothed before applying the threshold, many edges in
the binary thresholded image appear quite jagged. To combat this, a final round
of post-processing is applied. It consists of applying binary morphological
closing on each label to fill in dark spots, thus, rounding off the jagged edges.
The overall effect of this step is not major, but it can help. It is shown in
Figure \ref{fig:grad-thresh-post}.

\begin{figure}
    \begin{subfigure}[t]{0.45\textwidth}
        \frame{\includegraphics[width=\textwidth]{./resources/nuclei.png}}
        \caption{}
        \label{fig:nuclei}
    \end{subfigure}
    \begin{subfigure}[t]{0.45\textwidth}
        \frame{\includegraphics[width=\textwidth]{./resources/thresholded-gradient.png}}
        \caption{}
        \label{fig:grad-thresh}
    \end{subfigure}
    \begin{subfigure}[t]{0.45\textwidth}
        \frame{\includegraphics[width=\textwidth]{./resources/thresholded-gradient-post.png}}
        \caption{}
        \label{fig:grad-thresh-post}
    \end{subfigure}
    \begin{subfigure}[t]{0.45\textwidth}
        \frame{\includegraphics[width=\textwidth]{./resources/ground-truth-crop.png}}
        \caption{}
        \label{fig:ground-truth-crop}
    \end{subfigure}
    \caption{Steps inside the nuclei-stained channel compared with the ground
    truth (\subref{fig:ground-truth-crop}). In (\subref{fig:nuclei}) an image from the nuclei-stained channel is shown. In (\subref{fig:grad-thresh}) the result of applying the gradient threshold is displayed. Then in (\subref{fig:grad-thresh-post}) the previous image was post-processed using morphological closing.}
\end{figure}

An example of the final segmentation and the reference segmentation --- created by a
biologist --- is shown in Figure \ref{fig:ground-truth-crop}.

\clearpage

\chapter{Implementation}
\label{chp:implementation}
This chapter describes which libraries were used, lists the main parameters, and
details non-trivial parts of the implementation.

\section{Justification for Python and used libraries}
As part of the bachelor thesis, I developed a terminal-based program to run the
aforementioned algorithm on new data, and it is available as an attachment to the
thesis. It is written in Python 3.11 and can be run after installing the
required packages. Those are specified inside the included
\texttt{environment.yml} file. The environment does not need any additional setup
apart from installing the required packages.

I decided to use Python as my main language. Since it is an interpreted language
with dynamic typing, it allows for the rapid development of new prototypes and easy
experimentation. It also has a very minimal and expressive syntax, which makes
it easy to read and write code \cite{python-docs-tutorial}. Finally, it has a
vast ecosystem of tools and libraries, accelerating further development. 

A significant disadvantage of the fact that Python is an interpreted language
is its, in general, poor performance. This can be mitigated to a certain degree
by using libraries instead of the default Python's data structures. Many of
these libraries use some compiled language --- such as C, C++ or Fortran --- in
the background to provide great performance gains while not breaking
the expressivity and simplicity of Python's syntax.

For these reasons, the scientific stack based on Python was chosen instead of
other alternatives, such as the ImageJ/Fiji or OpenCV ecosystem.

The implementation heavily relies on the following packages to do the calculations:
\begin{itemize}
    \item{NumPy \cite{harris2020array} is an indispensable package, allowing
        Python to be utilized as a tool in scientific computing, and a
        significant part of the whole scientific computing environment in Python
        is based around this package. In its centre is an \texttt{ndarray} object
        representing a multi-dimensional matrix with a uniform data type. The
        NumPy team has managed to utilize vectorization to improve the
        performance of many common operations on matrices greatly. This is achieved by
        writing a lot of the code in a low-level language such as C, which is
        pre-compiled and optimized, and taking advantage of the well optimized
        Intel MKL or OpenBLAS linear algebra library. This is all done in the
        background, mostly invisible to the programmer \cite{numpyManual2022}.}
    \item{SciPy \cite{2020SciPy-NMeth} provides a vast set of functions and
        implemented algorithms for scientific computing. It has a wide variety of
        uses, such as signal analysis, linear algebra solving, or calculating the
        Fourier transforms of images. In the application, this library was not used
        directly, but it is required by \texttt{scikit-image}.}
    \item{scikit-image \cite{scikit-image}, building on the SciPy's
        multidimensional image processing, provides numerous functions to
        process image data. The application relies on its implementation of
        Laplacian and gradient operators, mathematical morphology, and the watershed
        algorithm. As it is an extension of SciPy, it fits nicely into the
        technology stack.}
    \item{Python Imaging Library (Pillow) is a Python package for reading image
        data in various formats, efficient representation of the data, and even
        some manipulating of the image data \cite{clarkc20102023}. The sole
        purpose of this package for us was to load TIFF files into memory
        and to store them in a format supported by the other libraries.}
    \item{matplotlib \cite{hunter2007} is a popular library for creating
        visualisation in Python. It supports drawing all sorts of plots, but the
        program uses it to show results after some of the steps visually.}
    \item{conda is cross-platform a package management system and environment
        management system. It allows the user to quickly set up and use different
        Python environments and it was used during the development of this program.
        \cite{conda-manual}}
\end{itemize}
These libraries provide a great environment for writing well-performing image
processing pipelines fast. All of these libraries have large communities
behind them, are continuously developed, and are open-sourced.

\section{Main parameters}
The user can tweak several parameters of the segmentation pipeline. For a full description of all
the parameters and tips on changing them, refer to Appendix \ref{app:running-the-program}.

The four main parameters are:
\begin{itemize}
    \item \texttt{zc\_sigma} (float): This number indicates the strength of the Gaussian smoothing
    applied before the zero-crossing algorithm is applied.
    \item \texttt{maximum\_projection\_from} (integer): From which layer should the maximum projection
    be calculated
    \item \texttt{maximum\_projection\_to} (integer): Along with the previous parameter, specify which layers are used in the maximum projection.
    \item \texttt{cell\_radius} (integer): This number should equal an average radius of a single cell in
    the image data.
\end{itemize}

\section{Non-trivial parts of the implementation}
Most of the steps inside the algorithm map nicely to available functions in the
libraries, for example the \texttt{scikit-image} library provides a function for
calculating the Laplacian of the image; hence, this section provides an overview of
the non-trivial steps which I had to implement.

\subsection{Zero-crossing algorithm}

\begin{figure}
    \begin{center}
        \includegraphics[width=0.3\linewidth]{resources/neighbourhood.png}
    \end{center}
    \caption{The neighbourhood of the black pixel consists of the grey pixels in
    4-connectedness. In 8-connectedness, the light grey pixels are part of the
    neighbourhood too}.
    \label{fig:neighbourhood}
\end{figure}

The zero-crossing algorithm starts by calculating the Laplacian using a function
included inside the \texttt{scikit-image} library. To
extract the zero-crossings, the program will iterate over the whole array
(excluding the 1-pixel border) and check the following condition. Given its position $(i, j)$
the algorithm checks if $\nabla^2 I(i, j) < 0$ and its neighbours have
non-negative Laplacian.  If this condition stands, then there is
a zero-crossing at position $(i, j)$. In this case, the neighbourhood is
calculated as 4-connectedness. That means that the neighbourhood of the image element
consists of the image elements on top, left, bottom, and right of the image element, as is
shown in Figure \ref{fig:neighbourhood}.

\subsection{Markers for watershed algorithm}
\label{sec-markers}
Extraction of the markers for the watershed function is not as trivial as it
might seem. The algorithm tries to solve this issue by taking the result of the
zero-crossing algorithm and filling in the individual cells. To accomplish this
a flood-filling algorithm is used.

\begin{figure}
    \begin{subfigure}[t]{0.48\textwidth}
        \centering
        \includegraphics[width=\textwidth]{resources/zc.png}
        \caption{}
        \label{fig:fill-map-zc}
    \end{subfigure}
    \begin{subfigure}[t]{0.48\textwidth}
        \centering
        \includegraphics[width=\textwidth]{resources/flood\_fill.png}
        \caption{}
        \label{fig:fill-map-flood-fill}
    \end{subfigure}
    \begin{subfigure}[t]{0.48\textwidth}
        \centering
        \includegraphics[width=\textwidth]{resources/mask-dist.png}
        \caption{}
        \label{fig:fill-map-mask}
    \end{subfigure}
    \caption{Segmenting cells after zero-crossing algorithm. In
    (\subref{fig:fill-map-zc}) the result of the zero-crossing algorithm is
    shown. Then the flood filling mechanism is run
    (\subref{fig:fill-map-flood-fill}) and finally the mask separating cells and
    background is shown (\subref{fig:fill-map-mask}).}
    \label{fig:flood_fill}
\end{figure}

We obtain a binary image from the zero-crossing algorithm. Right now there are
only two values, $0$ stands for the background, while $255$ stands for the found
edges. We assume that the user usually puts the whole specimen into view;
therefore we assume that the borders of the image belong to the background.

The program iterates over the borders of the image. If the current position is
a background pixel, a flood fill is ran from there with a new value of $128$.
After all the iterations are done, we are left with a grayscale
image with three possible values, $0$ represents the original background, $128$
represents all the space, which was filled during the last step, and $256$
indicates the original edges. The program has then managed to separate the interiors of
the cells and the background as can be seen in Figure
\ref{fig:fill-map-flood-fill}. Now we can use the newly acquired labelling to
create a segmentation of cells, as is seen in the in Figure
\ref{fig:fill-map-mask}.

\begin{figure}
    \begin{subfigure}[t]{0.48\textwidth}
        \centering
        \includegraphics[width=\textwidth]{resources/post-processed-mask.png}
        \caption{}
        \label{fig:dist-processed-mask}
    \end{subfigure}
    \begin{subfigure}[t]{0.48\textwidth}
        \centering
        \includegraphics[width=\textwidth]{resources/distance-map.png}
        \caption{}
        \label{fig:dist-map}
    \end{subfigure}
    \begin{subfigure}[t]{0.48\textwidth}
        \centering
        \includegraphics[width=\textwidth]{resources/labeled\_regions.png}
        \caption{}
        \label{fig:dist-regions}
    \end{subfigure}
    \caption{In (\subref{fig:dist-processed-mask}) the post-processed mask is
    shown. Then the distance transform is applied (\subref{fig:dist-map}), and
    finally the watershed algorithm is run (\subref{fig:dist-regions}).}
    \label{fig:distance_map}
\end{figure}

Before this can be applied in the next steps, some of the smaller cells clearly
outside of the specimen has to be removed. These are caused by the imperfections
in the denoising parts. Even after the maximum projection and the cell mask are
applied, in many cases, some residual noise can create small fake cells. To
remove them, a morphological opening is performed, and the result is shown in
Figure \ref{fig:dist-processed-mask}.

Now that the cells were segmented, the distance map is calculated. This is shown
in Figure \ref{fig:dist-map}. The distance map was
calculated using the \texttt{SciPy} library, and it uses Euclidean distance.

Next, to extract the markers, the peaks inside the distance map were found using
the \texttt{peak\_local\_max} function from the \texttt{scikit-image} library. The benefit of
this function is that the user can provide a minimum distance separating two
markers. Additionally, the function allows the user to specify a threshold, which
means that peaks with low intensity will be ignored, further improving the
labelling.

Finally, the watershed algorithm is run on the inverse of the distance map, with the results
visualised in Figure \ref{fig:dist-regions}.

\subsection{Maximum projection} 
\label{sec:max-proj}
To calculate the maximum projection of the data, the program iteratively uses
Gaussian smoothing to create multiple images. The first image is left unchanged.
The second image is created by applying the Gaussian filter to the first. The
third image is created by applying the Gaussian filter to the previous image
again and so on, forming multiple layers on top of the original image. Then we can
iterate over the original image and vertically project the biggest image intensity.
The result of doing this at every position is the maximum projection.

That concludes all the non-trivial parts of the program; the rest were a direct
application of mathematical functions in the aforementioned libraries.

\clearpage
\chapter{Evaluation}
\label{chp:evaluation}
This chapter deals with evaluating the results, common problems, and the
performance of the pipeline. Two evaluation measures are defined and analysed on the part of
the data set.

\section{Scoring of the algorithm}

We created two key measures to evaluate the results. To test how many cells were
detected, we devised a detection quality score using the F1 measure, and to test the
quality of the segmentation, a test based on the Jaccard score was created.
\begin{figure}
    \begin{center}
        \includegraphics{resources/inkscape/evaluation.png}
    \end{center}
    \caption{Simplified model of cells. Blue cells were predicted by the program
    while the black ones are cells found by the expert.}
    \label{fig:evaluation_basic}
\end{figure}
To talk about the scoring, we need to define a way how to pair the real cell
with the predicted segmentation. A simplified model of the problem is in Figure
\ref{fig:evaluation_basic}. The blue cells in the figure are the program's
prediction, while the black cells are real, created by an expert. As can be seen,
the predictions on the left are quite good as they cover the majority of the real
cell, while in the top right part, the program did not manage to properly
segment the cell, and in the bottom right, the program missed completely,
creating a new fictional cell, which was not found by the expert.

To take the analysis further, we need a way how to pair the relevant predictions
with the real segmentations. In Figure \ref{fig:evaluation_basic}, I only used two colours
to denote the predicted and real cells; however, in the real segmentation, each
cell has a unique label. To complicate things further, similar cells do not
have to have the same label in the predicted and real segmentation.

Solving this problem was done using a simple measure based on set operations.
Looking at each cell as a set of its pixels, we can talk about the intersection
of two cells. To illustrate, take a look at the cells in the top left corner of
the aforementioned figure. If we consider the two cells as sets of their
pixels, their intersection is equal to the whole blue set, and the union is
equal to the black set.

In general, by picking a single real cell and one predicted cell, we can construct
two sets containing their pixels. The set containing the pixels of the real cell is
denoted as $R_i$, while the second set containing the pixels of the predicted
cell is denoted as $P_j$. If the intersection is large enough, these two
cells are considered to be a pair and are later used in the performance measure. In full, the condition is
$$|R_i \cap P_j| > \alpha \cdot |R_i| \Leftrightarrow R_i \text{ and } P_j \text{ form a pair}$$
for some threshold $\alpha$. By default, this threshold is set to 0.5, ensuring
that every $R_i$ can have just one pairing member $S_i$.

\subsection{Segmentation accuracy}
\begin{figure}
    \begin{center}
        \includegraphics{resources/inkscape/evaluation_imprecise_segmentation.png}
    \end{center}
    \caption{Segmentation where the program labelled a lot of extra areas.}
    \label{fig:evaluation_imprecise}
\end{figure}
The first measure evaluates the accuracy of the predicted segmentation coincides
with the expert's segmentation. Just checking the size of the intersection is
not good enough, as if the program just created large blobs as in Figure
\ref{fig:evaluation_imprecise}, they would technically contain the whole cell;
however, such a result is not satisfactory because the segmentation contains a
a lot of areas around the cell.

To take this into account, the segmentation accuracy is evaluated using the
Jaccard index. The Jaccard index of two sets $A$ and $B$ is defined
as\cite{2020eelbode}:
$$J(A, B) = \frac{|A \cap B|}{|A \cup B|}.$$
The range of the Jaccard index is from 0 to 1 (included). If the two sets are
identical, then the Jaccard index is equal to 1, while if the two sets do not
have any identical elements, the index is equal to 0.

The benefit of this index is that, unlike plain intersections, the size of the two
sets are taken into account; therefore, the improper segmentation in
\ref{fig:evaluation_imprecise} would cause the index to be quite low since the
the denominator would be much larger than the numerator.

In the final evaluation, after the pairing, the Jaccard index of each pair of
cells will be calculated and averaged, leading to a metric of how well the
program segmented the data. Sometimes, the program's segmentation is not good enough,
and a couple of the real cells are not paired with any predicted cells. We decided
to score these cells as having zero Jaccard index. To provide a more meaningful analysis,
two different averages are calculated. If these missed areas are counted, then it
is denoted using $\overline{J}$. Otherwise, the symbol $\overline{J_S}$ is used.

\subsection{Detection performance}
\begin{figure}
    \begin{center}
        \includegraphics{resources/inkscape/evaluation_with_TP.png}
    \end{center}
    \caption{A simplified model of the cells, including which are considered to
    be true positives, false negatives, and false positives.}
    \label{fig:evaluation_with_TP}
\end{figure}

Apart from the quality of the segmentation, a useful measure is how many cells
the program was able to correctly find compared to how many were found by the expert. The
concept of pairs from the previous section is also used here to pair up the
predictions with the real cells. If a pair is found, then it is considered as a
true positive. If a cell found by an expert does not have a partner, then it is
considered to be a false negative, and if a cell found by the program does have
a corresponding cell found by an expert, it is considered to be a false
positive. In Figure \ref{fig:evaluation_with_TP}, this idea is visualised. To denote the number of
true positives TP is used. Similarly for false positives it is FP and for false
negatives it is FN.

Now that the notion of true positives, false positives, and false negatives have
been mapped to our results, the standard F1 score can be applied. It is defined as
the harmonic mean between recall and precision \cite{sklearn-f1score}, formally:
$$\text{recall} = \frac{\text{TP}}{\text{TP} + \text{FN}}$$
$$\text{precision} = \frac{\text{TP}}{\text{TP} + \text{FP}}$$
$$F_1 = \frac{2 \cdot \text{recall} \cdot \text{precision}}{\text{recall} +
\text{precision}}$$
Generally speaking, precision is a useful measure to see how many of the predicted
cells are actual cells. If our system has a low recall, that means that a
significant number of cells found by the expert were missed by our program. On
the other hand, precision talks about how many of the program's predictions are actual
cells. Then the F1 score is just a useful overall measure, accounting for
under-segmentation and over-segmentation in our data.

\subsection{Results}
When evaluating the program on the dataset, I took each image and gave myself around 2--5
minutes to find the best segmentation. Before looking at the $\overline{J_s}$ scores, I
visually inspected the results. Once I was satisfied and could not see any of the common
issues, I tweaked the parameters to find the biggest $\overline{J_s}$ score.
After that, I recorded the results and proceeded to the next image.

The final results are reported in Table \ref{tab:scores} and the used
configurations in Table \ref{table:configurations}. Every row begins with the
filename and provides various scores I achieved with the program. As can be seen
from the table, certain images got better results than others.

In general, the program has slightly higher recall than precision, meaning that the
program finds the majority of the cells marked by the expert; however, it also
finds many false positives. Those are cells that the program predicted but were
not found by the expert. Those are mostly caused by noise in the data.

Looking at the $\overline{J_S}$ scores, it is clear that if the program finds a
real cell, it segments it reasonably well. In the last row of Table \ref{tab:scores}, a
simple average of all the metrics are shown with the standard deviation. In
Figure \ref{fig:comparison-1} and in Figure \ref{fig:comparison-2} the
comparison between the program's prediction and the expert's labelling is shown.
In case of the second figure, the program has a lower $\bar{J}$ because several
of the cells were not found by the program, due to significant noise in the
membrane-stained channel. On the other hand, in Figure \ref{fig:comparison-1}
the effect of noise was more suppressed, and the program found more of the
cells, leading to higher $\bar{J}$.

\begin{figure}
    \begin{subfigure}[t]{0.45\textwidth}
        \centering
        \includegraphics[width=\textwidth]{./resources/t004-roi-10-375-predtiction.png}
        \caption{}
        \label{fig:comparison-pred-1}
    \end{subfigure}
    \begin{subfigure}[t]{0.45\textwidth}
        \centering
        \includegraphics[width=\textwidth]{./resources/t004-roi-10-375-truth.png}
        \caption{}
        \label{fig:comparison-truth-1}
    \end{subfigure}
    \caption{Prediction (\subref{fig:comparison-pred-1}) for the file
    \texttt{t004-roi-10} and the expert's labelling
    (\subref{fig:comparison-truth-1})}
    \label{fig:comparison-1}
\end{figure}
\begin{figure}
    \begin{subfigure}[t]{0.45\textwidth}
        \centering
        \includegraphics[width=\textwidth]{./resources/t008-roi-6-180-prediction.png}
        \caption{}
        \label{fig:comparison-pred-2}
    \end{subfigure}
    \begin{subfigure}[t]{0.45\textwidth}
        \centering
        \includegraphics[width=\textwidth]{./resources/t008-roi-6-180-truth.png}
        \caption{}
        \label{fig:comparison-truth-2}
    \end{subfigure}
    \caption{Prediction (\subref{fig:comparison-pred-2}) for the file
    \texttt{t004-roi-6} and the expert's labelling
    (\subref{fig:comparison-truth-2})}
    \label{fig:comparison-2}
\end{figure}

\begin{table}
    \hspace*{-0.5cm}\begin{tabular}{|| c|c|c|c|c|c|c ||}
        \toprule
        File name & $\overline{J}$ & $\overline{J_S}$ & \ensuremath{F_1} & Precision & Recall & Avg. time [s]\\
        \midrule
        t002-roi-1&0.36&0.57&0.54&0.50&0.64&4.93\\
        t002-roi-2&0.39&0.54&0.67&0.63&0.71&3.83\\
        t004-roi-4-68&0.39&0.61&0.44&0.33&0.64&7.11\\
        t004-roi-4-85&0.66&0.66&0.38&0.23&1.0&8.69\\
        t004-roi-6&0.33&0.61&0.58&0.6&0.54&2.49\\
        t004-roi-10&0.40&0.61&0.73&0.83&0.65&9.22\\
        t004-roi-11&0.43&0.60&0.63&0.6&0.72&4.89\\
        t006-roi-2&0.38&0.62&0.58&0.54&0.62&7.56\\
        t006-roi-5&0.29&0.63&0.42&0.39&0.46&6.27\\
        t008-roi-5&0.32&0.59&0.54&0.53&0.55&8.65\\
        t008-roi-6&0.20&0.54&0.41&0.44&0.38&3.78\\
        \midrule
        Averages&$0.38$&$0.61$&$0.54$&$0.51$&$0.63$&$6.13$\\
        Standard dev.&$0.11$&$0.037$&$0.12$&$0.16$&$0.16$&$2.29$\\
        \bottomrule
    \end{tabular}
    \caption{Evaluation of the files inside the data set.}
    \label{tab:scores}
\end{table}

\begin{table}
    \hspace*{-0.5cm}\begin{tabular}{|| c|c|c|c|c ||}
        \toprule
        File name & \texttt{max\_proj\_from} & \texttt{max\_proj\_to}& \texttt{zc\_sigma} & \texttt{cell\_size}\\
        \midrule
        t002-roi-1   &0&2&2.0&25\\
        t002-roi-2   &0&1&2.6&25\\
        t004-roi-4-68&0&2&3.0&25\\
        t004-roi-4-85&0&2&3.0&25\\
        t004-roi-6   &0&2&3.0&20\\
        t004-roi-10  &0&3&2.4&25\\
        t004-roi-11  &0&2&3.0&25\\
        t006-roi-2   &0&2&3.0&30\\
        t006-roi-5   &0&2&4.0&30\\
        t008-roi-5   &0&1&2.2&25\\
        t008-roi-6   &0&2&4.5&25\\
        \bottomrule
    \end{tabular}
    \caption{Segmentation settings for the files.}
    \label{table:configurations}
\end{table}

The last column in Table \ref{tab:scores} indicates the average time in seconds taken to run the program --- on my computer with the Ryzen 5700U notebook processor --- on
particular data. This was done using the \texttt{time} module from the standard
Python library and only the time taken to calculate was considered,
so the disk speed did not significantly affect the time. This was done over
four runs and averaged.

\section{Common problems}

Using the program on the dataset revealed two main typical issues with this
approach. As mentioned in the previous sections, many cells do not have a clear boundary
between them causing the program to segment them as one, or completely miss them.
This is evident in Figure \ref{fig:issues-joining}, where a lot of noise
causes the program to completely miss some of the cells. This can sometimes be
alleviated by lowering the \texttt{zc\_sigma} parameter, but overall, the method
fails in this area.

\begin{figure}
    \begin{subfigure}[t]{0.3\linewidth}
        \frame{\includegraphics[width=\textwidth]{./resources/issues-data.png}}
        \caption{}
        \label{fig:issues-joining-missing}
    \end{subfigure}
    \begin{subfigure}[t]{0.3\linewidth}
        \frame{\includegraphics[width=\textwidth]{./resources/issues-complicated-areas-filled.png}}
        \caption{}
        \label{fig:issues-joining-mask}
    \end{subfigure}
    \caption{A complicated area in the upper part of the image data (\subref{fig:issues-joining-missing}) causes the program to join multiple cells together in the segmentation (\subref{fig:issues-joining-mask}).}
    \label{fig:issues-joining}
\end{figure}

The second common issue happened during the marker creation for the watershed
algorithm. Here the biggest problem was that our assumption that the image's
outline was part of the background, causes the program to fill in the whole
specimen. This is shown in \ref{fig:issues-cropping}. This is because, in this
particular case, the division between the foreground and background was easily
visible and caused the program to create a connected contour around the
specimen. In this case, the solution is quite easy; after manually cropping the
bottom part
of the image, the program is able to properly segment the image.

\begin{figure}
    \begin{subfigure}[t]{0.3\textwidth}
        \frame{\includegraphics[width=\textwidth]{./resources/mask-failed.png}}
        \caption{Image.}
    \end{subfigure}
    \begin{subfigure}[t]{0.3\textwidth}
        \frame{\includegraphics[width=\textwidth]{./resources/issues-mask-croped.png}}
        \caption{Cropped image.}
    \end{subfigure}
    \caption{By cropping the image, the program was able to find better segmentation.}
    \label{fig:issues-cropping}
\end{figure}

\clearpage
\chapter{Conclusion}

The primary goal of this thesis was to develop a pipeline for the automatic
segmentation of membrane-stained and nuclei-stained image data from fluorescence microscopy.
Creating such tools makes it easier for researchers to analyse the behaviour of
cells by decreasing the time required to segment them manually. Due to the
unavailability of annotated data, techniques based on machine learning could not
be used; therefore, we settled on solving the segmentation using traditional
image processing techniques. The developed pipeline consists of two main steps.
In the first step, the Laplacian operator and zero-crossing algorithm is used to
find a coarse segmentation of the cells, while in the second step, this is
further refined using thresholding to find the nuclei-stained cells.

The developed program is fully functional and available as an attachment to the
thesis alongside the whole data set, including the ground truth created by the
expert and the already prepared membrane-stained and nuclei-stained data. The
configuration file is set up to reasonable defaults, which should provide fair
segmentation on most of the provided images.

The program provides a practical, easy-to-use terminal interface. In Appendix
\ref{app:running-the-program}, all the available parameters are documented, with tips
about which one to change in case of improper segmentation. The included zip
archive allows the user to test the program out of the box.

Lastly, the performance of the program and the quality of its segmentation has
been evaluated in Chapter \ref{chp:evaluation} using the F1 score and Jaccard
index. The program finds, on average, about half of the actual cells. Clearly
visible cells are nicely segmented, but the program often fails when some edges
are missing or there is a substantial unevenness in the signal heterogeneity. 

In the future, this work could be improved using additional steps in the
preprocessing stage. From our research into the topic, machine learning
techniques could be utilised to enhance the missing edges, or approximate
deconvolution algorithms might help with the denoising. However, those
techniques might take a considerable amount of reference data. Last but not
least, although the performance is reasonable, if the program was to be used for
the entire data, a just-in-time compiler like Numba or an ahead-of-time compiler
like Pythran could take the performance much further.

\printbibliography[heading=bibintoc] %% Print the bibliography.

\appendix %% Start the appendices.
\chapter{Running the program}
\label{app:running-the-program}
The program has two main modes of operation; those are segmentation and
evaluation. The first mode deals with loading the data and running the
segmentation algorithm, while the latter is used to score the predicted
segmentation using ground truth data. The user can choose the preferred
mode by running the program with the \texttt{-s} flag for segmentation and the
\texttt{-e} flag for the latter. Additionally, a helpful overview is
provided when the help flag \texttt{-h} is passed to the program.

The program starts by loading the configuration file \texttt{config.py}, inside
which is a dictionary object called \texttt{config}, and checks if all the types are
correct. This dictionary stores the
relevant settings to the segmentation and evaluation algorithm. The main parameters, which the
user might need to tweak to get better segmentation, are:
\begin{itemize}
    \item \texttt{zc\_sigma} (float): This number indicates the strength of the
        Gaussian smoothing before the zero-crossing algorithm is applied. When this
        number is too low, the program tends to find too many cells. If this number is 
        too large, then the program usually joins multiple cells together. The typical
        range is from 1.5 to 3.5.
    \item \texttt{maximum\_projection\_from} (integer): From which layer should
        the maximum projection be calculated. If you want to add more details, be sure
        to set this to 0 otherwise finer structures are lost during the maximum
        projection.
    \item \texttt{maximum\_projection\_to} (integer): Along with the previous
        parameter, specify which layers are used in the maximum projection. If the
        program creates several fictional cells inside a single cell, or is performing
        quite poorly in noisy areas, try increasing this number. Typically, there
        is not need to go above 4, and 2 is usually enough.
    \item \texttt{cell\_radius} (integer): This number should equal a roughly
        approximated radius of a single cell. This is mainly used during the
        thresholding in the membrane channel. If part of the specimen is
        classified as a background, try increasing the number.
    \item \texttt{axis\_relative\_sizes} (list of three floats): Specifies the
        ratios between dimensions of a single voxel.
    \item \texttt{show\_steps} (boolean): If this is equal to \texttt{True},
        the program will also plot partial results as the pipeline runs.
    \item \texttt{verbose} (boolean): If this flag equals True, the program
        will be more verbose during its runtime.
    \item \texttt{*\_folder} (string): These settings specify the relative paths
        to folders where data is loaded and stored.
\end{itemize}
The program also has a couple of extra tweakable parameters but in my experience,
they did not change the results significantly; therefore, I do not recommend
changing them in the first steps. Those parameters are:
\begin{itemize}
    \item \texttt{zc\_min\_distance} (float): This is passed to the
        \texttt{peak\_local\_max} function from \texttt{scikit-image} package.
        It changes how close can two markers in the watershed be. If you feel
        that the program is over-segmenting, increasing this number might help;
        however, changing the sigmas is a lot more significant.
    \item \texttt{pairs\_threshold} (float): This is used in the evaluation for
        determining if two sets can be considered pairs. 1 means that the sets
        must be identical, while 0.5 means that they must share at least half of
        the pixels.
    \item \texttt{nuclei\_sigma} (float): This parameter specifies the strength
        of the Gaussian smoothing in the nuclei-stained channel. Increasing this
        parameter can help with noise since the gradient threshold is based on the
        derivative of the image, but overall its effect is quite negligible on the
        final result.
    \item \texttt{segmentation\_postprocess\_radius} (integer): The radius of the
        disk-like structural element used in the post-processing of the
        segmentation. This can help to eliminate small regions which were
        incorrectly found due to noise. If the post-processing is deleting some of 
        the cells, decrease this number gradually. On the other hand, if the segmentation
        requires cleanup from many small areas, try increasing this number. Normally,
        this is a much smaller number than the \texttt{cell\_radius}.
\end{itemize}

When the user selects the segmentation option, the program looks at all TIFF
files inside folders \texttt{membrane\_data} and \texttt{nuclei\_data}, if two files have the
the same name, then the program pairs them together. Then it starts going over all
the pairs and running the segmentation algorithm, storing the resulting
segmentations inside the \texttt{segmentation\_output} folder.

If the user selects the evaluation option, the program looks into two different
folders --- \texttt{ground\_truth} and \texttt{prediction} --- pairs files with the
the same name and runs the evaluation algorithms, which are described in more
details in Chapter \ref{chp:evaluation}.

The code is split into the following files: 
\begin{itemize}
    \item \texttt{main.py}: The main application is located in the
        here; it parses the flags passed by the user, and loads
        up the configuration file. After that is done, the application starts
        either the segmentation or evaluation pipeline.
    \item \texttt{config.py}: This is the configuration file. It mainly includes
        the \texttt{config} dictionary, containing all the settings. It also
        checks the types of values inside the dictionary.
    \item \texttt{evaluation.py}: Here, all functions related to evaluating
        lies, e.g., finding pairs or calculating the Jaccard index.
    \item \texttt{filters.py}: This file contains implementations of various
        filters that the already mentioned libraries do not have implemented
        or which we decided to implement ourselves.
    \item \texttt{segmentation.py}: The code in this file is responsible for
        segmenting the raw data. It contains all the steps that have to be
        run, except for the few which are located in the \texttt{filters.py}
        file.
    \item \texttt{environment.yml}: Inside this file, all the required packages
        are specified, along with the Python version and which conda repository 
        channel to use while installing plugins.
    \item \texttt{tiff\_utils.py}: Miscellaneous code used
        to load the image data.
\end{itemize}


\end{document}