\documentclass[a4paper,oneside]{saip}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% defaults, keep as is, don't add anything that changes the document's defaults.
\usepackage{helvet}          % to have better matching fonts between Latex and Word (Title/Headings)
\usepackage{times}           % to have better matching fonts between Latex and Word (bulk text)
\usepackage{courierten}   % to have better matching fonts between Latex and Word (monospaced)
\usepackage[T1]{fontenc} % needed to deal with these fonts

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% common packages needed, add or comment as needed
\usepackage{graphicx}  %to include figures
\graphicspath{ {./img/} } 
\usepackage{url}
\usepackage{geometry} %	need: 
\geometry{%showframe, 
top=40.0mm, bottom=27.0mm, left=25.0mm, right=25.0mm}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% optional packages, you can add others as well for your own needs
%\usepackage{xcolor}
%\usepackage{algorithm2e}
%\usepackage[version=4]{mhchem} %for isotopes and all with \ce{^{15}_{8}O_{7}}
%\usepackage{amsmath}
%\usepackage{amssymb}
%\usepackage{stmaryrd}
%\usepackage{mathrsfs}
%\usepackage{mathabx}
%\usepackage{upgreek}
%\usepackage{comment}
%\usepackage{verbatim}
%\usepackage{}
\usepackage{siunitx}
\DeclareSIUnit\barn{b}
\usepackage{booktabs}
\usepackage{svg}
\usepackage{mathtools}
% \usepackage{lmodern}  % Enable modern font handling
\usepackage[noadjust]{cite}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DO NOT USE the hyperref package, your PDF will be included in a bigger document and the links will be broken.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}


\title{Deep Learning for High Throughput Decision Making on Diamond Content of PET Activated Kimberlite Rocks}

\author{Nicholas Connell$^{1}$, Timothy Brooks$^{2}$, Tiego Mpai$^{2}$, Thendo Nemakhavhani$^{2}$, Simon Connell$^{2}$ and Charis Harley$^{1}$}

\affil{$^1$Department of Electrical and Electronic Engineering Science, Faculty of Engineering and the Built Environment, University of Johannesburg, Johannesburg, South Africa}
\affil{$^2$Department of Mechanical Engineering Science, Faculty of Engineering and the Built Environment, University of Johannesburg, Johannesburg, South Africa}

%only put the email of the main corresponding author, not of all authors.
\email{\url{nickpconnell@gmail.com}}



\begin{abstract}
The application of computer-vision and deep-learning for the sorting of crushed rock in diamond mines is investigated. To this end vast quantities of positron emission tomographic images are synthesised to train convolutional neural networks as classifiers.
%Traditional diamond extraction methods employ multi-stage rock-crushing where stages are chosen to balance a mine's diamond-size distribution against the crush size at each stage so that an overall acceptable rate of diamond breakage is maintained. Positron emission tomography (PET) techniques can be used to reconstruct a density map of the distribution of PET isotopes within rocks. PET isotopes in diamond ores can be activated with a high energy gamma ray beam, and can penetrate to a depth of more than 30 cm. An upper bound of 10 to 15 cm on the crush size is given by the attenuation of 511 keV gammas in the rock. This decreases the rate of diamond breakage. Autonomous agents capable of on-the-fly decision making, these agents can then identify which rocks contain diamonds and extract them for further careful processing. The expected mass ratio of kimberlite to diamond is about a few billion to 1, thus a fairly accurate agent can reduce processing needs by at least 1000 times, we aim to reduce it by 10,000 times.
\end{abstract}


\section{Introduction}
The goal of the MinPET project is to develop, test, and validate a machine that can sort crushed rocks from a diamond mine. The machine is expected to be capable of sorting at least 500 tonnes of rock per hour. The sheer volume of rock that needs to be processed necessitates the use of autonomous AI agents. This rest of this paper describes the combination of synthetic data generation and deep-learning needed to create such agents and highlights the achievements of the MinPET research team 
\cite{ChinakaEM2016RSAa, CookMNH2016UGtc, cook_genetic, cook_overview, cook_thesis}.
% \cite{ChinakaEM2016RSAa}, \cite{CookMNH2016UGtc}, \cite{cook_genetic}, \cite{cook_overview}, \cite{cook_thesis}. 
Figure \ref{fig:run-of-mine-setup} depicts a hypothetical run-of-mine set up for the MinPET system. 
    
\begin{figure}
    \centering
    \includegraphics[width=0.5\linewidth]{img/minpetfig}
    \caption{Proposed run-of-mine MinPET system showing all critical functions of the system: activation, holding, detection, reconstruction, classification, and ejection}
    \label{fig:run-of-mine-setup}
\end{figure}

The scanning technique employed is positron emission tomography (PET). The rock is activated with \qty{40}{\MeV} gamma ray beam to induce the giant dipole resonance and transmute certain atoms into PET isotopes, namely carbon to carbon-11 so diamond content can be ascertained. This paper aims to show that the viability of training a convolutional neural network classify PET images of kimberlite rocks as containing diamond or barren.

%The MinPET machine offers three key benefits over traditional rock processing techniques: reduced diamond breakage, lowered costs to extract diamonds, and self-auditing. It lowers the chances of diamond being broken by stopping the crushing process earlier than traditional methods when the rock is approximately \qty{10}{\cm} in diameter. It lowers the extraction cost of diamonds because it only allows diamond containing rocks to proceed to the diamond extraction stage, where as traditional methods had to process all the rocks to find the diamonds. Lastly, the MinPET machine will be able to help audit the quantity and size of diamond contained in the rocks as the 3D PET image contains all of this information.

%%%%%%%%%%%%%%%% Problems with Diamond Mining %%%%%%%%%%%%%%%%%%%%
Traditional diamond mining sorting methods employ staged crushing, between which any revealed diamonds are retrieved. This requires vast amounts of energy and water, which is expensive and wasteful because deposits contain on average only a few carats of diamond per \qty{10000}{\kg} of rocks. Furthermore, it forces compromises to be made between profitability and acceptable diamond breakage; in mines known to produce diamonds larger than 100 carats these compromises often lead to the diamonds being broken into smaller less profitable pieces \cite{ANDRES201031}. Mines are often situated in energy insecure and water scarce countries where they negatively impact the local environment and put great stress on a host country's resources.

%%%%%%%%%%%%%%%% Benefits of Sensors %%%%%%%%%%%%%%%%%%%%
MinPET offers three key benefits over traditional rock processing techniques; reduced diamond breakage, lowered costs to extract diamonds, and lowered environmental impact by identifying and isolating diamond bearing rocks at an earlier stage in the crushing process. Immediate benefits include higher profits, lower energy costs, and lowered chance of damaging diamonds. The quantity of material requiring further processing will be greatly lowered because the majority of rocks are barren and will be discarded. The will make viable methods that are more environmentally friendly and are better at preserving diamond integrity \cite{ANDRES201031}. Additionally, the MinPET process will help audit the quantity and size of diamonds contained in the rocks as the PET image contains all of this information.

\section{Methodology}

The study consisted of four stages. First, the generation of synthetic data by Monte Carlo simulations using the Geant4 \cite{geant1} software. Second, the reconstruction of 3D images of the distribution of PET-isotopes using the Astra Toolbox Python library \cite{Astra1} and maximum likelihood expectation maximisation (MLEM). Third, training data was extracted from the images using custom Python code. Fourth, the AI was trained using an adapted version of a 3D ConvNet from the second place solution \cite{kaggle_2017} to the 2017 Kaggle competition written in Keras a Python library for machine learning.

\subsection{Simulation}

In total 3000 simulations were conducted to produce the training data and 3010 simulation were conducted to produce the test data. Each simulation contained one spherical kimberlite rock which housed one spherical diamond. The rock always had a diameter of \qty{100}{\mm}. The composition by mass fraction of the compounds of kimberlite used in the simulations is detailed in Table \ref{tab:kimb_comp}, which were determined for southern Africa \cite{kimb_comp, Tappe_Kimberlite}. For each training data simulation the diameter of the diamond in mm was sampled uniformly from the range $(0,15]$. The testing data was split into 10 groups each with different constant diamond diameter: \qty{1}{\mm}, ..., \qty{10}{\mm}. The diamond's location in every simulation was randomised by sampling each coordinate uniformly and resampled if the position placed any amount of the diamond outside of the rock.
    
\begin{table}[ht]
    \centering
        \label{tab:kimb_comp}
    \begin{tabular}{@{}c|cccccccccccc@{}}
         \toprule
         Compound & Si$_4$O$_{2}$ & MgO & CaO & CO$_2$ & H$_2$O & Fe$_2$O$_3$ & FeO\\
         Mass Fraction (\%) & 42.5562 & 28.6058 & 5.3233 & 5 & 5 & 4.6938 & 3.3138\\
         \midrule
         \midrule
         Compound & Al$_2$O$_3$ & TiO$_2$ & K$_2$O & P$_2$O$_5$ & MnO & Na$_2$O\\
         Mass Fraction (\%) & 2.5158 & 1.7979 & 0.6799 & 0.1914 & 0.1889 & 0.1309\\
         \bottomrule
    \end{tabular}
    \caption{Kimberlite Composition by Mass Fraction.}
\end{table}

The activity of the kimberlite and the diamond was determined by assuming a uniform electron flux of \qty{2.8086e13}{\per\centi\meter\squared} incident on a copper bremsstrahlung target, and incorporating a delay of 23 minutes to simulate the time spent in the hold hopper after irradiation and before detection. Table \ref{tab:pet_isotopes} summarises the possible pathways to produce PET isotopes, their half life in seconds, and the effective reaction cross-sections for a photon beam produced by an electron beam of \qty{40}{\mega\eV}.

\begin{table}[ht]
    \centering
        \label{tab:pet_isotopes}
    \begin{tabular}{@{}c|cccccc@{}}
        \toprule
        Reaction & $^{12}$C $\to$ $^{11}$C & $^{16}$O $\to$ $^{11}$C & $^{16}$O $\to$ $^{15}$O & $^{40}$Ca $\to$ $^{38}$K & $^{54}$Fe $\to$ $^{53}$Fe & $^{46}$Ti $\to$ $^{45}$Ti\\
        % \midrule
        Cross-Sec. (\unit{\milli\barn}) & 0.18791177 & 0.001010126 & 0.18038939 & 0.1121721 & 1.4198512 & 1.60765\\
        Half life (\unit{\second}) & 1220.04 & 1220.04 & 122.24 & 458.40 & 510.60 & 11088.00\\
        \bottomrule
    \end{tabular}
    \caption{PET-isotopes and their parent atoms with reaction cross-section at \qty{40}{\mega\eV} and isotope half life.}
\end{table}

The rock and the diamond are simulated passing through the detectors at \qty{1}{\meter\per\second}. There were two detector planes situated \qty{150}{mm} above and \qty{150}{mm} below the path of the rock. The detectors were rotated by \qty{1}{\degree} around the horizontal and vertical axes to take advantage of digital oversampling. The pixels of the detector were Bismuth Germanate (Bi$_4$Ge$_3$O$_{12}$). Pixels are grouped into photo-multiplier tubes (PMT), and PMTs were grouped into MinPET Units (MPU). Table \ref{tab:detector_measurements} summarises the geometry of the detectors' components by their side length, the gap between consecutive elements, and the grid they form. Each detector array was approximately \qty{1.6}{\meter} long in the rock's direction of travel and \qty{1.2}{\meter} wide.

\begin{table}[ht]
    \centering
        \label{tab:detector_measurements}
    \begin{tabular}{@{}cccc@{}}
         \toprule
         Component & Length (\unit{\mm}) & Gap (\unit{\mm}) & Grid (XxY)\\
         \midrule
         Pixel & 2.84 & 0.2 & 16x16\\
         PMT & 48.44 & 2.0 & 4x4\\
         MPU & 199.76 & 5.0 & 8x6\\
         \bottomrule
    \end{tabular}
    \caption{Detector component definitions.}
\end{table}

The output of simulations was a collection of coincident hit pairs. For two gamma ray hits to be considered a coincident they needed to come from different detectors, the difference in their hit times needed to be \qty{2}{\nano\second} or less, and there needed to be no other gamma rays within \qty{2}{\nano\second} of either of their hit times. This meant the final output file consisted of nine data points, two energy deposits by each hit, two sets of 3D coordinates, and one hit time.

\subsection{Image Reconstruction}

Each coincident hits pair can be used to calculate a line of response that connects the two hit locations. It can be assumed with some confidence that a positron will have annihilated somewhere along this line. Given that most positrons annihilate with an electron close to their parent atom one can begin to reconstruct the most likely distribution of PET isotopes in the rock given many lines of response.

First, a global translation was applied to all hit locations such that the reference changed from the detector frame to the moving rock frame. Equation \ref{eq:motion_correction} show how this was done, relying on three quantities: the velocity of the rock/belt $v_{belt}$, the hit time $ t_{hit}$, and the time at which the rock was centred between the detectors $t_{centre}$.

\begin{equation}
    \label{eq:motion_correction}
   x_{new} = x_{old} + v_{belt} (t_{centre} - t_{hit})
\end{equation}

Second, the line of response for each coincident pair was determined. The intersection of this line with the mid plane was computed, where the mid plane is the horizontal plane mid-way between the detectors. The angle that the line made with the mid plane was also computed in spherical coordinates. This resulted in four data points for each pair: the polar (measured from the vertical) and azimuthal angles, and the x and y intersection coordinates. Any lines of response with a polar angle greater than $\pi/4$ radians were discarded.

Third, the four dimensional data computed from the coincident pairs were binned to get a 4D rock sinogram. The bin size in the x and y axes was \qty{2}{\mm}. The bin size for the polar angle was $\pi/90$ radians. For the azimuthal bins, the bin size depended on which polar bin it is associated with. The number of azimuthal bin at  given polar bin $n_{\theta}$ is governed by equation \ref{eq:total_phi_bins}.

\begin{equation}
    \label{eq:total_phi_bins}
    N_{\phi} = \begin{cases}
            \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;1\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;, & \text{if } n_{\theta} = 0,\\
            \text{round} \left(180 \cdot \text{sin} \left( \frac{n_{\theta}}{45} \cdot \frac{\pi}{2}\right) \right), & \text{otherwise}.
        \end{cases}
\end{equation}

Finally, the 4D sinogram was given to an MLEM algorithm that computed a density map of PET isotopes. The code for the MLEM algorithm used the Astra Toolbox. The output was a 150x150x150 array where a given value represented the estimated density of PET-isotopes for the voxel. Each voxel had a side length of \qty{2}{\mm}. Figure \ref{fig:pet_img} contains an example of a reconstructed PET image of a \qty{10}{\mm} diameter diamond in a rock.

    
\begin{figure}[ht]
    \centering
    \includegraphics[width=0.5\linewidth]{img/PET_image.png}
    \caption{Example of a reconstructed image showing a \qty{10}{mm} diameter diamond (red) in a \qty{100}{mm} diameter kimberlite rock (green) surrounded by noise (blue).}
    \label{fig:pet_img}
\end{figure}

\subsection{Convolutional Neural Network Architecture}
The convolutional neural network (CNN) architecture presented in Table \ref{tab:cnn_arch} was adapted from a 2017 solution to a Kaggle competition on classifying lung tissue and predicting nodule malignancy \cite{kaggle_2017}. The domains of diamond classification and lung tissue classification overlap nicely, given that both have a lot of background noise and the target can be of any size or shape. Furthermore, CNNs are a foundational aspect of deep learning for computer vision and have been shown to excel at image classification tasks.
    
% \begin{figure}[ht]
%     \centering
%     \includegraphics[width=1\linewidth]{img/cnn_architecture_stringified.png}
%     \caption{High-level summary of the CNN classifier architecture.}
%     \label{fig:cnn_arch}
% \end{figure}

\begin{table}[ht]
    \centering
        \label{tab:cnn_arch}
    \begin{tabular}{@{}ccccc@{}}
         \toprule
         Layer & Window & Activation & Output shape & Features\\
         \midrule
         Input &  &  & 21x21x21 & 1\\
         Average & 2x2x2 &  & 11x11x11 & 1\\
         Convolution & 3x3x3 & ReLU & 11x11x11 & 64\\
         Max pool & 2x2x2 &  & 10x10x10 & 64\\
         Convolution & 3x3x3 & ReLU & 10x10x10 & 128\\
         Max pool & 2x2x2 &  & 9x9x9 & 128\\
         Convolution & 3x3x3 & ReLU & 9x9x9 & 256\\
         Convolution & 3x3x3 & ReLU & 9x9x9 & 256\\
         Max pool & 2x2x2 &  & 8x8x8 & 256\\
         Convolution & 3x3x3 & ReLU & 8x8x8 & 512\\
         Convolution & 3x3x3 & ReLU & 8x8x8 & 512\\
         Max pool & 2x2x2 &  & 7x7x7 & 512\\
         Convolution & 2x2x2 & ReLU & 6x6x6 & 64\\
         Convolution & 6x6x6 & Sigmoid & 1x1x1 & 1\\
         \bottomrule
    \end{tabular}
    \caption{High-level summary of the CNN classifier architecture.}
\end{table}

The layers listed in Table \ref{tab:cnn_arch} can be divided into seven stages. The first is the input stage which immediately downsamples the input via averaging so that the training time is faster and more manageable. The next four stages are variations on convolutions followed by a max pool. The first of these captures information about patterns in the data, the subsequent stages build upon the previous ones to create and monitor ever more complex patterns and features. The number of features doubles between the stages to capture the exponential growth of pattern complexity and combinations. The sixth stage is a bottleneck which imposes regularisation constraints on the network and forces it to discard useless features. The final stage is the prediction, which uses a sigmoid activation function because it outputs a value between 0 and 1, allowing us to interpret the output as the CNN's confidence that a sample of rock contains a diamond.

\subsection{Training, Testing, and Data preparation}

The reconstructed PET images were sorted into three groups according to the size of the diamond in the simulation. The diamond diameters in the first group (DG1) ranged from \qty{0}{\mm} to \qty{5}{\mm}, in the second group (DG2) they ranged from \qty{5}{\mm} to \qty{10}{\mm}, and in the third group (DG3) they ranged from \qty{10}{\mm} to \qty{15}{\mm}. The are seven possible ways to combine these groups, each of which will be used to train a CNN. The images were de-noised with Gaussian blurring using a 1 voxel wide kernel. Then the values in the images were rescaled to the interval $[0,1]$ using the global minimum and maximum values of both the training and testing datasets.

Second, each image was subsampled to get 20 diamond samples and 20 barren samples. Each sample is a cube with a side length of 21 voxels. The cubes were sampled by randomising the location of the central voxel. For the 20 diamond cubes the diamond was always fully contained within the cube. While for the 20 barren cubes their central voxel was always in the kimberlite sphere and they were as far away as possible from the diamond. This meant that the 3000 training simulations became 60,000 diamond cubes and 60,000 barren cubes. The test data underwent the same sub sampling treatment as the training data. Giving rise to 6020 diamond cubes and 6020 barren for each of the diamond sizes.

\section{Results and Discussion}

In total, seven CNNs were trained on the seven combinations of the three training dataset groups (DG1, DG2, and DG3). The results are recorded below in Figure \ref{fig:recall_precision_graphs}. These results show that none of the combinations dominate in terms of performance in all four metrics. However, CNNs that were exposed to \qty{0}{\mm}-\qty{5}{\mm} (DG1), generally have higher recall on diamonds less than \qty{7}{\mm} in diameter, but at the same time they let through many false positives. Conversely, the CNNs that weren't exposed to DG1 perform worse on diamonds less than \qty{7}{\mm} in diameter, but almost never let through false positives. This discrepancy is likely down to the cubes containing diamonds less than \qty{2}{\mm} in diameter beginning to fade into the background. This would confuse the CNN training process giving rise to a low threshold for classifying a cube as containing a diamond. More experimentation is needed to verify which subset of diamond training samples allows a CNN to classify smaller diamonds without increasing the false positive rate.

Image quality may also be limiting the performance of the CNNs, the MLEM algorithm does not provide image quality improvements after only one iteration. Improving the image quality may lead to substantial gains in discriminatory power, which could be achieved by improving the MLEM algorithm used.

\begin{figure}[ht]
    \centering
    \includesvg[width=0.8\linewidth]{img/Recall_Precision_FPR_ROCAUC.svg}
    \caption{Plots showing test data performance of seven CNNs: Recall (upper left), ROC area under the curve (upper right), Precision (lower left), and False positive rate (lower right).}
    \label{fig:recall_precision_graphs}
\end{figure}

Out of the seven CNNs, the orange one, trained on DG2 and DG3, is the best performing model. It dominates the pink and brown models because it saw twice as many examples during training. It also maintains a very low false positive rate and a high precision. However, it is not the climax of this study. Importantly it captures only 50\% of diamonds at the \qty{5}{\mm} size, which are a critical component of any mine's revenue source. Moreover, the training data can be improved by increasing the complexity of rock and diamond shapes as well as the number of rocks and the size of cuboid the CNN is trained on. The aim is for a CNN look at \qty{100}{\cm}x\qty{30}{\cm}x\qty{10}{\cm} belt of rocks and decide whether there are any diamonds. Furthermore, the CNN does not tell where a rock is which is crucial information for the ejection system.

\section{Conclusion}

The results provide compelling evidence that a CNN can be trained to tell the difference between diamond containing rocks and barren rocks. The limit of the models trained in this study appears to be related to the identification of diamonds of diameter between \qty{5}{\mm} and \qty{6}{\mm}. Further work needs to be carried out to investigate how to lower this limit without increasing the false positive rate. Steps in this direction include improving the image quality and identifying and removing confounding training examples. These can be achieved by improving the MLEM algorithm. Moreover, the structure and hyperparameters of the CNN were not fully investigated, and this should be addressed in future work, for example by using a large cube size.

More work also needs to be done to increase the realism of the simulations. This can be achieved via increasing the number and shape variations of the rocks and diamonds. Modelling a non-uniform activation field will also help increase the realism and thus reliability of the simulations and any results derived from them.

%Bibliography
%make sure you fill the .bib properly, e.g. edit "biblio.bib" and add all your references in there.
\bibliographystyle{IEEEtran}
\bibliography{references}

\end{document}