Docteur

Olumide Okubadejo

Thèse Dirigée Par

Professeur, UGA Gioacchino Viggiani

M Hugues

Rapporteur M Sébastien

Brisard Ifsttar -Marne-La-Vallée

Rapporteur M Emmanuel Roubin

GOUILLART Mme Emmanuelle

M Pierre-

Yves Coulon

M Gioacchino Viggiani

THÈSE Pour obtenir le grade de

progressive breakage under mechanical loading with image analysis of 3D+t tomography images

Thèse soutenue publiquement le 27

Résumé

La rupture du grain dans les matériaux granulaires a été relativement difficile à calculer et à caractériser dans les images de tomographie. c'est parce qu'il est complexe de formuler un algorithme pour la caractérisation des grains qui bougent et se cassent.

Dans cette thèse, nous mettons en évidence des approches informatiques qui améliorent la compréhension des phénomènes de rupture et de broyage dans les matériaux granulaires. En raison de l'interdépendance de la précision de la segmentation et de la capacité de calcul de la rupture, nous commençons par examiner les techniques d'élimination du bruit dans les matériaux granulaires. Les techniques d'élimination du bruit sont analysées à l'aide d'un ensemble de matériaux auxquels elles s'appliquent. Deuxièmement, nous passons d'une approche morphologique par bassin versant à la segmentation des géomatériaux à une approche hiérarchique qui permet de mieux capturer les informations a priori à partir des sources de données.

Les moyens par lesquels la segmentation spécifique au contexte ou à l'image peut être réalisée sont itérés. Troisièmement, nous présentons un modèle pour capturer les ruptures d'images statiques; sans considération de motion. Enfin, nous présentons des modèles spatio-temporels qui suivent l'évolution de la casse dans les images de matériaux granulaires.

I would equally thank my supervisors at the 3SR laboratory Edward Ando and Cino Viggiani. Thanks for your patience and time in introducing me to the field of Geomechanics. Thank you Eddy for being the perfect beginning to a new field, a new life and a Ph.D. in France. Thank you Cino for pushing me to the end and being more than a friend.

One day, they were born, they were seeds Another day, walked the torrents, high seas On the third, by fire; reborn. They are trees.

Settling in France and understanding France couldn't have been possible without my Friends Dorjan, Ilaria, Jeanne, Mathias, Maddi, Martha, Mohammad and Thanos, otherwise known as the seeds. Thanks for the guidance, friendship, love, drunk discussions and most importantly, memories.

Hitherto

Perfection, we cannot define.

Alberto. "Perfetto"

I would like to thank Alberto Terzolo, the Italian friend I met in France. He carried me on his back when I couldn't walk and spoke on my behalf when I could not talk. I cannot begin to recount the depth of your help. I can say categorically that you are the perfect friend that no one deserves. I would like to acknowledge and thank Dr. Chloe Arson and Pei Wang for the interesting collaboration during the course of the Ph.D.

Girlfriend of the past

May this ring last.

Thank you, Folake.

Family is first? Aghast.

My family; first of all friends

Last line of defense.

Thanks to my Parents Oluwole and Oyefunke. Thanks to my siblings Olubunmi, Oluyemi and Olutomi. You, my family were my last line of defense.

List of Tables

List of Symbols

Problem definition

Grain breakage is an important phenomenon in a large number of industrial problems. These industrial processes include processes that subject granular materials to high stresses. Such industrial processes include the installation of supporting piles into the soil, petroleum extraction, and the installation of wind turbines. The study of grain breakage is also essential when grain breakage is the final goal of a process, e.g., crushing of extracted ores from mining, or particles from pharmaceuticals. Multi-scale models for granular media require quantitative measurements coming from experiments. For a complex phenomenon such as grain breakage, a series of high-resolution tomography images are essential. The ability to quantify critical variables related to grain breakage, e.g., particle size, and shape evolution, is required for modeling. The quantification of such variables is not trivial. We can obtain high-resolution images of grain breakage by X-ray imaging or Neutron imaging. These images mostly consist of a sequence of 3-dimensional sub-images,

showing the temporal change in grains. Although these images are of high resolution, they sometimes contain noise. These noisy components degrade the accuracy of critical variable computations.

In this thesis, we approach the different elements necessary for the computation of variables relating to breakage. We start from the removal of noise in images for which such computation is required, to the identification of granular material and subsequently, temporal tracking of grains even when they break. We outline various methods for different parts.

Challenges associated with the problem

We summarize the steps associated with computing the values of parameters that relate to breakage in these images as;

• obtaining the tomography image

• removal of noise-related artifacts

• identification of individual grains

• Tracking the displacement, rotation, and breakage of grains across the temporal axis

All of these sub-problems need to handled in the same sequence as presented due to the dependence of the next step on the preceding step. The propagation of errors from one step to another is a challenge.

During the identification of grains, it is vital to identify the particles of a broken grain as particles belonging to a single grain. Such identification enables us to calculate the values of grain parameters such as grain volume without loss of precision.

It also allows us to track the evolution of each grain post-breakage without using its unique label ID. However, it is also necessary to identify such grains as broken to characterize its evolution leading to breakage. Thus a grain broken into fragments should implicitly be characterized as broken but its fragments should retain the information of belonging to the original particle.

Overview of thesis

Chapter 2

In Chapter 2, we present the basics of x-ray tomography and methods by which we generate X-ray tomography images of geomaterials. We also present an overview of the problems associated with x-ray tomography image generation. We also present an introduction to breakage of geomaterials; the types of breakage and methods by which they occur. This chapter serves as an introduction to the concepts used throughout this thesis.

Chapter 3

In chapter 3, we present the datasets used during this work. We explain the methods by which some of the naturally occurring materials imaged are formed, the effective content of these materials, and the conditions under which they were imaged. We also present an introduction to artificially generated spheres used extensively during the thesis.

Chapter 4

In chapter 4, we present methods by which we quantify noise in tomography images of granular materials. First, we review the method for the computation of noise presented in [START_REF] Liu | Noise estimation from a single image[END_REF]. Then we present a neural network based substitute that achieves better noise estimation. Our method is based on the encoder structure commonly used in neural networks.

Chapter 5

In chapter 5, we present various edge preserving methods for noise removal present brief observations on the use of each method. In chapter 6, we present a novel method by which granular material properties such as the nature of granular contact can provide better segmentation accuracy. The contact model is used to inform the segmentation and serves a prior. This prior based hierarchical segmentation biases the segmentation towards extracting grains that conform to the prior model.

Chapter 7

In chapter 7, We formulate a method by which we can add texture to artificial grains (code generated) and break these grains. Using our formulation, we study the parameters that affect the computation of breakage outlined in [START_REF] Karatza | Study of temporal and spatial evolution of deformation and breakage of dry granular materials using x-ray computed tomography and the discrete element method[END_REF]. We propose methods for the detection and tracking of broken grains. We compare the segmentation between two successive time-steps and use the disparity in results to account for breakage.

Chapter 8

In chapter 8, We propose a method for the detection and tracking of broken grains.

We compare the segmentation between two successive time-steps and use the disparity in results to account for breakage.

Chapter 9

Here we present our significant conclusions from this work, limitations, and future perspectives.

Chapter 2

Grain breakage: A short literature review 2.1 X-ray basics X-ray CT is also used in the non-destructive volumetric study of rare specimens, e.g., fossils, meteorites, 3D measurement of fluid flow fields, porosity, microporosity, fracture extent and roughness, 3D fabric determination (foliations, shape preferred orientations, network properties). It is suitable for use in many fields and in many contexts. Its extensive use is because it does not destroy the material with which it interacts.

Wilhelm Rontgen discovered X-rays while experimenting with cathode rays in a glass tube [START_REF] Robert | Squire's fundamentals of radiology[END_REF]. The discovery of x-rays motivated research into its use and its behavior [START_REF] Eric | The birth and early infancy of X-rays[END_REF], [START_REF] Fr Arutyunian | The Compton effect on relativistic electrons and the possibility of obtaining high energy beams[END_REF]. One such experiment was done by Arthur Compton [START_REF] Fr Arutyunian | The Compton effect on relativistic electrons and the possibility of obtaining high energy beams[END_REF], in which he directed x-rays to a metal surface. He discovered that the flux of x-rays decreased upon scattering. He observed that metal surface ejected electrons when x-rays hit the surface. He also observed an increase in the wavelength of the beam, inferring a reduction in x-ray flux. The absorption of sufficient x-ray led to the release of free electrons. As x-rays interact with materials, their change in intensity correlates to Figure 2.1: A typical x-ray system for imaging materials the atomic number of the material. [START_REF] Fr Arutyunian | The Compton effect on relativistic electrons and the possibility of obtaining high energy beams[END_REF], [START_REF] Ando | Experimental investigation of microstructural changes in deforming granular media using x-ray tomography[END_REF].

In order to reconstruct a 3D volume of x-ray attenuation, we rotate the sample or the source-detector system. The source emits x-rays, and the detector captures the x-rays after material interaction. This is shown in figure 2.1. The rotation of the sample or the source-detector takes place at angle intervals. At every rotation angle, the detector captures the attenuated beam resulting in image projections. These projections capture material property variance at these angular intervals. A flat source and detector, emitting x-rays and recapturing them respectively, can do so in two dimensions. We obtain a 3-dimensional image representation by "reconstructing" the projections of the attenuation field into a single map of attenuation.

To explain the process of reconstruction, we introduce the projection-slice theorem in 1.

Theorem 1 (Projection-Slice Theorem ). In N dimensions, the projection-slice theorem states that the Fourier transform of the projection of an N -dimensional function f (r) onto an m-dimensional linear submanifold is equal to an m-dimensional slice of the N -dimensional Fourier transform of that function consisting of an mdimensional linear submanifold through the origin in the Fourier space which is parallel to the projection submanifold.

The projection-slice theorem in 1 postulates that a perfect reconstruction requires an infinite amount of projections. In reality, an infinite number of projections will take an infinite amount of time to acquire. Therefore, perfect reconstruction is impossible. Apart from reconstruction errors due to the impossibility of perfect reconstruction, noise sources can affect the quality of reconstruction. These sources range from human errors to machine and instrumentation errors. Together these form the basis for noise in reconstructed images.

The mathematics of reconstruction

Definition 1 (Beer-Lambert Law). For a monochromatic, non-refractive zero width x-ray beam that traverses a homogenous material along a distance z from the origin, the intensity J(z) is given by

J(z) = J 0 (z) exp -τ (z)dz (2.1)
where -τ (z)dz is the coefficient of attenuation being integrated through a variation in z. This highlights x-rays and x-ray attenuation energy as being a function of position.

We describe the total attenuation p at position x and y as

Γ(r) = ln J J 0 -τ (z)dz (2.2)
Equation 2.2 states that the ration of the input ray intensity over the output ray intensity after a logarithm operation represents the line integral of the attenuation coefficients along the x-ray path. The problem of reconstruction can thus be stated as: given the line integrals of an object, how can we calculate its attenuation distribution?

Based on the polar geometry, r is equal to x cos θ + y sin θ, when the projection is at the angle θ

This thus simplifies to equation 2.3 Which describes the radon transformation of an unknown density function f (x, y).

Γ θ (r) = ∞ -∞ ∞ -∞ f (x, y)δ(x cos θ + y sin θ -r)dxdy (2.3)
As with Fourier transformations, the inverse transform of the Radon equation should theoretically yield the original image.

f (x, y) = 1 2π π 0 g θ (r)(x cos θ + y sin θ)dθ (2.4)
Due to the limited number of projections, several algorithms have been designed, based on the Radon transform to perform the reconstruction projections into tomography images. We examine some of these reconstruction algorithms.

Fourier domain reconstruction

Consider a sample f (x, y) with x-ray sampling at equally spaced angular intervals.

Fourier domain reconstruction transforms each projection using the discrete Fourier transform. Its output is the corresponding frequency spectra on a polar raster. Gaps (sparsity) created in the frequency domain on the raster are fixed, using various interpolation techniques. See [START_REF] Peters | Principles of Computerised Tomography[END_REF], [START_REF] Magnusson | Linogram and other direct Fourier methods for tomographic reconstruction[END_REF].

The inverse discrete Fourier transform of the raster yields the image approximation. Sparsity in the frequency domain is a cause for high noise levels in reconstruction by this method. Algorithmic improvements have been made to change sparsity and thus reduce image noise levels.

A concentric square raster can be generated by changing the angle between projections as follows

θ = R 0 max{| cos θ|, | sin θ|}
The use of the concentric square raster has the advantage of generating a uniformly spaced rectangular grid in the frequency domain and inherently reducing interpolation errors and in effect, reconstruction noise.

Filtered Back-projection reconstruction

The filtered back-projection algorithm is one of the most commonly used algorithms for tomographic image reconstruction. There are two steps in the filtered backprojection algorithm. The first step is filtering the projections by a kernel whose frequency response is |ω|. This can be visualized as a simple weighting of each projection in the frequency domain. The next step is the back-projection step, which consists of projecting back each view through the line corresponding to the direction in which the projection data was collected. The derivation of the filtered back-projection algorithm can be seen in [START_REF] Hsieh | Computed tomography: principles, design, artifacts, and recent advances[END_REF]. The choice of |ω| corresponds to the ramp filter. Practically, the ramp filter amplifies statistical noise in projection data.

Iterative reconstruction

Iterative reconstruction techniques are a set of techniques in which the image reconstruction is obtained by solving a set of linear equations. These methods have gained traction over the years because of the ease with which apriori information can be introduced. See [START_REF] Jeffrey | Penalized weighted least-squares image reconstruction for positron emission tomography[END_REF], [START_REF] Gabor | Fundamentals of computerized tomography: image reconstruction from projections[END_REF], [START_REF] Wang | High performance model based image reconstruction[END_REF], [START_REF] Adler | Learned primal-dual reconstruction[END_REF], [START_REF] Hammernik | Learning a variational network for reconstruction of accelerated MRI data[END_REF] for detailed insight into these set of algorithms.

Noise and artifacts in tomography images

We evaluate reconstruction algorithms based on how close the reconstructed image is to reality. We define reality as the true attenuation coefficients [START_REF] Hsieh | Computed tomography: principles, design, artifacts, and recent advances[END_REF]. Severe deviation from reality might occlude image features and reduce the capacity for both observational inference and quantitative inference. The deviation from reality in the resulting images manifests as inconsistencies in the reconstructed image. These inconsistencies can be a result of faulty data collection process, instrumentation, or system design. These usually result in artifacts, blur, or noise (random sources).

The combination of artifacts, blur, and noise degrades the quality of a CT image.

Noise sources

Noise impairs the visibility of low contrast objects [START_REF] Sprawls | AAPM tutorial. CT image detail and noise[END_REF]. It is characterized as the addition of unwanted information or data to pixels (picture elements). If the aggregate information contained in a group of picture elements is not significantly distinct from their environment (high contrast), noise erodes all or a part of such aggregate information [START_REF] Sprawls | AAPM tutorial. CT image detail and noise[END_REF], [START_REF] Behrendt | Image fusion in dual energy computed tomography: effect on contrast enhancement, signal-to-noise ratio and image quality in computed tomography angiography[END_REF]. The addition of noise is mostly due to the erroneous measurement of pixel attenuation characteristics. Such faulty measurements can be due to equipment faults, human error, or a variety of other sources. Noise erodes image structures and often limits the ability to visualize image structures.

Due to the variety of inputs by which noise can affect an image, it is challenging to mitigate during scanning and reconstruction [START_REF] Stephen O Rice | Mathematical analysis of random noise[END_REF]. It is thus modeled in literature as an additive Gaussian with a mean of 0 (additive white Gaussian) [START_REF] Rosenblatt | A central limit theorem and a strong mixing condition[END_REF], [START_REF] Cc Heyde | Central limit theorem[END_REF].

Artifacts

Different kinds of artifacts can be introduced into the resulting image during scanning or reconstruction. These include ring artifacts, beam hardening, aliasing arti-facts such as partial volume effect.

Ring artifacts [START_REF] Kyriakou | Ring artifact correction for high-resolution micro CT[END_REF], [START_REF] Münch | Stripe and ring artifact removal with combined wavelet-Fourier filtering[END_REF], [START_REF] Prell | Comparison of ring artifact correction methods for flat-detector CT[END_REF] are introduced when there is faulty calibration in the scanner detector. This faulty calibration gives a consistently erroneous reading at each angular position, resulting in a circular artifact.

Beam hardening on the other is caused by selective attenuation in x-rays with polychromatic energies [START_REF] Rodney | Beam hardening in x-ray reconstructive tomography[END_REF]. They can also be manifest as streaking artifact [START_REF] Rodney | Beam hardening in x-ray reconstructive tomography[END_REF], [START_REF] Richard | A pattern recognition method for the removal of streaking artifact in computed tomography[END_REF] or the cupping artifact [START_REF] Rodney | Beam hardening in x-ray reconstructive tomography[END_REF], [START_REF] Mc Altunbas | A post-reconstruction method to correct cupping artifacts in cone beam breast computed tomography[END_REF], [START_REF] Yang | Cupping artifact correction and automated classification for high-resolution dedicated breast CT images[END_REF]. In very different cross sections, dark bands or streaks can appear between two dense objects in an image. They occur because the portion of the beam that passes through one of the objects at certain positions is hardened less than when it passes through both objects at other positions. The resulting artifact is called the streaking artifact.

X rays passing through the middle portion of a uniform cylindrical phantom are hardened more than those passing through the edges because they are passing through more material. As the beam becomes harder, the rate at which it is attenuated decreases, so the beam is more intense when it reaches the detectors.

Therefore, the resultant attenuation profile differs from the ideal profile that would be obtained without beam hardening. The resulting artifact is the cupping artifact.

Beam artifacts can be reduced using adequate beam filtration [START_REF] Hsieh | Computed tomography: principles, design, artifacts, and recent advances[END_REF].

Partial-Volume Effects Partial volume effect results as errors due to discretization. Because each pixel in a computed tomography image represents the attenuation properties of a specific material volume, if that volume is comprised of several substances, then the resulting computed tomography value represents some average of their properties. This is termed the partial-volume effect. Furthermore, because of the inherent resolution limitations of X-ray computed tomography, all material boundaries are blurred to some extent, and thus the material in any one voxel can affect computed tomography values of surrounding voxels [START_REF] Andrew L Alexander | Analysis of partial volume effects in diffusion-tensor MRI[END_REF]. The effect of partial volume is thus; material edges are not represented as hard edges in the images, but as soft ramp edges, which is a linear combination of material and void.

Blur

Blur can be described as the smoothing of pixel information. It can also be described as the inadequacy of information in picture elements. Blur can be introduced during scanning when a large spot size and a small detector size is used [START_REF] Sprawls | AAPM tutorial. CT image detail and noise[END_REF], [START_REF] Leitgeb | Performance of fourier domain vs. time domain optical coherence tomography[END_REF]. The reconstructed image is blurred if the voxel size is smaller than the characteristic sharpness of the beam [START_REF] Boas | CT artifacts: causes and reduction techniques[END_REF]. It is vital to set reconstruction parameters in a way that does not result in a blurred image. Many techniques can be used to reduce blur in tomography images such as the point spread function [START_REF] Myasnikova | Estimation of errors in gene expression data introduced by diffractive blurring of confocal images[END_REF], iterative Richardson-Lucy algorithm, and the Laplacian sharpening filter.

Grain breakage

Breakage can be observed in many engineering and geological processes and structures. These includes earth dams [START_REF] Kenneth | Compressibility and crushing of granular soil in anisotropic triaxial compression[END_REF], pile driving [START_REF] Yasufuku | Pile end-bearing capacity in crushable sands[END_REF], railway ballasts [START_REF] Lim | Discrete element modelling of railway ballast[END_REF], rapid penetration [START_REF] Omidvar | Response of granular media to rapid penetration[END_REF] and shearing of granular faults gouge [START_REF] Abe | Grain fracture in 3D numerical simulations of granular shear[END_REF]. Many factors are responsible for the breakage of granular materials. Among these factors are loading mode, microstructural heterogeneities, stresses, and temperature. The high number of influences causing grains to break makes breakage a complex phenomenon to study. When we subject a granular assembly to high confining pressure, the mean stress within the assembly rises. When the mean stress gets high, breakage occurs.

In scientific research, it has been observed that under high effective stresses, breakage of particles or particle failure might occur. Particle breakage can change the physical and mechanical behavior of a granular assembly [START_REF] Jaeger | Failure of rocks under tensile conditions[END_REF], [START_REF] Bobby | Crushing of soil particles[END_REF], [START_REF] Poul V Lade | Significance of particle crushing in granular materials[END_REF], [START_REF] Mcdowell | The fractal crushing of granular materials[END_REF], [START_REF] Mr Coop | Particle breakage during shearing of a carbonate sand[END_REF], [START_REF] Yang | Sand grain crushing and interface shearing during displacement pile installation in sand[END_REF]. Particle breakage in an assembly affects the grain size distribution [START_REF] Mcdowell | The fractal crushing of granular materials[END_REF]. The grain size distribution (GSD) describes the relative sizes of grains present in an assembly. Particle breakage also affects the shear strength [START_REF] Daouadji | An elastoplastic model for granular materials taking into account grain breakage[END_REF], [START_REF] Lobo-Guerrero | Visualization of crushing evolution in granular materials under compression using DEM[END_REF]. Shear strength is the resistance of a material to a load that can cause sliding failure. The breaking of particles can also affect the solid fraction [START_REF] Md Bolton | Micro-and macro-mechanical behaviour of DEM crushable materials[END_REF], which is the fraction of the assembly volume that is solid. It can also affect the yield surface [START_REF] Miura | Stress-strain characteristics of sand in a particle-crushing region[END_REF].

Parameters influencing breakage

Shape of material

The effect of the shape of a material relative to soil behavior was summarised in [START_REF] Santamarina | Soil behaviour: The role of particle shape[END_REF] using three principles. These are sphericity, angularity, and roughness. Sphericity influences fabric anisotropy by increasing stiffness and residual friction angle. Increased angularity, on the other hand, can cause a decrease in small-strain stiffness and an increase in high-strain strength. In a bid to quantify the effect of particle shape on particle crushing behavior, Afshar et al. in [START_REF] Afshar | Impact of particle shape on breakage of recycled construction and demolition aggregates[END_REF] built a DEM model to simulate different shaped recycled construction and demolition materials. The effect of shape on particle crushing was evaluated. Granular materials with shapes closer to spheres were evaluated to be more resistant to crushing.

Behaviour of material

Different mechanical properties can influence the deformation behavior of materials. These properties include; tensile strength [START_REF] Roberts | The compressibility of sand[END_REF] and the shear strength [START_REF] Bono | DEM of triaxial tests on crushable cemented sand[END_REF] of the material. In [START_REF] Billam | Some aspects of the behaviour of granular materials at high pressures[END_REF], Billam et al. observed that the crushing force is dependent on the nature of the material. Lo et al. also observed in [START_REF] Kwan Yee Lo | Response of particulate materials at high pressures[END_REF] the dependence of particle crushing on the particle mineralogy. However, in Yamamuro et al. [START_REF] Yamamuro | Drained sand behavior in axisymmetric tests at high pressures[END_REF], it is observed that the fundamental influencing properties of granular material relative to crushing is the inherent strength and the effective stress rate.

Grain size distribution

The effect of grain size distribution has previously been analyzed from two perspectives. Marsal [START_REF]Discussion of shear strength[END_REF] indicated that a large particle size could lead to large contact normals, which in turn amplifies the possibility of breakage. The probability of large particles undergoing breakage is increased due to an increase in the possibility of inclusions in large-sized granular materials.

Hardin [START_REF] Bobby | Crushing of soil particles[END_REF] analyzed data of several granular materials with different sized particles under loading. He used b p to describe the potential for breakage of a particle whose size is D.

b p = log 10 D D min f or D ≥ D min (2.5) b p = 0 f or D ≥ D m in (2.6)
D min refers to the minimum size of particle which can crush for any given particle. This value is usually evaluated through laboratory tests.

Where b p represents the potential for breakage for a given size fraction, the potential for breakage for the material is as given using equation 2.7.

B p = 1 0 b p df (2.7)
The amount of crushing is represented as

B t = 1 0 b p0 -b pl df (2.8)
Figure 2.2: Relationship between total breakage and breakage potential [START_REF] Bobby | Crushing of soil particles[END_REF] b p0 equals the initial value of b p and b pl is the value of b p after loading. Measurements in [START_REF] Bobby | Crushing of soil particles[END_REF] show that the value B t is directly proportional to B p .

Thus relative breakage is characterized as

B r = B t B p (2.9)

Relative humidity

Much research has focused primarily on breakage as a function of the mechanical and behavioral properties of granular materials. However, breakage can be influenced by external factors such as relative humidity [START_REF] Luciano | Suction effects on rockfill compressibility[END_REF], [START_REF] Chester | Subcritical creep compaction of quartz sand at diagenetic conditions: Effects of water and grain size[END_REF], [START_REF] Rh Brzesowsky | Compaction creep of sands due to time-dependent grain failure: Effects of chemical environment, applied stress, and grain size[END_REF]. In these works, [START_REF] Luciano | Suction effects on rockfill compressibility[END_REF], [START_REF] Chester | Subcritical creep compaction of quartz sand at diagenetic conditions: Effects of water and grain size[END_REF], [START_REF] Rh Brzesowsky | Compaction creep of sands due to time-dependent grain failure: Effects of chemical environment, applied stress, and grain size[END_REF], it is shown that exposure to water can lead to breakage by lowering the yield strength. Chester et al. in [START_REF] Chester | Subcritical creep compaction of quartz sand at diagenetic conditions: Effects of water and grain size[END_REF] show that the injection of water into dry samples can lead to volumetric collapse. Thus, in the presence of water, the required stress for the critical failure of a specimen through fracture is reduced. This effect of water vapor on fracture has also been observed in [START_REF] Sm Wiederhorn | Influence of water vapor on crack propagation in soda-lime glass[END_REF], [START_REF] Kean | Subcritical crack propagation in rocks: theory, experimental results and applications[END_REF], [START_REF] Nara | Influence of relative humidity on fracture toughness of rock: implications for subcritical crack growth[END_REF].

Influence of breakage on material properties

Breakage can lead to changes in an assembly of grains. Three parameters that can be affected significantly are the shape of the particles and the grain size distribution.

Particle shape

Angular particles cause stress concentrations and thus crush under the influence of a load. Thus as breakage progress, the mean shape of particles become less angular.

According to Ueng et al. in [START_REF] Ueng | Energy aspects of particle breakage in drained shear of sands[END_REF], less angular particles could lead to a reduction of the friction angle and thus induce a decrease in the shearing strength.

Grain size distribution

The breaking of particles effectively leads to a reduction in the size of the grain and thus influence the grain size distribution of the grain assembly. Einav in [START_REF] Einav | Breakage mechanics-part I: theory[END_REF], [START_REF] Einav | Breakage mechanics-Part II: Modelling granular materials[END_REF], describes breakage by a change in the surface of the grain size distribution curve. Changes in grain sizes result in changes to the grain size distribution curve.

According to Einav [START_REF] Einav | Breakage mechanics-part I: theory[END_REF], [START_REF] Einav | Breakage mechanics-Part II: Modelling granular materials[END_REF], the surface of the curve before loading is higher than the surface of the curve after loading which is in turn higher than the surface of a size distribution state called the ultimate distribution. The ultimate distribution is a grain state where further breakage can no longer occur.

Breakage types

Breakage is a way to dissipate an energy build up within a material. It is an irreversible change of shape or size of a particle. When cracks appear within materials, the localized stresses within a grain are modified. With an increase in the build-up of stresses, cracks propagate, resulting in breakage [START_REF] Roylance | Engineering viscoelasticity[END_REF].

Daouadji et al. in [START_REF] Daouadji | An enhanced constitutive model for crushable granular materials[END_REF] identified three breakage modes in granular materials.

These are fracture, attrition, and abrasion, as can be seen in figure 2.4. The difference between these breakage modes is characterized by a difference in the shape and size of the particles. During breakage, some parametric changes such as changes in shape and grain size distribution can be observed, according to [START_REF] Daouadji | An enhanced constitutive model for crushable granular materials[END_REF].

Single particle breakage tests

Single particle breakage tests study the effects and characteristics of single grains undergoing breakage. Studying one particle at a time allows fine-grained control and easy tracking of broken particles. Although some results obtained can be extrapolated to understand breakage in grain assemblies, it is not representative of contact topologies and coordination numbers found in the grain assemblies. In [START_REF] Gr Mcdowellć | On the micromechanics of crushable aggregates[END_REF],

McDowell observed that particle breakage was the principal source of plastic volumetric compression on the normal compression line. He relates the micromechanics of grain fracture to the macroscopic deformation of crushable aggregates by performing one-dimensional compression tests. In one-dimensional compression tests, a vertical load is applied with restraints applied along the lateral axis. This vertical load results in vertical deformation without lateral deformation. McDowell [START_REF] Gr Mcdowellć | On the micromechanics of crushable aggregates[END_REF] observed high yield stress to indicate the start of particle fracture. With the same test framework, [START_REF] Yoshimoto | Evaluation of shear strength and mechanical properties of granulated coal ash based on single particle strength[END_REF] observed that when crushing strength increases, the shear stiffness increased, and the volumetric strain decreased.

Other types of single particle crushing tests have also been used to evaluate particle crushing properties. Salami et al. in [START_REF] Salami | An experimental study on the influence of the coordination number on grain crushing[END_REF] performed several multi-point crushing tests on cylindrical samples to investigate the effect of coordination number on crushing properties. Salami et al. [START_REF] Salami | An experimental study on the influence of the coordination number on grain crushing[END_REF] demonstrated that different contact orientations can influence the fragmentation of a single particle and that the critical force was highly dependent on the contact area. Nakata et al. in [START_REF] Nakata | One-dimensional compression behaviour of uniformly graded sand related to single particle crushing strength[END_REF] evaluated the relationship between crushing characteristics and one-dimensional compression test.

Nakata et al. observed that particle stress grows as void ratio increases, and the coordination number decreases.

Zhao et al. in [START_REF] Zhao | An investigation of single sand particle fracture using X-ray micro-tomography[END_REF] studied particle crushing behavior using a computed tomography scan of four Leighton Buzzard sand and four highly decomposed granite particles. They demonstrated that complicated morphology and microstructure pattern often results in intricate fracture patterns. They observed that simpler morphology often yields in breakage that is parallel to the loading direction. They also observed that particle fragments often converge to the same sphericity and convexity levels.

Todisco et al. [START_REF] Mc Todisco | Multiple contact compression tests on sand particles[END_REF] observed from single particle compression tests that altering the coordination number affected breakage. They observed that a coordination number greater than four generally resulted in breakage by splitting. Abrasion and attrition showed in figure 2.4 were observed to occur with coordination number equal to four.

Multi particle breakage tests

A change in the grain size distribution of an assembly is usually indicative of breakage [START_REF] Bobby | Crushing of soil particles[END_REF]. It not so much only the breakage of single particles, but the effect of particle breakage on a grain assembly. However, the bridge between the breakage of a single particle and full-scale measurements in an assembly is not trivial [START_REF] Guida | Breakage mechanisms of highly porous particles in 1D compression revealed by X-ray tomography[END_REF].

Karatza et al. in [START_REF] Karatza | Study of temporal and spatial evolution of deformation and breakage of dry granular materials using x-ray computed tomography and the discrete element method[END_REF] studied the evolution of grain size distribution and other breakage parameters in Zeolite using X-ray tomography imaging. They imaged the deformation and reconstructed the projections into a 4-dimensional stack of 3-dimensional images. To track grains through time, identified grain labels in a 3D image are correlated to the next image using digital image/volume correlation described in [START_REF] Stephen | Discrete and continuum analysis of localised deformation in sand using X-ray [mu] CT and volumetric digital image correlation[END_REF], [START_REF] Andò | Experimental micromechanics: grain-scale observation of sand deformation[END_REF], [START_REF] Tudisco | TomoWarp2: a local digital volume correlation code[END_REF]. The degree of match of a grain label to the next timestep (deformed image) is quantified using the normalized correlation coefficient [START_REF] John | Fast template matching[END_REF]. A normalized correlation coefficient of the maximum value of 1.0 is reflective of a perfectly correlated label. To identify broken grains, Karatza et al. in [START_REF] Karatza | Study of temporal and spatial evolution of deformation and breakage of dry granular materials using x-ray computed tomography and the discrete element method[END_REF] imposed a threshold of 0.94 on the normalized correlation coefficient. They characterized grains with a normalized correlation coefficient value lower value than 0.94 as broken.

An analysis of the extracted grains confirms the relationship between breakage type and coordination number initially observed in [START_REF] Mc Todisco | Multiple contact compression tests on sand particles[END_REF].

Guida et al in [START_REF] Guida | Breakage mechanisms of highly porous particles in 1D compression revealed by X-ray tomography[END_REF], adopted a similar set-up as [START_REF] Karatza | Study of temporal and spatial evolution of deformation and breakage of dry granular materials using x-ray computed tomography and the discrete element method[END_REF] to study breakage mechanisms in glsleca. As opposed to sand grains, LECA particles break at relatively low-stress levels [START_REF] Guida | Breakage mechanisms of highly porous particles in 1D compression revealed by X-ray tomography[END_REF]. They also have high porosity.

Hurley et al. in [START_REF] Hurley | In situ grain fracture mechanics during uniaxial compaction of granular solids[END_REF] performed an extensive study on the mechanics of breakage by combining 3D x-ray tomography and 3D x-ray diffraction. X-ray tomography gave insight into the granular and material properties, while x-ray diffraction study, allowed access to localized force statistics. Both systems provided combined access to a variety of local and global measurements on the sample, including grain kinematics, contacts, average intra-granular stresses, and inter-particle forces.

Conclusion

In this chapter, we introduce the basics of x-ray computed tomography image generation. We show how samples are imaged using an x-ray scanner. The sample or the source-detector on the scanner is rotated at constant angles to generate projections for each angle. We highlight how these projections can be combined into a 3-dimensional image using reconstruction techniques. Further, we highlight the problems with reconstruction. Next, this chapter reviews previous research relating to breakage in granular materials. We identify the primary breakage modes in granular materials under loading. These are splitting, abrasion, and attrition, as shown in figure 2.4. Next, we highlight research studying the breakage of single grains under load. We indicate that although single particle breakage gives fine-grained control over the process (tracking of broken particles), results obtained are difficult to extrapolate to grain assemblies. In single particle crushing, breakage is shown to be dependent on mechanical properties such as yield stress, shear stiffness, and volumetric strain. Multi-particle tests highlight the effect of localized breakage on the global behavior of an assembly of granular materials.We highlight research that shows that breakage affects global shape behavior and size distribution.

Chapter 3 Datasets

Introduction

In this chapter, we present the datasets used in this work. The datasets presented are granular materials scanned and reconstructed at the laboratoire 3SR, Grenoble. The

x-ray scanner in Laboratoire 3SR (Grenoble) is a built-to-specification laboratory

x-ray tomograph supplied by RX-Solutions (Annecy). The cabin in which the x-ray scanner is housed is lead-lined. The x-ray source in this tomograph emits a conebeam in the direction of the detector. By implication, as one translates an observed specimen from the detector towards the source, the specimen is enlarged on the detector due to geometric magnification. [START_REF] Ando | Experimental investigation of microstructural changes in deforming granular media using x-ray tomography[END_REF].

Some specifics on the construction and mounting of this scanner makes it easy to carry out both triaxial and ooedometric tests. These include:

• Translation trolley

The trolley, on which the rotation stage can translate in the direction of the beam, allows for easy control of the zoom level for the specimen being scanned.

It is made of solid tell component and can thus carry high loads.

• Room for experiments The x-ray cabin is of considerably larger width and height than in an ordinary x-ray scanner. This considerably large size means that large amounts of experimental equipment can be mounted inside the cabin.

• Connection to outside There is the easy passage of tubes in and out of the cabin, resulting in less exposure to x-ray radiation.

• Translatable source and detector The source and detector can be moved synchronously, which allows experimental installations of different heights to be scanned. All of the datasets which are subsequently presented were scanned using this

x-ray scanner at the Laboratoire 3SR.

Datasets

Sapphire Spheres

Sapphire spheres is a synthetically produced material. This dataset's underlying sample had almost perfect spheres of 5 different sizes. The sizes of the grains range from 300 µm to 600 µm in increments of 100 µm. The manufacturer (Sandoz Fils S.A.

in Switzerland) grade of the spheres is between grade 5 and 10, meaning sphericity tolerances between 0.125 µm to 0.250 µm and diameter tolerances of ±1 µm to 2.5 µm.

In the corresponding tomography generated image, the pixel size is 7.50 µm/px.

A slice of sapphire spheres is shown in figure 3.2. (Cu = 1.07), with minimum and maximum particle sizes of 1.09mm and 1.50mm, respectively. Each zeolite particle has a density of 2.18gcm3 and an estimated crush strength of 15N (based on the manufacturer's specifications). The zeolite granules are produced from zeolite powder (2µm diameter) with a mineral binder like clay.

Henceforth, zeolite will refer to the granules and not the powder. Zeolites are water-insoluble, isotropic and highly hygroscopic and hence very sensitive to moisture and temperature changes [START_REF] Karatza | Study of temporal and spatial evolution of deformation and breakage of dry granular materials using x-ray computed tomography and the discrete element method[END_REF]. A slice and 3-dimensional specimen of zeolite are shown in figure 3.3.

Leighton Buzzard

Highly Decomposed Granite

Highly decomposed granite is derived from granite via its weathering to the point that the parent material readily fractures into smaller pieces of weaker rock. It has practical uses that include its incorporation into paving and driveway materials, residential gardening materials in arid environments, as well as various types of walkways and heavy-use paths in parks. Highly decomposed granite is a nonhomogenous material, and it is elongated and angular. Images were obtained during an Oedometric compression test performed at the Laboratoire 3SR.

Caicos Ooid

Caicos ooid comes from Ambergis Shoal in the Caicos platform in the British West Indies. Caicos ooids grow in marine environments and are formed through material accretion. This mode of formation implies that they start from a seed point, which can be small quartz or calcite and continuously grow over time [START_REF] Ando | Experimental investigation of microstructural changes in deforming granular media using x-ray tomography[END_REF]. The grains tend to grow by physical attachment as they roll or through precipitation over its surface. This usually generates well-rounded grains. Images for Caicos Ooids were obtained from a triaxial compression test. Evident is a significant amount of blur.

Also, visual noise such as the ring artifact is present in the image. A slice of this image is shown in figure 3.6 The silica particles making up these clumps have been crushed in-situ from larger rocks. Since particles have been generated by crushing, they are generally angular.

Hostun sand grains are also visibly angular in the x-ray tomography images shown in Figure 3.7 [START_REF] Ando | Experimental investigation of microstructural changes in deforming granular media using x-ray tomography[END_REF].

Kalisphera

Kalisphera is an artificial sphere generation tool used to produce 3D raster images of spheres, including their partial volume effect. The ability to generate these spheres,

gives the capacity to quantitatively evaluate the metrological performance of our image-based measurement techniques with reference to a ground truth. [START_REF] Tengattini | Kalisphera: an analytical tool to reproduce the partial volume effect of spheres imaged in 3D[END_REF]. In this work, kalisphera is mostly to gain quantitative insight into the performance of developed algorithms due to the access it gives to the ground truth. An example of the result of this data generation tool is shown in figure 3.8.

Conclusion

Here we gave insight into the datasets used in this work. We briefly present details about naturally occurring grains, their shape, the conditions under which they were Quantifying Noise

Introduction

Tomography imaging has been used extensively for the characterization and understanding of geomaterials. Images are first scanned using an x-ray scanner, and the resulting radiographs are used to reconstruct a 3D image of the imaged sample.

However, it is difficult to obtain a perfect reconstruction void of noise due to the numerous noise sources present due to imaging and reconstruction. Due to the variety of noise inputs, noise in tomography images of geomaterials is usually approximated as a zero mean Gaussian. Although there has been research into the effect of noise on reconstruction accuracy, there has been limited research into quantifying the noise in these reconstructed images. We propose that the quantification of noise in these images might aid in setting the hyper-parameters of noise removal algorithms. Liu et al. in [START_REF] Liu | Noise estimation from a single image[END_REF] proposed a model-based approach to quantifying noise in images. It quantifies the noise as the standard deviation value of the most occurring mean-standard deviation pair computed within pixel neighborhoods. We propose an alternative noise quantification method based on the encoder model of neural networks. The proposed method would not require a quantification on the entire image but can make noise inference from part of the image as small as a single slice.

Background

Noise can be formulated as the Minkowski distance between the true reconstruction and the generated reconstruction. However, in reality, this is near impossible to compute. This is due to the ill-posed formulation where a perfect reconstruction is necessary to compute noise. To quantify noise (approximate) noise, without reference to the perfect reconstruction, [START_REF] Liu | Noise estimation from a single image[END_REF] proposed a solution based on a statistical analysis of the pixel values. In [START_REF] Liu | Noise estimation from a single image[END_REF], the method proposed computes the mean and standard deviation within the neighborhood of each pixel and then formulates noise as the standard deviation in the most occurring mean-standard deviation pair.

Thus the histogram counts the observations that fall into disjoint mean-standard deviation categories. Each category is defined using a bin. This method gives the estimation of the noise in the image. In this formulation, there is an implicit assumption of noise to be a zero-mean additive Gaussian and thus is the reason it can be characterized using a scalar standard deviation parameter.

Consider an noisy reconstructed image I. Its equivalent hypothetical noiseless image is denoted as I ∞ and it noise free estimation as I e . if

I e ≈ I ∞
then equation can be rewritten as

σ 2 = (|I| p -|I e | p ) 1/p (4.1)
where I e = f (I) A potential solution is to partition the image into constant volume regions with no notion of homogeneity or heterogeneity (neighborhoods). The mean and standard deviation(which we classify as noise) is thus computed on these neighborhoods and by an aggregation function projected to the entire image. This resembles an intuitive explanation for the formulation developed in [START_REF] Liu | Noise estimation from a single image[END_REF].

A
We proceed to validate the effectiveness of this formulation for the computation of noise using Kalisphera. In [START_REF] Tengattini | Kalisphera: an analytical tool to reproduce the partial volume effect of spheres imaged in 3D[END_REF], the tool Kalisphera is provided to aid in the creation of artificial 3D spheres. The effect due to Partial Volume is simulated on the created spheres to mimic this effect as obtained in Natural Images nearly. Since additive white Gaussian noise can be added numerically to an otherwise perfect artificial assembly of grains, we use the defined formulation to compute the noise and compare it with the ground truth. Due to the number of inclusions in Zeolite, as the number of bins is increased, the value of noise increases. This is a result of these inclusions creating high frequencies in the histogram, reducing the averaging effect of bins by increasing the number of bins, highlights these high frequencies.

Images, visually hypothesized to contain blur are deemed less noisy. In images like hostun sand having structured noise with predictable repeatability like the ring artifact, there is, little effect on the noise values.

In the natural images, it is difficult to evaluate the accuracy of the noise values computed since we do not know the true noiseless image.

A deep neural network approach

Let I be the tomography image for which we want to estimate noise. Assuming a linear mixing model, mathematically the observed image can be represented as Noise in tomography due to the multiplicity of inputs is usually represented as Gaussian of zero mean.

I = I ∞ + n (4.2)
f (n) = 1 √ 2πσ 2 exp - n 2σ 2 (4.3)
As such, the distribution n is characterizable by its variance. The task of estimating the noise in I is thus one of finding the variance.

It is easy to find the variance σ 2 if the noise n is separated from the mixed model I. However, this is mostly not the case since I ∞ is unknown. Since I is the only known value, the task of estimating the noise n, parameterizable by its estimated variance σ2 is defined as

σ2 = ρ(I) (4.4)
The function ρ(I) in reality can be difficult to define. As such we define ρ(I) as a neural network. The task is thus to train a neural network, parameterized by θ such that the loss L. We define the loss as the difference between the predicted noise noise and the true noise as shown in equation 4.5

L = || σ2 -σ 2 || 1 (4.5)
Where σ2 is the estimated noise and σ 2 is the true noise value. Once minimized on the training set (which are images for which σ 2 is known), it can then be used to estimate the noise on images for which the input noise is not known.

Architecture details

The architecture used takes inspiration from the encoder model described in [START_REF] Bengio | Learning deep architectures for AI[END_REF].

We show how the encoder model can be used to quantify noise in tomography images of geomaterials. Our neural network model consist of two major repeating layer; the convolutional layer and the downsampling layer.

Convolution Layer

A Convolutional Neural Network (ConvNet/CNN) is a Deep Learning algorithm which can take in an input image, assign importance (learnable weights and biases)

to various aspects/objects in the image and be able to differentiate one from the other. As such, the convolutional layer learns relevant features in the image useful for the computation of noise. The convolutional layer consists of several filters that can be convolved with the input to yield these relevant features. While in

I σ2

Conv.

Downsample

Conv.

Downsample

Conv.

Downsample

Conv.

Downsample

Conv.

Downsample

Conv.

Downsample

Conv.

Downsample

Conv. The convolutional layer requires four parameters; the number of filters, the kernel size, padding, and the stride.

Averaging

Downsampling Layer

The downsampling layer is also a convolutional layer. However, it has a stride of 2. This implies that it learns a downsampling of the input data that encodes the best representation of the input data. As regards the task of noise removal, this layer removes redundancy from the learned representation of the convolutional layer above it. This ensures that only data points that characterize the noise are kept. It also learns a high-level embedding of the convolutional layer above it. This implies that it compresses the information learned in the convolutional layer above.

Structure

In the structure as seen in figure 4.2 and table 4.2, the convolutional and downsampling layers are alternated. The downsampling layer reduces the resolution by half.

This layer learns to extract only the most significant features from the previous layer that contributes to the network ability to minimize the loss. The convolution layer learns to extract features that can be used to characterize noise. Alternating these layers imply that we learn localized noise values at a given layer across the image, then keep only the values considered as important. As compared with the modelbased approach, we hypothesize that the encoder based approach learns a noise distribution within the image, where local noise values are not equally weighted.

For instance, we hypothesize that in regions with partial volume effect, it learns to weight their contributions less.

At the lowest level, we average the values together using an averaging layer.

The output of this layer is the estimated noise σ2 Sample mini-batch of data pairs (I, σ 2 )

3: σ2 = ρ (I; θ E ) 4: L = || σ2 -σ 2 || 1 5: θ ← θ -λ∇ θ L 6: end for
Training of the proposed network allows the model to learn a model that separates the noise component in the tomography image and quantify the variance of the separated noise distribution. This computed variance serves as an estimation for the noise. However, training such model on a 3-dimensional image can be problematic due to the size of these images. Another problem with training such a network is the need for data for which the variance of the noise distribution is known. In reality, this is often difficult to obtain or estimate. To overcome the problem of datasets with a known noise distribution, we create a kalisphera dataset and add a 0 centered Gaussian noise distribution. Kalisphera is an analytical sphere generation tool that can generate a dataset of artificial spherical grains [START_REF] Tengattini | Kalisphera: an analytical tool to reproduce the partial volume effect of spheres imaged in 3D[END_REF]. The generated dataset has certain realistic components obtainable in real dataset such partial volume effect and blur. However, to create such dataset, it is necessary to have to radii and grain positions that correspond to a mechanically stable configuration. We obtain these configurations from Discrete element modeling (DEM) experiments performed in [START_REF] Wiebicke | On the metrology of interparticle contacts in sand from x-ray tomography images[END_REF].

Evaluation

Implementation details

The outlined models were trained using a kalisphera dataset of 7, 500 × 500 × 500 All trainings were performed with the ADAM optimizer [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF], employing a learning rate of 0.0002 and a batch size of 32 samples per GPU. 2 RTX2080Ti GPUs were used. Training was automatically stopped when no improvement was observed on the validation set in the preceding 30 epochs, and the learned parameters leading to the lowest validation loss were selected.

Results

Artificial dataset

The noise estimation using the trained Encoder based neural network shows stronger correlation to the input noise. This can be seen in figure 4.3. At low noise realizations (between 0 and 0.025), the model-based approach often struggles. However, the encoder based approach to determining noise performs well at high and low noise realizations. There is a possibility, however, that the encoder based model overfit to the kalisphera data and cannot generalize well to other datasets. Overfitting is often when the trained model cannot generalize beyond the data on which it is trained.

The significance of overfitting is explored in our evaluation on natural datasets.

Natural dataset

Natural occurring datasets were evaluated using the trained model. Since these datasets were of different sizes (per slice) than the one the model was trained on, two or more slices were often concatenated. We argue that the concatenation approach does not distort the noise distribution as it is expected to be the same from slice to slice. We averaged the noise predictions for the different concatenated slices.

Although it is difficult to evaluate these results due to the fact there is no ground truth for comparison; there exists some trends we consider interesting.

Conclusion

We reviewed the technique in use for the computation of noise in tomography imaging. Using kalisphera, we trained a neural network encoder model to compute noise on tomography images. The computed noise on kalisphera correlates more linearly with the input noise than the model-based approach. This chapter is significant on two accounts;

• The computes the input noise in Kalisphera with better accuracy than the model-based approach.

• It introduces the use of kalisphera data as a means by which complicated "unknown" models can be approximated using neural networks with realistic data for which we have ground truth.

Chapter 5

Image Filtering

Introduction

Noise filtering is an important step in the identification and tracking of granular materials as the process by which we remove noise and artifacts. In effect, noise filtering is an important step for the characterization of granular material. In this chapter, we introduce filters used during this work to process our datasets. We briefly comment on individual filters and the results we obtain when they are used.

A review of methods for image filtering

Linear filters and Convolution

Filters transform an input signal into an output signal. A simple class of filters is the linear filter. By definition, linear filters possess the characteristic, linearity. Linearity encapsulates two fundamental principles; additivity and homogeneity. A noise filter which adheres to these axioms of homogeneity and additivity is termed a linear filter.

Images are discrete; thus, filter implementations are discrete. Discrete linear filters can also be space invariant. This implies that, at each spatial position, the filter composition is constant. As a result of space invariance, the filter transformation operation can be represented by a convolution [START_REF] Ahmed I Zayed | A convolution and product theorem for the fractional Fourier transform[END_REF], [START_REF] Katznelson | An introduction to harmonic analysis[END_REF].

(K σ * I)(x) := R 3 K σ (x -y)I(y)dy (5.1)
Noise removal filters based on these principles are examined with some detail.

Mean filter

The transformation operator also referred to as kernel for the mean filter is a boxcar function. Thus, the output of the transformation point is the average of the value of points, in the interval defined by the boxcar function [START_REF] Zeki | A vision of the brain[END_REF].

K σ = (b -a)V f (a, b; x) = V H(x -a) -H(x -b) (5.2)
On an image, each pixel is replaced by the average of its' neighbors. The neighbors are defined using a neighborhood scheme; defined by a notion of close distance pixels. The neighborhood function is usually a square neighborhood of a given width.

Gaussian filter

The Gaussian filter has a Gaussian kernel. As a result, a point transformation is a weighted average of windowed pixels. The window is defined by the width of the filter. This is, in turn, defined as the scale of the Gaussian function. The weighting decreases with an increase in spatial distance to the central pixel [START_REF] Zeki | A vision of the brain[END_REF].

K σ = 1 2πσ e ||x-y|| 2σ 2
(5.3)

Limitations of linear filters

Edges are very important in tomography image processing because they mostly embed shape characteristics and physical structure properties of granular materials.

These boundaries, hence, characterize the localized granular materials within the image scene. Linear filters do not preserve edges. They filter across edges, without a delimiting notion of edges. The blurring of edges is an undesirable effect, as the erosion of structural properties skews future computations to be made.

Another undesirable quality is edge shifting. Since a Gaussian filter replaces a voxel value by the weighted combination of it's neighbor voxels, there is an edge shift. This implies that the resulting location of the edge is not the exact location of it. This affects the accurate quantification of properties such as the coordination number. It is thus expedient to examine a class of filters that preserve edge structure whilst removing noise.

Edge aware filters

Edge-based filters take extra parameterization; either in terms of image specific edge quantification or other apriori information. This extra parameterization makes them nonlinear and in effect not linear space invariant. Essentially, this means they can hardly be formulated as a convolution operation except approximations are made [START_REF] Paris | A fast approximation of the bilateral filter using a signal processing approach[END_REF]. However, the principle of pixel neighborhood is fundamental to their applicability.

Anisotropic diffusion

Anisotropic diffusion is a partial differential equation (PDE) based filtering technique that models image filtering as a diffusion process. The diffusion equation, a generalization of the heat equation describes changes to image voxel density in images undergoing diffusion over time. This model lays a foundation for the scale-space theory as described in [START_REF] Weickert | Anisotropic diffusion in image processing[END_REF]. It was first introduced by [START_REF] Perona | Scale-space and edge detection using anisotropic diffusion[END_REF].

∂I ∂t = ∂ 2 I ∂x 2 + ∂ 2 I ∂y 2 + ∂ 2 I ∂z 2 (5.4)
∂I ∂t = div c(x, y, z, t)∇I = ∇c.∇I + c(c, y, z, t)∆I (5.5) c(.) is a function that determines the diffusion coefficient. The degree to which diffusion occurs in time is dependent on the diffusion coefficient. Therefore for an edge preserving filter, a diffusion coefficient that is modeled to discourage diffusion at edges and encourage diffusion at flat zones will preserve edge structures.

It models the flux to be dependent on edge strength and thus avoid material changes in such regions, in effect, smoothing out homogeneous regions. This ensures that at points where the edge strength is high, diffusion is minimal. [START_REF] Perona | Scale-space and edge detection using anisotropic diffusion[END_REF] proposed that the following equations be used to model an anisotropic behavior on the diffusion based on edge strength

c( ∇I ) = e ( ∇I /K) 2 and c( ∇I ) = 1 1 + ∇I K
This has proven to be sufficient in [START_REF] Sheppard | Techniques for image enhancement and segmentation of tomographic images of porous materials[END_REF], [START_REF] Andrew M Druckrey | 3D characterization of sand particle-to-particle contact and morphology[END_REF], [START_REF] Druckrey | Research Paper[END_REF] for the removal of noise in tomography images.

Total Variation filtering

Total variation filtering was first introduced by [START_REF] Leonid I Rudin | Nonlinear total variation based noise removal algorithms[END_REF]. It models noise filtering as an energy minimisation problem.

I n = I + σ n
where σ n is white additive Gaussian noise.

To recover the original image from the noisy image, [START_REF] Leonid I Rudin | Nonlinear total variation based noise removal algorithms[END_REF] proposed this minimization problem.

E(.) = min I I n -I 2 2 +λ DI 1 (5.6)
Variations and stable solutions to the minimization has been proposed by [START_REF] Leonid I Rudin | Nonlinear total variation based noise removal algorithms[END_REF],

[59], [START_REF] Chambolle | Total variation minimization and a class of binary MRF models[END_REF], [START_REF] Beck | Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems[END_REF], [START_REF] Zhu | Duality-based algorithms for total-variation-regularized image restoration[END_REF] Median filtering

A median filter cannot be expressed as a convolution. This is because the median filter is based on order statistics of an image patch, and the resulting output voxel of a median filter is not a combination of other voxels within a given window.

This technique has been used extensively to process tomography images [START_REF] Zelelew | A volumetrics thresholding algorithm for processing asphalt concrete X-ray CT images[END_REF], [START_REF] Muter | Improved segmentation of Xray tomography data from porous rocks using a dual filtering approach[END_REF]. It computes a filtered value for every voxel in the image by replacing the voxel intensity under consideration, with the median intensity within its neighborhood.

Median filtering is well suited for the removal of point noise (salt and pepper noise) that arise due to magnetic effects. Noise due to magnetic effects, tend to generate voxel intensities with values near the extremities of the intensity spectrum. This characteristic is exploited by median filtering in its discrimination against isolated extreme intensity differences within a neighborhood. In discrete systems like images, this median can be computed from the histogram h(x, i) centred around position x [START_REF] Ma | Constant time weighted median filtering for stereo matching and beyond[END_REF] This filter, by implication, works best in low noise situations with less than severe distortions to the voxel values.

Bilateral filtering

Bilateral filter is an edge-aware filter which is computed as a weighted average of neighbor voxels in a defined neighborhood. It is similar to Gaussian filtering. The difference exists in the way the neighborhood weights are computed. The bilateral filter's weighting does not only take into account the spatial distance of the neighbors but also the intensity difference.

The origin and formulation of bilateral filtering have mostly been jointly attributed to [START_REF] Volker | Non-linear gaussian filters performing edge preserving diffusion[END_REF], [START_REF] Stephen | SUSAN-a new approach to low level image processing[END_REF], [START_REF] Tomasi | Bilateral filtering for gray and color images[END_REF], [START_REF] Yaroslavsky | Digital picture processing: an introduction[END_REF]. We denote the bilateral as BF [.]

G σ = 1 2πσ e ||x-y|| 2σ 2 B[I] = 1 W q∈S G σs ( p -q )G σr |I p -I q | (5.7)
The normalization parameter ensures that the weights in a given window sum up to 1. The parameter σ s controls the width of the spatial Gaussian and thus determines the weight G σs . G σs assigns weights to voxels based on spatial distance.

A large σ s will thus smooth large features while a small σ s will restrict the filtering to small features. The parameter σ r works in the same way to determine the value of G σr and thus control how voxel differences affect the central voxel. As σ r increases, the behavior of the bilateral filter starts to approach that of a Gaussian filter. This is so because a large σ r relaxes the restriction on the effect of allowable voxel difference.

[72] theorised the relationship between noise σ n and σ r to be σ r = 1.95σ n [START_REF] Sochen | Diffusions and confusions in signal and image processing[END_REF], [START_REF] Barash | Fundamental relationship between bilateral filtering, adaptive smoothing, and the nonlinear diffusion equation[END_REF], [START_REF] Barash | A common framework for nonlinear diffusion, adaptive smoothing, bilateral filtering and mean shift[END_REF], [START_REF] Paris | A fast approximation of the bilateral filter using a signal processing approach[END_REF] have characterised the bilateral filter, its effects and failures, sufficiently. One of the highlighted problems is the huge computational complexity in computing the bilateral transformation of an image.

Guided filtering

The guided image filtering process was designed by [START_REF] He | Guided image filtering[END_REF] to solve some of the highlighted problems. It assumes a linear model. Thus, It maps linearly, an input image, to an output image using a guidance image and a set of weights. The input image can also be the guidance image.

The primary assumption of the guided filter is a local linear model between the guidance image I and the output image Q. It is assumed that Q is a linear transform of I in a window w k centered around the voxel k

Q i = a k I i + b k ∀i ∈ w k (5.8)
(a k , b k ) are linear coefficients, which are constant within a given window w k . This linear model ensures that an edge in I will be an edge in the output.

∆Q = a∆I

The implementation seeks a model that minimizes the distance between the filter output and its input. Mathematically, minimizing a cost function such as 5.9

E(a k , b k ) = i∈w k a k I i + b k -P 2 + a 2 k (5.9)
is a regularising parameter which constrains the model parameter a k from becoming too large.

Solving for a k and b k , we obtain,

a k = 1 n i∈w k I i P i -µ I k µ P k σ 2 + b k = µ P k -a k µ I k
where µ P k is the mean of the kth window in image P, µ I k is the mean of the kth window in the image I and σ 2 k is the variance of the window in the image I. To understand this filter intuitively, we take the simplest case where the guidance image is equal to the image for which filtering is required.

In that case,

a k = σ 2 σ 2 + b k = µ k (1 -a k )
This highlights two dominating processes, a k highlights the edge, while b k performs an edge-aware mean filtering of the image, slowing down to 0 Where the a k approaches 1. The relative importance of the edge is controlled using . The parameter thus performs discrimination between edge and noise.

The output image is thus Guided filtering decouples edge structure other image structures, filter these image structure and couple the structures back into an output image.

Q = āI + b (5.

Brief comments on the application of noise filters to our datasets

Bilateral and total variation filters remove noise extensively. However, they also remove the texture/signal. These filters remove pixel variance and produce piecewise smooth regions. They are useful in applications where texture is not required or when the next operation is a thresholding operation. It also makes them useful for the removal of a ring artifact.

The total variation method controls fidelity by using a weighting term. Of note, is, this weighting term also controls the edge fidelity. This highlights the tendency of total variation denoising to introduce blur into an image, with an increase in weights. This is verified in C.5. We observe that as the weight term is increased, edge shifting increases. Although total variation filtering is capable of filtering high noise contents, care must be taken to factor in the consequent effect on the edges.

Edges are an important part of a qualitative characterization. Humans identify three features as edges: a sharp, step-function like intensity change, a ridge, gradient change, or both [START_REF] Paris | Bilateral filtering: Theory and applications[END_REF]. The bilateral filter is particularly good at preserving step-like edges because the range term averages together all similar values within the neighborhood space domain and also assigns tiny weights to different values on the opposite side of the step. However, due to the effects of partial volume, tomography images possess gradient like edges, whose gradients are steep in high-resolution images and less steep in low-resolution images. This can also result in edge shifting or dislocation and should be avoided where edge fidelity is important. Extensions have been proposed in literature to solve this problem; trilateral filter [START_REF] Choudhury | The trilateral filter for high contrast images and meshes[END_REF] , joint bilateral filter [START_REF] Elad | On the origin of the bilateral filter and ways to improve it[END_REF]. However, the guided filter also solves this problem due to its gradient preserving property discussed in [START_REF] He | Guided image filtering[END_REF]. In figure C.4, the image has a low resolution.

Thus, the combination of partial volume and its low resolution results in bilateral filtering, filtering across edges. Guided filter, however, retains edge boundaries. In Appendix A, we compare the runtimes of noise removal algorithms. In Appendix B, we show the denoising capacity of individual filters as noise in increased in the dataset and in appendix C, we show results obtained by the use of each denoising algorithm.

Conclusions

We summarize our observation from the use of denoising filters;

• In the case of low-resolution images, bilateral filtering leads to edge shifts that distort the structure of grains.

• Total variation denoising and bilateral filtering filter towards a piecewise constant image and thus are good for the removal of noise in images where texture is not important.

• Guided filtering is less severe and conserves structural relationships and thus should be used in cases where texture, contact, and grain edges are of importance.

Chapter 6

Contact based hierarchical segmentation

Introduction

X-ray tomography is used in the study of granular materials and has enabled significant findings. Such studies include the study of particle morphology and its evolution as a result of mechanical loading [START_REF] Stephen | Discrete and continuum analysis of localised deformation in sand using X-ray [mu] CT and volumetric digital image correlation[END_REF], [START_REF] Ando | Experimental investigation of microstructural changes in deforming granular media using x-ray tomography[END_REF]. The accuracy of measured micro-structural properties is sometimes limited by the image analysis methods used on the resulting image volume. Image segmentation methods have been used extensively for grain labeling on the acquired grayscale volume images (in which each voxel is associated with a scalar value) [START_REF] Stephen | Discrete and continuum analysis of localised deformation in sand using X-ray [mu] CT and volumetric digital image correlation[END_REF], [START_REF] Ando | Experimental investigation of microstructural changes in deforming granular media using x-ray tomography[END_REF]. The most notable of segmentation procedures used in tomography imaging for granular materials is morphological watershed [START_REF] Andò | Experimental micromechanics: grain-scale observation of sand deformation[END_REF], [START_REF] Tariel | 3D Multi-scale segmentation of granular materials[END_REF], since these materials can often be imaged with sufficient contrast so that image thresholding can meaningfully be applied to identify the solid The scale of interest for analysis can be the grain scale, or it can be a finer (micro) scale. An example of the use of fine scale segmentation might be to extract the texture patterns in individual grains. Usually, we require a different segmenta-tion for the study of granular materials at each scale. The quality of segmentation is usually subjective. This subjectivity is because segmentation quality relies on the extraction of the required image features. On a given image, this can differ from task to task. This task dependence makes a single-scale segmentation (one resulting segmentation), difficult to evaluate. Morphological watershed is a way by which we get single scale segmentation. In a morphological watershed segmentation technique, the image is transformed into a topological surface (distance map or gradient) [START_REF] Beucher | The morphological approach to segmentation: the watershed transformation[END_REF]. The morphological watershed algorithm floods the topography, starting from pre-defined marker positions until the marker basins meet on watershed lines [START_REF] Beucher | The morphological approach to segmentation: the watershed transformation[END_REF].

Hierarchical segmentation, has become a major trend [START_REF] Deng | Imagenet: A large-scale hierarchical image database[END_REF], [START_REF] Arbelaez | Contour detection and hierarchical image segmentation[END_REF], [START_REF] Ren | Image segmentation by cascaded region agglomeration[END_REF], [START_REF] Arbeláez | Multiscale combinatorial grouping[END_REF] due to its multiscale solution. A hierarchical segmentation is a segmentation technique that results in a sequence of segmentation maps; from a fine scale to a coarse scale (i.e., small to large regions) [START_REF] Guigues | Scale-sets image analysis[END_REF], [START_REF] Jamil | A hierarchical image segmentation algorithm based on an observation scale[END_REF], [START_REF] Santana | Evaluation of Combinations of Watershed Hierarchies[END_REF]. Multiscale segmentation, in this context, is a subset of the study of scale-spaces (algorithmic solutions that span across multiple scales) seen extensively in literature [START_REF] Andrew P Witkin | Scale-space filtering[END_REF]. Consequently, the resulting hierarchical segmentation is not a single partitioning of the image pixels into sets but rather, a multiscale structure comprised of segmentation (partition sets) at increasing scales. Many studies have formulated watershed based on graphs as a hierarchical scheme [START_REF] Soille | Constrained connectivity for hierarchical image partitioning and simplification[END_REF], [START_REF] Cousty | Watershed cuts: Minimum spanning forests and the drop of water principle[END_REF]. This implies that the segmentation of an image can result in a solution at each scale of the image.

The finest scale, i.e., the one with the partition with smaller regions, of the hierarchy can be an initial segmentation or the image pixel set. Such initial segmentation is usually an oversegmentation of the image. The oversegmented image is based on a defined minima set. A set of minimum can be chosen randomly or using a local-minima generation procedure [START_REF] Beucher | Maxima and minima: a review[END_REF]. The minima set can be used along with morphological watershed to generate an initial segmentation.

A sequence of hierarchical segmentation maps can be obtained from such minima based oversegmentation by the successive removal of a minimum. At each iteration, the minimum which ranks lowest is removed. The ranking is generated by an attribute ordering function (extinction function). Attribute ordering functions rank minima based on defined attributes on the resulting region [START_REF] Guigues | Scale-sets image analysis[END_REF], [START_REF] Santana | Evaluation of Combinations of Watershed Hierarchies[END_REF] such as volume or area. The attribute ordering function maps a given segment/region resulting from flooding/segmentation using each minimum from the minima set to a scalar value. The removal of a minimum implies that the image region it generated be merged to the most similar region to it. A similarity measure computes region closeness [START_REF] Jamil | A hierarchical image segmentation algorithm based on an observation scale[END_REF].

Geometric properties of granular materials have been studied widely in literature since some geometrical features can control mechanical, e.g., particle shape, and inter-particle contacts. Granular contact is the characterization of the mode by which two grains are in contact. Some studies of granular contact relate a contact to the size and shape of grains in the assembly [START_REF] Azéma | Force chains and contact network topology in sheared packings of elongated particles[END_REF]. These studies imply that the mode by which grains are in contact is dependent on the size and shape of the grains.

Where grains are similar in terms of size and shape, the nature contact is identical.

Although the nature of contact is known (can be inferred from the physical grains being imaged), it is seldom put into account during the hierarchical segmentation process, and as such, it cannot contribute towards generating meaningful segments.

Thus we propose a method by which hierarchical segmentation can be applied to yield segmentation maps that conform to known contact types. We derive a hierarchical segmentation process that takes known contact types into account. The contact model is added to ensure that our hierarchical segmentation conserves image regions that fit the model.

In this chapter:

• We investigate the use of a contact model in hierarchical segmentation of granular materials. We claim that this leads to better and informative segmentation of granular materials.

• We show how generated datasets using kalisphera [START_REF] Tengattini | Kalisphera: an analytical tool to reproduce the partial volume effect of spheres imaged in 3D[END_REF] can offer insight into algorithmic performance. Kalisphera is a modelling tool by which we can generate images of spherical granular materials. Noise and blur are added to generated data in the style of [START_REF] Wiebicke | On the metrology of interparticle contacts in sand from x-ray tomography images[END_REF].

• We show also that with fractured or broken grains, our formulation can represent them as a single grain. It does this provided there is enough contact between broken particles.

Our chapter organization is as follows; First, we review the fundamentals. Here we review the fundamentals of hierarchies, attribute filtering and saliency maps.

Then we present the proposed granular contact model and show how it can be used with hierarchical segmentation. Next, we present our experimental procedure and set-up. Next, we evaluate the results. We compare the proposed method against hierarchical schemes that do not incorporate a contact model. We also examine different dataset types and scenarios where grains break. To close the chapter, we present our major conclusions.

Definitions

Graphs and notations

A 3-dimensional tomography image can be represented as a weighted digraph G = (V , E, w), whose vertices V are image pixels in the image I ⊂ Z 3 , and edges E are adjacent pixel pairs as defined by an adjacency relation A [START_REF] Soille | Constrained connectivity for hierarchical image partitioning and simplification[END_REF], [START_REF] Santana | Evaluation of Combinations of Watershed Hierarchies[END_REF]. A pixel x is adjacent to a pixel y, if x is in the neighbourhood of y. An edge between x and y is denoted by e x,y .

Partitions and hierarchy

A partition P of a finite vertex set V is a set of nonempty disjoint subsets of V whose union is V . Any element of the partition P of V is called a region of P and represents a region of connected pixels/superpixels in the image. Any pixel x which is an element of the set V , uniquely belongs to an element (region) of P [START_REF] Soille | Constrained connectivity for hierarchical image partitioning and simplification[END_REF], [START_REF] Santana | Evaluation of Combinations of Watershed Hierarchies[END_REF].

This unique relationship is denoted as [P ] x . Given two different partitions P and P of a set V , we say that P is a refinement of P if any region of P is included in a region of P [START_REF] Soille | Constrained connectivity for hierarchical image partitioning and simplification[END_REF], [START_REF] Santana | Evaluation of Combinations of Watershed Hierarchies[END_REF].

A hierarchy (on V ) is a sequence H = (P 0 , . . . , P l ) of indexed partitions of P such that P i-1 is a refinement of P i , for any i ∈ 1, . . . , l. The integer l is called the depth of H.

Saliency maps and ultrametric maps

The cut of P (for the graph G) denoted by φ(P), is the set of edges of the graph G whose two vertices belong to different regions of P. The saliency map of H is a mapping Φ(H) from the edge set E to 0, . . . , l such that each edge is represented by the maximum partition depth λ in which it belongs to the cut set [START_REF] Arbelaez | Boundary extraction in natural images using ultrametric contour maps[END_REF], [START_REF] Arbelaez | From contours to regions: An empirical evaluation[END_REF]. An Ultrametric Contour Map (UCM) is an image representation of a saliency where pixel values corresponds to the mapped values in the range 0, . . . , l map [START_REF] Arbelaez | Boundary extraction in natural images using ultrametric contour maps[END_REF], [START_REF] Arbelaez | From contours to regions: An empirical evaluation[END_REF].

Minimum spanning tree (MST)

The Minimum Spanning Tree (MST) is a subgraph T of G such T is connected, acyclic and includes all the vertices in G. For a graph to be connected implies that for any two nodes in the graph, there is a path between them. Acyclic implies that there are no circular paths. A spanning tree thus represents different edge structure on a connected graph. A graph G can thus have multiple spanning trees. The weight of a spanning tree is the sum of all its connecting edge weights. Thus a minimum spanning tree is a spanning tree T with minimum possible weight.

Hierarchical Segmentation

For each 3D image, a fine partition set can be the 3D pixel points or can be produced by an initial segmentation (as in a set of superpixels). Figure 6.5 shows the Minimum Spanning Forest (MSF)/ fine partition. This fine partitioning contains all the image edges. A similarity measure is defined between adjacent regions of this fine partition set. The superpixels are the nodes/vertexes of the graph structure. Adjacent pixels are connected together by edges with weight w. w is computed according to a dissimilarity measure which we will define.

Starting with a base minima set, the lowest ranked minimum is progressively removed from the minima set according to an attribute ranking function (extinction function). This results in an indexed hierarchy of partitions (H, λ), with H a hierarchy of partitions and λ : H → R + . λ is a function that maps each successive partition to scalar values in an increasing level set. When a minimum is removed, the image region associated with it is merged to the most similar adjacent region (defined by the dissimilarity function). The creation of partition sets, as a function of decreasing minima set, results in a hierarchy [START_REF] Santana | Evaluation of Combinations of Watershed Hierarchies[END_REF].

The hierarchy can be made to emphasize an image statistic. This is achieved by choosing an attribute ranking function that ranks regions based on the chosen statistic. For instance to emphasize objects of uniform size we use the area attribute function. At every level set of the hierarchy (partition set), each region results from a minimum in the level minima set.

On the resulting hierarchy, the saliency map and consequently the UCM is computed [START_REF] Arbelaez | Boundary extraction in natural images using ultrametric contour maps[END_REF]. The UCM is a mapping from the saliency map to an image. Every node's pixel value is mapped to its saliency value in the resulting image. The resulting image shows the persistence of edges in the scale space. Stronger edges persist in the scale-space across multiple resolutions. They are thus more likely to appear in more image partitions across the indexed hierarchy.

Constrained connected components (α-ω)

Extracting relevant partitions at different resolutions in the hierarchy has been actively researched over time. A segmentation tree was used and optimized in [START_REF] Steven | Picture segmentation by a tree traversal algorithm[END_REF]. A stepwise optimization was adopted in [START_REF] Goldberg | Hierarchical segmentation using a composite criterion for remotely sensed imagery[END_REF]. [START_REF] Oj Morris | Graph theory for image analysis: an approach based on the shortest spanning tree[END_REF] extracted relevant nodal connections using a shortest spanning tree segmentation. The pyramid of the region adjacency graph was adopted in [START_REF] Chassery | Géométrie discréte: en analyse d'images[END_REF]. A graph weighted hierarchy was used in [START_REF] Haxhimusa | Evaluating hierarchical graph-based segmentation[END_REF], [START_REF] Scharcanski | Paired graph hierarchies for color image segmentation and representation[END_REF].

Guigues et. al. in [START_REF] Guigues | Scale-sets image analysis[END_REF] introduced hierarchical cuts as a function of minimizing an energy criterion. However, in [START_REF] Soille | Constrained connectivity for hierarchical image partitioning and simplification[END_REF] the idea of using constrained connectivity was introduced. This formulates segmentation at a given resolution as one of finding the right maximal-connected components threshold α on the ultrametric map.

A connectivity criterion is used to extract partitions at a given scale or resolution. The connectivity criterion stipulates that two pixels are connected if there is a path P between them and the difference between successive elements in the path does not exceed a threshold. This was introduced in [START_REF] Nagao | Region extraction and shape analysis in aerial photographs[END_REF] and is linked to single-linkage clustering and minimum spanning trees in [START_REF] John | Minimum spanning trees and single linkage cluster analysis[END_REF].

A path P between two pixels x and y in G or I is a sequence of n > 1 pixels (x = p 1 , . . . , p n = y). Thus two pixels x and y of an image I are α-connected if there exists a path P between them such that range of intensity values between two successive pixels does not exceed a local range parameter α. The α connected components of a pixel x is set of pixels that are α connected to the pixel. It is denoted as CC α (x). On the hierarchical scale-space, varying the value of α varies the resolution of extraction.

An important property of α-connected components is that they form an ordered sequence such that

CC α (x) ⊆ CC α (x) (6.1)
for all α ≥ α.

This implies that the regions that result from CC α (x) is a refinement of CC α (x).

A problem with this extraction strategy is that it suffers from the chaining effect [START_REF] Hambrusch | Parallel algorithms for grayscale image component labeling on a mesh-connected computer[END_REF], [START_REF] Soille | Constrained connectivity for hierarchical image partitioning and simplification[END_REF]. Chaining effect is best explained in 6.1, where although α is 5, the overall range is 9. This is because α is satisfied locally along the path but not globally. In [START_REF] Hambrusch | Parallel algorithms for grayscale image component labeling on a mesh-connected computer[END_REF], an additional parameter ω is introduced to limit the overall variation.

According to [START_REF] Wang | On parameter-dependent connected components of gray images[END_REF], two pixels x, y of image I are (α, ω)-connected if there is a path P between them such that range of successive points is less than α and the overall range is less than ω. However, as is evident and also noted in [START_REF] Wang | On parameter-dependent connected components of gray images[END_REF], [START_REF] Soille | Constrained connectivity for hierarchical image partitioning and simplification[END_REF],

this connectivity relation, although it is symmetric and reflexive, it is not transitive. This definition hence is not optimal. To satisfy these properties, [START_REF] Soille | Constrained connectivity for hierarchical image partitioning and simplification[END_REF], defined the α, ω connectivity as

CC α,ω (x) = max CC αi (x)|α i ≤ α, R CC αi (x) ≤ ω (6.2)
where R is a range function. This implies that the overall range of the connected components must be less or equal to ω.

Minima Ranking

A linear function f : r i → R maps each region r i in the partition to a scalar value.

The functional mapping is based on a defined criterion (volume, area, dynamics) on the region corresponding to each minimum. The scalar attribute values are sorted and the lowest ranked minimum is removed. Its corresponding region is merged with the closest region.

Here we add a contact function to the attribute ranking function. This enables valid granular contacts to persist in the hierarchy space. A single scale segmentation is extracted from the contour maps using the α -ω method described in [START_REF] Soille | Constrained connectivity for hierarchical image partitioning and simplification[END_REF]. We use this method due to fact that it y i corresponds to the target label. p and q correspond to scalar constants. When p is small and finite and q is 1, the formulation corresponds to the random walker formulation. However when p is small and finite and q is 2, graph cuts. As p → ∞, variations in q correspond to a family of watershed variations referred to as power watershed.

Proposed contact based attribute function

The contact-based attribute function requires a separation between foreground and background. This is so because contact interaction between two grains is different from contact interaction between grain and void (background). This implicitly means that we would consider 2-phase materials only. Thus we first separate the background. This separation between background and foreground can be achieved using image thresholding methods such as Otsu thresholding [START_REF] Otsu | A threshold selection method from gray-level histograms[END_REF].

Recall that a cut of a partition set P, is a subset of the edge set E such that both vertexes of an edge, e x,y ∈ E, v x and v y belong to adjacent regions. Consider two adjacent regions in a partition P, r i and r j . The contact set C is thus a partitioning of the cut set, such that ∀c i,j ∈ C, c i,j is a set containing edges between two adjacent regions r i , r j in contact. A region can be composed of multiple vertices and thus we can have multiple edges between two regions.

We propose a contact based attribute function to highlight two contact properties necessary for a contact model in a spherical granular material assembly.

• Length With the background node excluded, the total contact surface area per region surface area should be minimized. Let C i be a subset of C containing all the contacts of a given region r i and all of its adjacent regions in contact.

K(i) = 1 |C i | |c i,j | ∀c i,j ∈ C i (6.4)
For two spheres in contact, the contact surface relative to the surface area should be minimal. Minimizing grain to grain contact length constrains the segmentation to small contact lengths. Large contact length in a spherical grain assembly mostly signifies over-segmentation.

• Flatness Real contact surfaces are mostly approximated by a plane. As such, we penalize contacts for which an approximation by a plane induces large deviations from the plane center (outliers). We compute how closely a contact surface can be approximated by a plane, by looking at the minimum of the contact bounding box dimension. This approximation reflects the prevailing shape of the contact.

We find the enclosing cuboid (w, h, d), aligned with the image axis, for the set of points. The minimal enclosing cuboid is computed as the bounding box on the contact points which is the bounding box on the convex hull of the set of points. We denote flatness of a contact c i,j ∈ C i as;

u(c i,j ) = min(w, h, d)
Thus for a contact set of a region (minima), the flatness score is denoted as, the sum of values u(c i,j ) over the contact set, normalized by the cardinality of the contact set |{u(c i,j ), u(c i,k ), ...}|. This flatness measure of region r is denoted by U (r).

U (i) = 1 |C i | u(c i,j ) ∀c i,j ∈ C i (6.5)
In reality, this approximation might not be sufficient for contacts of grain with more complex shapes. This is because complex (angular) shapes have a larger variance in the type of contact and can thus have larger contact areas and flatness areas. We show in figure 6.4 a wrongly segmented grain having both long contact lengths and curved contacts.

We combine these properties into an affine objective function, allowing larger penalization. This implies that each parameter serves as a scale for the other, thus leading both properties having a larger combined effect in the resulting scalar value for each region. The overall function to be minimized is, calculates the average contact length. Thus its parametric maximum is the surface area of the grain. The parametric limit of U (i) is the radius of a grain. Assuming, a spherical model, we normalize K to refer to its radius, thus ensuring the same scale for both independent measures. This results in equation 6.6.

O(r) = 1
U (i) K(i) (6.6)
Contact attribute descriptor is a soft-feature. This implies, that although it can improve the accuracy of descriptor statistics such volume or dynamics, it is not sufficiently reliable as a stand-alone descriptor for an attribute function as it is non-increasing. In other words, the contact attribute function should be combined with other descriptor functions to improve the accuracy of segmentation. The cost function space of contact-alone based attribute function is more difficult to navigate due to its non-linearity and non-convexity i.e., it does not satisfy the increasingness criterion as defined in [START_REF] Guigues | Scale-sets image analysis[END_REF]. As such it is not monotonic across scales. However, it can aid other attribute functions towards having a well-defined minimum in the cost function space. To combine the contact function with other attribute functions, we design a minima ranking function as seen in algorithm 3. A feature such as volume is computed on the minima set. Each minimum is given a rank. The contact index is also computed with each minimum earning a rank. Both ranks are combined using Where N c is the contact ranking, N m is the ranking from an attribute like volume, area or dynamics and ψ is a weight. The linear combination of the attribute filter with the contact filter yields a smoothing on the space of attribute filter choices.

N r = N m + ψN c
This implies that the contact attribute imposes a bias towards our contact model especially when regions/minima are closely ranked. This is why we mostly adopt ψ ≤ 1.0 As such the overall combination should remain a monotonically increasing or weakly monotonically increasing function. UpdateMSF(minima, segments)

UpdateHierarchy(segments) end while

This combined value N r defines a score for each minimum. We iteratively remove the least scoring minimum. The iterative removal of minimums results in a hierarchy with the saliency of its features (edges) showing how much they appear within the hierarchy. As discussed, the saliency map can be converted into an image representation called the ultrametric contour map. To extract a single segmentation at a given scale, we use the (α, ω) technique to extract clusters. Where α = ω and α, ω ∈ R. Starting from a given pixel point x we find the CC α,ω of that point which yields all pixels points belonging to the cluster. This is iteratively done until all pixels in the UCM belongs to a cluster. The clusters are independent of the starting pixel. The extracted clusters represent the segmentation.

Experimental study

Datasets

Synthetic dataset

We assessed our method by using synthetic images of spheres. The edges of the spheres were analytically represented as partially filled pixels. These partially filled voxels simulate partial volume effect seen in natural images. Kalipshera was used to generate these synthetic spheres [START_REF] Tengattini | Kalisphera: an analytical tool to reproduce the partial volume effect of spheres imaged in 3D[END_REF]. A data set used to build the kalisphera assembly was obtained from [START_REF] Wiebicke | On the metrology of interparticle contacts in sand from x-ray tomography images[END_REF]. It consists of a time-series of 7 mechanically stable configurations of 5522 spheres undergoing compression. These configurations were obtained using a discrete element modeling (DEM) simulation. The data set consists of 5522 granular positions and radii. These positions and radii are then rendered into a 3D image using kalisphera, having only to choose the pixel size relative to the sphere size. The pixel size is 15e -6 m/pixel. The size distribution of the grains is shown in figure 6.5e.

The result is 7 volumes of 500×500×500 pixels with 5522 grains with the grains displaced(translated) within the timesteps. The dynamic range of each image is 8bit, with the background (pores) at 64 and the foreground (grains) at 192. 10% standard deviation of the maximum image value(192) of noise is added to each image. An example of kalisphera generated grains is shown in figure 6.6.

Real dataset

Our method is also evaluated on a granular material imaged using x-ray tomography imaging. Sapphire spheres which are shown in 6.7 is a spherically shaped material.

The imaged sample consists of regular spheres with physical diameters ranging from 300mm to 800mm. The sample image used is a 256 × 256 × 256 volume containing 109 grains.

The algorithm is also evaluated on Leighton Buzzard sand. Leighton Buzzard sand is angular in shape and thus deviates from the spherical model iterated. This dataset is also a temporal dataset. Implying that it contains several 3D images with grains undergoing compression. During compression, grains also break. This in effect provides more size variations along the temporal dimension. We evaluate four temporal sample images with a size of 300 × 300 × 300. A sample image is shown in figure 6.7.

The algorithm is evaluated on images of Zeolite. This dataset contains 5 temporal images. They are spherical in shape and posses visible inclusions in the grains.

The variance between the size of the grains is not wide and thus the sizes can be said to uniform. This is shown in figure 6.7.

Assessment method

To account for performance, we use a supervised assessment strategy outlined in [START_REF] Perret | Evaluation of Hierarchical Watersheds[END_REF]. The likelihood of extracting optimal segmentation is evaluated with respect to the ground truth. The quality of the extracted segmentation is measured using the Bidirectional Consistency Error (BCE) [START_REF] Royal Martin | An empirical approach to grouping and segmentation[END_REF].

Bidirectional Consistency Error (BCE)

The bidirectional consistency error defined in [START_REF] Royal Martin | An empirical approach to grouping and segmentation[END_REF] is used to measure the fidelity of the obtained segmentation. This criterion is defined between 0.0 and 1.0 with 1, indicating a perfect match.

A segmentation map is a partition set P which composed of regions r i . Thus we denote the resulting segmentation as P g and the reference segmentation as P r .

BCE(P g , P r ) = 1 n n i=1
max E(P g , P r , r i ), E(P r , P g , r i ) (

where r i is a region in the partition P g and E is a function computed as;

E(P 1 , P 2 , r i ) = |R(P 1 , r i )|\R(P 2 , r i ) |R(P 1 , r i )| (6.8)
|.| is the cardinality of the set, R(.) is a function that calculates the corresponding region in the segmentation map/ partition set that the region in consideration belongs to. \ denotes the set difference.

Hierarchy Evaluation

The accuracy potential of a hierarchy is computed as the curve BCE vs the fragmentation level. Since BCE measures how closely a segmentation matches the ground truth, the measure of BCE against fragmentation level, shows closely our method approaches an optimal BCE. This implies that it measures how the merge decisions affect the overall optimization.

Fragmentation level is measured as the expected number of segments over the actual number of segments in the reference segmentation map.

f ragmentation = |P g | |P r |
|P r | refer to the number of segments in the ground truth image while |P g | is the number of segments gotten. This measures the tendency of a hierarchical process to make the right merging decisions leading towards the best possible segmentation; close enough to the ground truth segmentation.

Experimental set-up

Our hierarchical watershed is compared against a morphological watershed on the kalisphera dataset. For hierarchical segmentation, a single scale segmentation is extracted by varying the α, ω parameter as iterated in [START_REF] Soille | Constrained connectivity for hierarchical image partitioning and simplification[END_REF]. The values of α, ω is progressively increased (α and ω always have the same value), till the image is a single label or region. For each extracted segmentation, we compute the BCE score and the fragmentation level.

Two variants of the morphological watershed are considered; gradient-based and distance map-based. In the gradient approach, we compute the gradient magnitude of the 3D image. The minima on the gradient magnitude surface are computed and both are used by the watershed transform to generate a label image. The distance map approach starts by thresholding the image using Otsu thresholding [START_REF] Otsu | A threshold selection method from gray-level histograms[END_REF].

A distance map is computed on the threshold image using the Euclidean distance map function. The minima of the inverted distance map are used by the watershed function to generate also a corresponding label image. The connectivity used in both gradient and distance map computations is the 26-neighborhood connection (fully connected). To have the same basis for comparison, we use a masking layer given by the thresholding (otsu thresholding)to separate foreground from background in morphological approaches. This implies that we use the separated background/foreground used in the hierarchical schemes as an overlay for the images, thus having the same starting point which is a properly segmented background.

We also verify the effect of combining a contact based function with known attribute functions. Thus we compare the result accuracy between when these attribute functions are used alone and when they are combined with a contact model.

We use the combination strategy outlined in equation 6. thus leads to a cascade of better decision making steps. This, in turn, leads to image structures persisting across many scales. This proves a significant improvement over dynamics or volume alone as a metric for attribute filtering.

Comparison against morphological watershed

The segmentation algorithm predominantly used in geomechanics is the morphological watershed algorithm. To have an overview of the improvements our model provides, we compare the proposed method against it. We aim to quantitatively evaluate the performance and merits of our method. Our method is compared to two morphological watershed variants; gradient and distance map based. We compare these methods on the Kalisphera dataset. We iterate through the hierarchy (varying α, ω), calculating the best BCE value from the hierarchy for each α-ω pair.

We compare this BCE against BCE values obtained using morphological watershed variants. Results are shown in table 6.11b. Results of best BCE values against morphological and power watershed averaged over 7 time-steps is shown in 6.9.

We observe that our proposed method performs better than all the examined methods. It is also shown that the gradient-based morphological watershed performs the least. It is argued that the amplification of noise by gradients mostly leads against hierarchical based methods using the defined contact model to poor segmentation of geomaterials. Gradient induced problems coupled with partial volume effect is argued to be the reason for its poor performance. This has informed the trend of the use of distance maps in the processing of images of granular materials. The distance map based watershed is also affected by noise. But, it is only affected as long as the image noise is large enough to distort the thresholding of such image. However, it still leads to better performance than the gradientbased approach. Our dynamics hierarchical segmentation plus contact model has the best performance. This is so because it has significant descriptive power as it uses both the topological landscape and contact prior. This leads to significantly better performance. We argue that more descriptive statistics leads to better segmentation.

We also compare our method against the power watershed algorithm. We vary p and set p = 20. We observe better performance than the power watershed approach.

We explain this variance in subsection 6.5.4.

Evaluation of unsupervised hierarchies on Sapphire Sphere grains

In analyzing the effect and usefulness of the proposed method, we evaluate its performance on synthetic materials; Sapphire spheres imaged using x-ray tomography at the Laboratoire 3SR. Due to the subjective nature of ground-truths on this dataset,

we qualitatively analyze the results obtained to visually understand the effect of the proposed method. In the ultrametric contour map (UCM), we observe in figure 6.13 that weak gradients at grain contact points are more resistant to being merged prematurely, thus resulting in higher saliency values. This reduces the possibility of undersegmentation. This is shown in figure 6. [START_REF] Peters | Principles of Computerised Tomography[END_REF].

This shows that the segmentation merging process evolves in such a way as to find a contact configuration that is realistic in terms of the proposed model. It is, therefore, less attracted to oversegmentation or grain configurations with contacts that do not fit the defined model.

To further verify this phenomenon, we observe the segmentation of two grains in contact. As shown in figure 6.12, we observe that when a dynamics contact extinction is applied, the best obtainable segmentation under segments both grains, merging them as one due to the weak contact between them. However, a combination with the attribute contact model yields a clear separation between both grains in the extracted segmentation. The combination clearly improves the ability to segment grains in contact with weak delineating gradients. These types of gradients are common in 3D x-ray tomography images due to low imaging resolution and the effects of partial volume.

Comparison against power watershed

To validate the performance of the proposed method against graph-based techniques, we compare its performance against power watershed.

We applied the power watershed algorithm to sapphire spheres and zeolite grain assemblies. We observe that provided the gradient magnitudes within a grain are less than the gradient magnitude between void and grain, there exists a configuration in p and q for which power watershed performs just as well as the contact based segmentation. This can be seen in figure 6. [START_REF] Rosenblatt | A central limit theorem and a strong mixing condition[END_REF] We argue that this so because as p → ∞, the formulation becomes biased towards the preservation of edges and is less perturbed by weak edges. The justification for the use of our contact based method in light of the performance of the power watershed family would simply be the ease by which we can navigate the scale-space. This implies that we are able to navigate through a variety of segmentation using contact based segmentation and choose the optimal segmentation for the specified task.

Preliminary conclusions

• The proposed method outperforms morphological watershed variants in the segmentation of spherical and near-spherical grains.

• The proposed method biases the merging towards conservation of contacts as the edge weight criterion using gradients might not be enough due to poor imaging resolutions or partial volume effect.

• In the observed images there is a configuration in power watershed (p, q) that gives almost similar results in contact preservation.

Qualitative evaluation on Leighton Buzzard sand

The contact interaction in the Leighton buzzard assembly deviates from the characteristic contact interaction observed with non-spherical data due to its angular shape. As already stated, the contact mode is highly dependent on the shape of the material. In assemblies where there is no consistently defined shape such as angular materials, there exists a wide variance in the contact types observed. Some of the contact types observed in the assembly are shown in figure 6.15. The complexity of finding an optimal labeling is augmented by the fact that some grains are fractured or broken. In figure 6.15a, we observe a flat contact where the contact Surface between two grains is larger than that obtained as with spherical grains. Such large contact area violates the assumptions for our formulation in which the area of contact is minimized. As such this would most likely result in undersegmentation. The contact type is shown in 6.15b resulting from the fracture of the grain violates this assumption even more. Broken grains can even violate the Flatness assumption depending on the mode of fracture. Fracture modes would be explored in detail in the succeeding chapters.

Our algorithm is applied to 3D tomography images of Leighton Buzzard sand with different ψ values.Recall that ψ is the weighting signifying the contribution of the contact extinction. Contact extinction weighted by ψ was paired with the dynamics extinction since we have established that a combination of both yields better descriptive statistics. In figure 6.16, we observe that to optimally segment these contact types, it is necessary to limit ψ to low values. This is consistent as we expect a huge dependence on contact extinction to drive the segmentation towards merging flat contacts and in effect lead to undersegmentation. However, this behavior is necessary for grains where breakage occurs (see figure 6.17), as this drives the segmentation to identify broken grains as a single grain as opposed to identifying it as multiple grains. The segmentation of multiple contact points also violates the assumptions made in our formulation due to the fact that multiple points of contact can affect the shape of contact. The resulting contact bounding box for which we deem flatness necessary can thus violate the flatness criterion. This is shown in figure 6.18. Here we observe that a higher ψ value results in undersegmentation of these grains.

Getting our algorithm to work at the sample scale is more complicated due to the complex interdependence between shape, contact and breakage and how they relate to the ψ parameter. We show that it can be carefully tuned to capture broken grains, flat contacts, and multi-point contacts when they occur. In figure 6.19, we observe that the application of our contact model results in the segmentation of

Qualitative evaluation on zeolite undergoing breakage

With the understanding that our formulation augments the accuracy of segmentation and also the probability of associating broken grains as a single grain (this is important to track the evolution of breakage), we apply the formulation to a nearspherical natural occurring grain undergoing "simple" fracture. Simple fracture in this regard is characterized as the splitting of grain into two or three sub-particles, with contact between the sub-particles. As zeolite grains are near-spherical, we evaluate how well our algorithm is able to take advantage of the sphericity assumption and segment a broken grain without assigning different labels to the sub-particles.

We observe in figure 6.20 that with a ψ = 0.38, a couple of broken grains previously assigned labels as different grains are now assigned labels depicting them as a single grain. In the zeolite image examined, the degree to which broken grains are captured is augmented from 22.4% to 68%. However, it falls apart when fracture patterns do not yield particles that are in contact or that have a significant contact area.

Conclusion

The combination of a contact model with dynamics and volume extinction functions is shown to increase the accuracy of grain segmentation in 3D tomography images.

saliency values of these contact regions are augmented, thus making them more resistant to merging. Although the contact prior is simple and can be said to apply to spherical grains, the improved result validates the combination with a contact model yields better segmentation results in multi-scale segmentation

We recall some of the conclusions made in this chapter

• Our method outperforms morphological watershed in the detection and preservation of valid granular contacts, thereby reducing oversegmentation.

• Provided there is no strong texture or noise within granular materials, there is a parameterization for the power watershed family that would yield contact preserving segmentation.

• Using the described method on angular materials with variance in contact topologies might result in sub-optimal results due to the assumption of a spherical contact.

• In spherical grains or near spherical grains with fracture or breakage, our method detects fractured pieces like a grain provided there is contact between the grain particles.

Chapter 7

A study on threshold based breakage detection

Introduction

In this chapter, we examine methods by which we detect broken grains in 3D images of granular materials. The ability to detect grains undergoing breakage is relevant for parametric studies on granular materials. These parametric studies establish links between theoretical models in geomechanics and experimental models by validating or invalidating theoretical claims. As such, the ability to detect grains that break is necessary to understand theoretical models that characterize breakage and breakage of grains in grain assemblies. An example of studies that benefit from the ability to detect broken grains is the study of the interaction of granular materials under high loads.

The detection of broken grains is not trivial. It is so because effectively characterizing breakage in an assembly of granular materials involves processing the resulting 3D sequence of images, labeling the grains, correlating these grains across time-steps and attributing grain particles to parent grains when breakage occurs.

This sequence of steps implies that we perform both spatial correlation and temporal correlation of particles.

In [START_REF] Guida | Breakage mechanisms of highly porous particles in 1D compression revealed by X-ray tomography[END_REF], [START_REF] Karatza | Study of temporal and spatial evolution of deformation and breakage of dry granular materials using x-ray computed tomography and the discrete element method[END_REF], methods were formulated to detect grains; however, these methods have their failings. In this chapter, we highlight the problems associated with these approach to breakage detection proposed in [START_REF] Guida | Breakage mechanisms of highly porous particles in 1D compression revealed by X-ray tomography[END_REF], [START_REF] Karatza | Study of temporal and spatial evolution of deformation and breakage of dry granular materials using x-ray computed tomography and the discrete element method[END_REF] by showing the large number of influences on choosing the right threshold for the isolation of broken grains.

The studies in [START_REF] Guida | Breakage mechanisms of highly porous particles in 1D compression revealed by X-ray tomography[END_REF], [START_REF] Karatza | Study of temporal and spatial evolution of deformation and breakage of dry granular materials using x-ray computed tomography and the discrete element method[END_REF] highlight a process for the detection of breakage in 3D granular assemblies. The method presented detects breakage by examining the value of the normalized correlation coefficient of correlated grains between a time-step where breakage does not occur and a time-step where breakage occurs. Correlation in this context refers to the algorithmic procedure in matching identified grains by their texture a subsequent time-step. In a time-series (stack of 3D tomography images), a previous time-step image is segmented using any of the segmentation methods, and a correlation is performed. The normalized correlation coefficient is a scalar value that specifies the degree of matching for an identified grain. It is usually between 0 and 1 with 1 signifying a perfectly correlated grain. The normalized correlation coefficient of grain can be computed using equation 7.1.

m n [f (m + i, n + j) -f ][g(m, n) -ḡ] m n [f (m, n) -f ] 2 [g(m, n) -ḡ] 2 (7.1)
In the method, they impose a threshold on the normalized correlation coefficient values. Consequently, grains with values below the threshold are characterized as broken. In [START_REF] Karatza | Study of temporal and spatial evolution of deformation and breakage of dry granular materials using x-ray computed tomography and the discrete element method[END_REF] a threshold of 0.96 on 3D tomography images of zeolite grains while [START_REF] Guida | Breakage mechanisms of highly porous particles in 1D compression revealed by X-ray tomography[END_REF] imposed a threshold of 0.94.

One problem reported as a result of this method is the arbitrary method in choosing the threshold value. In [START_REF] Guida | Breakage mechanisms of highly porous particles in 1D compression revealed by X-ray tomography[END_REF], it is shown that a slightly higher threshold value misses out on detecting some broken grains, and a slightly lower threshold detects texture variance as breakage. It is thus imperative to tune this parameter per breakage detection operation. Per-image tuning can be expensive in terms of memory and time since there is no guide to choosing the right threshold. In this chapter, we attempt to find variables by which guide the choice of a suitable threshold for the method described in [START_REF] Guida | Breakage mechanisms of highly porous particles in 1D compression revealed by X-ray tomography[END_REF]. We hypothesize that there is a relationship between the threshold value choice and parametric variables in the image such as random noise, breakage gap, internal porosity, and the number of broken particles. Thus we attempt to find relationships between these variables and the task of choosing the right threshold for a given detection operation. Our analysis is done using kalisphera grains.

Analytical study

To study the effects of the outline parameters on the optimal threshold choice, it is important to be able to control the way grains breakage. In multiple particle tests, it is almost impossible to define breakage patterns. Although it is less difficult when performing single particle breakage tests, it still is not trivial. To study these

Fracturing kalisphera grains

Creating a kalisphera grain using the kalisphera tool is done by defining a grain radius and grain centre. The kalisphera tool then generates the grain with the given parameters. Breaking such a grain into two parts can be simple to achieve analytically. However, complex breakage patterns are not. To break such grain, we approximate breakage as comprising a breakage centre, breakage gap size and rotation (breakage orientation).

The breakage center is the approximate center of the breakage plane to be imposed. The breakage gap size is the distance between two particles are breakage and the breakage orientation is the orientation of the plane of breakage. Implicitly we approximate that every complex breakage pattern is a sequence of splitting events.

Thus to impose a complex pattern, we recursively split the grain and its particles.

The result of every breakage operation is a left and a right component which can further be broken. Based on this representation, we stack breakage operations on a grain and its sub-particles recursively and in effect can create complex breakage patterns. These patterns are as shown in figures 7.1a, 7.1c, 7.1e.

Imposing texture on kalisphera grains

The localized texture is also a characteristic of grains, which is useful in correlation operations. We impose local texture on kalisphera grains by adding Gaussian noise to particles separated using the breakage algorithm. To add local texture to a region of the grain without breaking the grain, we set all breakage gaps to 0. Textured regions are identified by applying the successive splitting operation, however, without the breakage gap. We show an example of a locally textured grain in 7.2.

Results and discussion

Size ratio of breakage particles

To validate the effect of the size ratio between the particles and the normalized correlation coefficient, we created a 3D kalisphera grain of 30 pixels in radius. A corresponding temporal displaced image was created with a displacement of 1 pixel in the x,y, and z axes. The displaced kalisphera grain is artificially breakaged into To validate this hypothesis, we create a kalisphera grain with localized texture on the smaller sized particle. Then we observe the values of locally textured grains vs.

when the texture is uniformly distributed (as in zero texture) as shown in figure 7.2.

The results in figure 7.4, show a difference in the variation of normalized correlation when there is localized texture, and when there is no localized texture. Unique texture patterns contribute significantly to the ability of correlation to recognize a grain. When grains with these landmark patterns break, the effect on the correlation coefficient is stronger than in homogeneous or textureless grains. This phenomenon is the case also in grains with internal porosity as the internal porosity is viewed as texture. This reason can be the probable reason why highly porous Leica grains in [START_REF] Guida | Breakage mechanisms of highly porous particles in 1D compression revealed by X-ray tomography[END_REF] need a significantly lower threshold than Zeolite grains in [START_REF] Karatza | Study of temporal and spatial evolution of deformation and breakage of dry granular materials using x-ray computed tomography and the discrete element method[END_REF].

Number of breakage particles

To account for the effect of the number of breakage particles on the normalized correlation coefficient, a 3D kalisphera grain with a radius of 30 pixels was artificially broken into 2, 3 and 4 particles. A uniform displacement of 1 in all three axes is imposed and each particle is translated. In the case of the image volume with 2 particles, we progressively increase the breakage gap from 1 pixel to 10 pixels. For the other volumes with more than 2 particles, the breakage gap corresponding to the first splitting event is progressively widened. The resulting graph in 7.5 highlights the dependence of the normalized correlation coefficient on the number of particles.

The splitting events resulting in more than 2 particles correlate poorly proportion to the number of particles into which they are split. In effect, grains with higher number of particles have lower normalized correlation coefficient.

Grain size

We also investigate for the effect of grain size on the resulting normalized correlation coefficient. 3D Kalisphera grains with radius ranging from 20 to 90 pixels are created and split along the centre and the breakage gap is progressively widened. It is size. We argue that this is due to the fact that smaller sized grains have less effective texture which is a basis for correlation. This highlights a difficulty in optimally setting the threshold parameter, as grain assemblies usually contain grains which vary in size.

breakage gap

In figures 7.5, 7.6, 7.3, it is observed in all experimental procedures that an increase in breakage gap, always yields a corresponding decrease in the normalized correlation coefficient. This implies, that in any breakage scenario (splitting, attrition or abrasion), an increase in the gap between particles induces a decrease in the normalized correlation coefficient by a factor that corresponds to the dominating environmental factor (grain size, breakage centre and number of resulting particles).

Conclusions

The combined effects of breakage gap, texture, grain size and breakage particle size makes choosing a suitable threshold very difficult. This difficulty is amplified by the fact that grain assemblies posses texture variations and size variations. These variations go on to affect the process of breakage thereby resulting in multiple breakage scenarios with variational differences that hitherto cannot fully be captured by a scalar threshold parameter. In cases where the assembly is homogeneous in size, shape or texture, the ability to capture fully all breakage cases using a threshold parameter exists.

Furthermore, the threshold parameter based detection is not well suited for progressive breakage since the original unbroken state should be correlated to the broken state. As breakage approaches crushing and particles disperse, breakage becomes difficult to characterize.

Chapter 8

Detecting Breakage To formalize the proposed method, we define the parameters involved. Let a tomography time-series be a set of images denoted as S containing images denoted as S i , with i representing the position in the time-series such that i ∈ [1 . . . T ]. T is the sequence length of the image sequence S i.e., the number of images in the time-series. Let a segmentation of S i be a mapping G(.) such that S i is transformed into a labelled(segmented) image L i such that each pixel maps into a label in the label set. Here we define a correlation C to be a mapping from one labeled image L x to another labeled image L y such that every label in L x is associated to a label in L y . Thus the method of back-correlation for the detection of breakage maps L i+1

to L i such that when two labels in L i+1 map to the same label in L i , the grain is considered broken. As is observed, this is dependent on the quality of segmentation.

We denote the sequence of labeled images as L.

However, this method of breakage detection can be computationally intensive since it performs a correlation on each patch. This is sometimes more than is necessary, especially in cases where a few grains break. Therefore we propose an approximation that uses minima to detect broken grains.

Minima based back-correlation

The detection of breakage requires that the particles of a broken grain are correlated in time to a point where the grain in unbroken and that these particles are independently identified to enable their tracking as they undergo more changes in time.

In the minima based detection, we find a coupling between temporal minima and in effect, correlate grains from their unbroken to their broken states. Consider an image sequence S. For S i ∈ S, there is a set of minima M i which maps to a labelling L i . We define a new correlation function C m which maps the set of minima M i+1 to set of pixels M corr i+1 in S i . If S i has been sufficiently labelled and transformed to L i , we propose that M corr i+1 can inform the labelling of S i+1 since its labels are known due to the labelling of L i . A grain is thus broken where two minimum couple to the same label value. This is approximately equal to back-correlation based detection but with correlation done using only the minima and not the entire particle.

Combined back-correlation

The minima based approach bases the correlation of the entire particle or grain on a point. This is prone to errors in cases with high breakage events. As such, we combine both schemes to characterize the evolving breakage of a sample. In the first temporal half, breakage is detected using the minima based approach and in the last temporal half by full back-correlation approach.

Algorithm 4 Combined back-correlation based detection

Visualising Broken grains

Breakage capture refers to the representation of a grain, although broken by one label. However, to detect the grain as a broken grain, it is necessary to map the initial representation of a broken grain (by multiple labels) to its unique representation as one label in the captured or revised state. As such, we correlate its representation before the application of the back correlation algorithm to its representation after the breakage detection algorithm has been applied. Where two particles in the precorrection state map to one particle after its application, we characterize such a grain as broken or "further-broken" depending on the timeline of progressive breakage which we aim to characterize. An example of such characterization is shown in figure 8.3.

Rectifying segmentation using spatiotemporal correlation

A grain is characterized as broken when two segments (particles) in the next temporal image correlate to one label in the previous time-step. Although this captures most breakage events, it also characterizes oversegmentation as breakage in cases where the grain in the next time-step is over-segmented. This amounts to some false positives in the detection of breakage. One way to avoid this is to create better quality images through pre-processing (filtering). We, however, try to detect these cases and use it to perform an informed segmentation and re-labeling of grains. For the informed segmentation procedure, we detect oversegmentation by looking at breakage events and characterizing the contact between the particles of potentially broken grains. Where a contact is characterised as being the result of oversegmentation, the particles are joined together.

In order to perform this characterisation, let the length of the contact between two labels L i L j be |C ij |. Where C xy is the set of boundary points. Also let the three principal axes be p x , p y and p z . Oversegmentation is defined here as

|C ij | ≥ min(p x , p y , p z ).γ (8.1)
Where γ is a scalar parameter which should at least be greater than 0.5, we set gamma to be 0.75. This is found to be sufficient in distinguishing actual breakage cases from oversegmentation.

Experimental Setup

The combined backward correlation breakage detection algorithm is applied to a 3D time sequence of zeolite grains undergoing deformation and breakage and Leighton

Buzzard sand under compression. The zeolite images were obtained from experiments performed in [START_REF] Karatza | Study of temporal and spatial evolution of deformation and breakage of dry granular materials using x-ray computed tomography and the discrete element method[END_REF]. The grains are near spherical in shape and posses impu- The images are filtered using a bilateral filter to remove noise components. In the course of our experiments, the segmentation algorithm used is a morphological watershed, and image correlation was performed using SPAM The grains are first separated using morphological watershed; then broken grains are detected using the combined back correlation algorithm with the first half of samples matched using their minima and the remaining time-steps (more breakage events) correlated using the full back-correlation approach. The labels of successive time-steps are corrected to match the ancestor labels.

We also apply the combined backward detection to 3D samples of Leighton Buzzard sand varying in time and undergoing progressive breakage. The sample

Results and Discussion

Between time-steps 1 and 2 of the zeolite dataset we observe mainly splitting and abrasion breakage events. These breakage patterns are, however, simple enough to be sufficiently captured by the back-correlation algorithm, as shown in 8.6. We characterize the successful capture of a broken grain as when it is assigned the same label as its unbroken parent grain. These breakage patterns are relatively easy to capture due to the compact and non-dispersed breakage modes. This means that the resulting broken particles are either close together as in the parent grain or partially connected. Also, as there are limited events of breakage simultaneously occurring, it is even simpler to identify breakage events. The initial time-steps of Leighton Buzzard sand also present simple breakage cases owing to their large sizes due to the resolution at which they are imaged. Figure 8.7 shows the capture of these breakage cases.

As we progress over time, the breakage modes increase in complexity. In figure 8.8 we analyse the results of detecting breakage on the succeeding time-step i.e from time-step 2 to time-step 3 of zeolite dataset. More complicated breakage patterns are observed which are correspondingly more difficult to detect and capture. As such false positives are thus more evident in the results, we show. We argue that as progressive breakage approaches its ultimate distribution (see [START_REF] Einav | Breakage mechanics-part I: theory[END_REF]), the particle space becomes denser and as such grayscale pixel variance approaches 0. This implies that as crushing occurs within an assembly of grains, the unique granular texture This texture dependence is evident in the result difference between the Leighton buzzard sample and the zeolite sample. The reduced resolution at which zeolite is processed (and the number of grains) ensures that the accuracy of capture is less when grains appear crushed. Thus our method performs better in processing Leighton Buzzard sand than zeolite as we go further in the breakage timeline.

This correlates well with results obtained from texture experiments performed on Kalisphera grains in the previous chapter. Evident also in the succeeding time-step 8.9 is a slight increase in the number of false positives. Leighton Buzzard sand, however, has a higher volume per grain.

As such and consistent with observations on kalisphera grain experiments, size and texture aid correlation. We observe that these elements (size and texture) determine the quality of breakage detection.

Segmentation Rectification

Using the defined rectification algorithm in subsection 8.1.4 we show that spatiotemporal information from the back-correlation algorithm can be used to detect oversegmentation in temporal images.

In figure 8.14, we show a Leighton buzzard grain oversegmented initially due to a hole in the grain. The oversegmentation is rectified during the breakage detection process. Although the process of segmentation rectification reduces the error rate of detection and increases the quality of segmentation, it is computationally intensive due to the computation of an Eigen decomposition on each grain. Although iterative approximations to the Eigen-decomposition exist using the power law [START_REF] Chatelin | Eigenvalues of Matrices: Revised Edition[END_REF], we do not research further.

Chapter 9

Conclusions and Perspectives

The principal objective of this Ph.D. thesis was to explore techniques for the identification and tracking of granular materials undergoing breakage. The task of grain identification and tracking is directly hinged on the quality of data on which algorithms are built. Notably affecting the quality of data is noise. A large part of this thesis is devoted to studying the quantification and the removal of noise to augment the quality of data on which grains are identified and tracked. First, we tackle the task of computing the amount of noise in a tomography image. We reviewed the technique in use for the computation of noise in tomography imaging. Using kalisphera, we trained a neural network encoder model to compute noise on tomography images. The computed noise on kalisphera correlates more linearly with the input noise than the model-based approach. This chapter is significant on two accounts;

• The computes the input noise in Kalisphera with better accuracy than the model based approach.

• It introduces the use of kalisphera data as a means by which complicated "unknown" models can be approximated using neural networks with realistic data which we have a ground truth.

This chapter also goes further to present the failings of our trained encoder. Due to the fact that it was trained using spherical data, it does not generalize to nonspherical data. Means by which this can be mitigated include the DEM approach where complex grains can represented as clusters of spheres. A wider variety of shape factors can in effect be simulated.

In chapter 5, we present grain filtering algorithms we used and how they can be applied to granular materials. We briefly discussed our observations in using each of the filtering algorithms.

In chapter 6, we explore the segmentation (identification) of granular materials. We introduce hierarchical segmentation for granular materials as against the commonly used morphological watershed. We then propose a method by which contact structure information can be used as apriori information to guide the hierarchical segmentation process. This chapter shows that the introduction of a simple spherical contact model drives the segmentation to extract a labeling that is more mechanically consistent with expected contact. This in effect reduces the amount of oversegmentation in tomography images of granular materials. Subsequently we

show that such contact model used for the segmentation of zeolite grains (near spherical) induces the segmentation to label broken grains as a single provided they are in contact. These are reasons why this is significant;

• In a single image, with no temporal context, the task of determining if a grain (whole) is broken is much less trivial than determining if a broken particle is a part of a broken grain. The former can be established by looking at the grain datapoints for a plane of breakage. It can also be achieved by checking for a disparity in segmentation when a morphological watershed is used and when a contact based hierarchical segmentation is used.

• In a temporal image with grains undergoing displacement, rotation and breakage, it is easier to separate (in an image) a clustered broken grain, than to aggregate broken particles into one frame. The context for separating a grain is the grain, however the context for aggregating particles into a grain is the entire spatial image.

However the conceived model for biasing the segmentation towards labeling grains in view of ensuring realistic contacts is spherical and thus works best on spherical or near-spherical grains. Future work in this regard would be to take into account more complicated contact types and shapes.

Chapter 7 studies the only known method by which breakage is detected in granular materials. The method reported in [START_REF] Guida | Breakage mechanisms of highly porous particles in 1D compression revealed by X-ray tomography[END_REF], [START_REF] Karatza | Study of temporal and spatial evolution of deformation and breakage of dry granular materials using x-ray computed tomography and the discrete element method[END_REF] imposes a threshold on the correlation coefficient (degree to which temporally displaced grains match). Grains that fall below the imposed threshold are characterized as broken. A problem that exists with this method to breakage characterization is how to the set the threshold parameter. In this chapter we list potential influences on the threshold parameter and study the effects of these influences with single particle characterizations using kalisphera. To study the effect of these influences we implement a recursive method by which these analytically derived spheres can be broken and by which localized texture ( only a section of the grain) can be imposed. In this study the following Structured noise and texture contribute to an inherent non linear noise response.

And in these cases and at high noise, the performance of filters with linear descriptors degrade. It is noted however, that the bilateral filter performs consistently across all datasets.

conclusions

• Not all noise in images can be characterised as a gaussian.

• Filters behave differently in different texture scenarios.

• More experimentation is needed to validate the induced assumptions.

8. 1 . 4

 14 Rectifying segmentation using spatiotemporal correlation . . . 8.2 Experimental Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.3 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . 9 Conclusions and Perspectives A Algorithmic speed and efficiency of noise removal algorithms B Noise response in natural datasets C Denoising of dataset images List of Figures 1.1 Problem statement figure . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 A typical x-ray system for imaging materials . . . . . . . . . . . . . . 2.2 Relationship between total breakage and breakage potential [21] . . . 2.3 Relative breakage according to Einav in [77] . . . . . . . . . . . . . . 2.4 Different breakage types [94] . . . . . . . . . . . . . . . . . . . . . . . 2.5 Detecting breakage in zeolite images of granular materials undergoing compression [141] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.6 Detecting breakage in images of light-expanded clay aggregate (LECA) undergoing compression [138] . . . . . . . . . . . . . . . . . . . . . . 3.1 Laboratoire 3SR x-ray scanner . . . . . . . . . . . . . . . . . . . . . . 3.2 Dataset slices(SSFK) . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Dataset slices(Zeolite) . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Dataset slices(Leighton buzzard) . . . . . . . . . . . . . . . . . . . . . 3.5 Dataset slices(Highly Decomposed Granite) . . . . . . . . . . . . . . 3.6 Dataset slices(Caicos Ooid) . . . . . . . . . . . . . . . . . . . . . . . 3.7 Dataset slices(Hostun sand) . . . . . . . . . . . . . . . . . . . . . . . 3.8 Dataset slices(Kalipshera) . . . . . . . . . . . . . . . . . . . . . . . . 3.9 Dataset slices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 Progressively adding noise to kalisphera images . . . . . . . . . . . . 4.2 Model diagram for the noise quantifying encoder. Conv. refers to convolution. . . 4.3 Noise estimation in Kalisphera datasets using a Model based approach and a Neural network based approach . . . . . . . . . . . . . . . . . . 6.1 7 × 7 image and its partitions into α-connected components [80] . . . 6.2 7 × 7 image and its partitions into α -ω-connected components [80] . 6.3 Kalisphera grain experiments . . . . . . . . . . . . . . . . . . . . . . 6.4 Oversegmentation of two zeolite grains using regular hierarchical watershed . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.5 kalisphera image of artificially generated grains . . . . . . . . . . . . 6.6 The first three timesteps of the generated kalisphera dataset shown in a, b, and c. The corresponding images with noise added is shown in d, e, and f. The images are thresholded using otsu threshold(g, h, i) and an oversegmentation map is generated for each . . . . . . . . . 6.7 Slice and 3D image showing Leighton buzzard sand depicted as LB, cropped cross section of Sapphire spheres depicted as SS and cropped Zeolite. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.8 Fragmentation plots for volume and dynamics . . . . . . . . . . . . . 6.9 BCE-Fragmentation curves showing the performance against watershed variants. The results are averaged over 7 Kalisphera Images . . . 6.10 BCE-Fragmentation curves averaged over 7 Kalisphera images . . . . 6.11 Natural dataset . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.12 Two sapphire grain segmentation . . . . . . . . . . . . . . . . . . . . 6.13 Volume attribute function + contact model . . . . . . . . . . . . . . . 6.15 Contact types in Leighton Buzzard . . . . . . . . . . . . . . . . . . . 6.16 Varying ψ for flat contacts . . . . . . . . . . . . . . . . . . . . . . . . 6.17 Varying ψ for breakage contacts . . . . . . . . . . . . . . . . . . . . . 6.18 Varying ψ for breakage contacts . . . . . . . . . . . . . . . . . . . . . 6.19 Varying ψ for breakage contacts . . . . . . . . . . . . . . . . . . . . . 6.20 Breakage capture in zeolite grains with ψ = 0.38 . . . . . . . . . . . . 7.1 More complex breakage operation through the stacking of three breakage modes. Both left and right are further breakaged . . . . . . . . . 7.2 Textured kalisphera grain . . . . . . . . . . . . . . . . . . . . . . . . 7.3 Normalized correlation coefficient as the breakage centre is varied along the Y-axis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.4 Normalized correlation coefficient as the breakage centre is for textured and non-textured grains . . . . . . . . . . . . . . . . . . . . . . 7.5 Normalized correlation coefficient for broken kalisphera grains, where the number of particles is 2, 3 and 4 respectively. . . . . . . . . . . . 7.6 Normalized correlation coefficient for kalisphera grains of different sizes100 8.1 Slices for 3D image sequence showing Leighton buzzard sand undergoing progressive breakage from time-step 1 to time-step 7 . . . . . . 8.2 3D image sequence showing Leighton buzzard sand undergoing progressive breakage from time-step 1 to time-step 7 . . . . . . . . . . . 8.3 Leighton buzzard sand grain undergoing breakage and further breakage under progressive loading . . . . . . . . . . . . . . . . . . . . . . 8.4 Detecting breakage in Leighton buzzard sand . . . . . . . . . . . . . . 8.5 Detecting breakage in zeolite grains . . . . . . . . . . . . . . . . . . . 8.6 Breakage detection from time-step 1 to time-step 2 using the back correlation algorithm for zeolite grains undergoing compression and progressive breakage . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.7 Breakage detection from time-step 2 to time-step 3 using the back correlation detection algorithm in Leighton buzzard sand grains undergoing compression and progressive breakage . . . . . . . . . . . . . 8.8 Breakage detection from time-step 2 to time-step 3 using the back correlation algorithm in zeolite grains undergoing compression and progressive breakage . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.9 Breakage detection from time-step 3 to time-step 4 using the back correlation algorithm in zeolite grains undergoing compression and progressive breakage . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.10 Breakage detection from time-step 2 to time-step 3 using the back correlation algorithm in zeolite grains undergoing compression and progressive breakage . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.11 Breakage detection from time-step 3 to time-step 4 using the back correlation detection algorithm in Leighton buzzard sand grains undergoing compression and progressive breakage . . . . . . . . . . . . . 8.12 Breakage detection from time-step 5 to time-step 6 using the back correlation detection algorithm in Leighton buzzard sand grains undergoing compression and progressive breakage . . . . . . . . . . . . . 8.13 Breakage detection from time-step 6 to time-step 7 using the back correlation detection algorithm in Leighton buzzard sand grains undergoing compression and progressive breakage . . . . . . . . . . . . . 8.14 Breakage detection from time-step 3 to time-step 4 using the back correlation detection algorithm in Leighton buzzard sand grains undergoing compression and progressive breakage . . . . . . . . . . . . . A.1 Runtime for noise filtering algorithms for image volumes . . . . . . . A.2 Runtime for noise filtering algorithms for image volumes . . . . . . . A.3 Runtime for guided filter on image volumes . . . . . . . . . . . . . . . B.1 Noise regression plots . . . . . . . . . . . . . . . . . . . . . . . . . . . B.2 Noise plots (bar) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.3 Noise plots (Graph) . . . . . . . . . . . . . . . . . . . . . . . . . . . . C.1 Noise removal effect of different filters on sapphires spheres . . . . . . C.2 Noise removal effect of different filters on zeolite . . . . . . . . . . . . C.3 Noise removal effect of different filters on highly decomposable granite C.4 Noise removal effect of different filters on hostun sand . . . . . . . . . C.5 Blur effects of total variation filtering . . . . . . . . . . . . . . . . . .

4. 1

 1 Calculated noise for dataset images . . . . . . . . . . . . . . . . . . . 4.2 Architecture details for the encoder. k refers to kernel size for the convolution, s to stride, p to padding. . . . . . . . . . . . . . . . . . . 4.3 Average noise on the natural datasets. All results use the same encoder and the model based approach. S.Spheres refers to Sapphire spheres, H. sand refers to hostun sand, C. ooids refers to Caicos ooids, L. Buzzard to Leighton buzzard, H.D Granite to Highly decomposable granite. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.1 Change in average grayscale variance computed over grains undergoing progressive breakage . . . . . . . . . . . . . . . . . . . . . . . . .

  Abbreviation Description τ Linear attenuation co-efficient CC Connected Components C A set of Contacts z Depth e A graph edge E Set of edges G Graph representation H Hierarchy of partitions I Image representation σ 2 Noise P A Partition in an Hierarchy P A path between two pixels J Projection intensity r Labeled Image region T The subgraph of a graph G Γ Total attenuation v A vertex in a graph G V Set of vertices of a graph

Figure 1 . 1 :

 11 Figure 1.1: Breakage across the temporal stack in zeolite grains undergoing oedometric compression

2. 1 . 1

 11 From x-rays to images X-ray Computed Tomography (CT) is a non-destructive technique for visualizing the internal features of solid objects, and for obtaining digital information on their 3D geometry and properties. It has been used extensively in experimental geomechanics to image geomaterials. Its use helps to understand the material composition and its fundamental behavioral responses.

Figure 2 . 3 :

 23 Figure 2.3: Relative breakage according to Einav in [77]

Figure 2 . 4 :

 24 Figure 2.4: Different breakage types[START_REF] Daouadji | An enhanced constitutive model for crushable granular materials[END_REF]

Figure 2 . 5 :Figure 2 . 6 :

 2526 Figure 2.5: Detecting breakage in zeolite images of granular materials undergoing compression[START_REF] Karatza | Study of temporal and spatial evolution of deformation and breakage of dry granular materials using x-ray computed tomography and the discrete element method[END_REF]

Figure 3 . 1 :

 31 Figure 3.1: Laboratoire 3SR x-ray scanner

Figure 3 . 2 :Figure 3 . 3 :

 3233 Figure 3.2: Dataset slices: Image slices shown from 3D images of sapphire spheres

Figure 3 . 4 :

 34 Figure 3.4: Dataset slices: Image slices shown from 3D images of Leighton buzzard

Figure 3 . 5 :Figure 3 . 6 :

 3536 Figure 3.5: Dataset slices: Image slices shown from 3D images of Highly Decomposed Granite

3. 2 . 6 Figure 3 . 7 :

 2637 Figure 3.7: Dataset slices: Image slices shown from 3D images of Hostun sand

Figure 3 . 8 :

 38 Figure 3.8: Dataset slices: Image slices shown from 3D images of Kalisphera

Figure 3 . 9 :

 39 Figure 3.9: Dataset slices: Image slices shown from 3D images of sapphire spheres, zeolite, highly decomposable granite, caicos ooids, leighton buzzard and hostun sand imaged under different mechanical conditions

  simple definition of f (I) can be as a transformation of I by the mean gray value of the entire image. This implies that the noiseless representation is estimated to be the mean value constant image I µ . In any image with distinct foreground and background, this does not work. It induces the loss of texture and latent content in the image. However its failure gives insight into potential formulations of the function f (I).

  (a) Kalisphera (b) Kalisphera(Noise=0.1) (c) Noise plot

Figure 4 . 1 :

 41 Figure 4.1: Noise is progressively added to the Kalisphera image and calculated using the computational model for noise computation

Figure 4 . 2 :

 42 Figure 4.2: Model diagram for the noise quantifying encoder. Conv. refers to convolution.

4. 4

 4 Training Procedure Algorithm 1 Training procedure 1: for b = 1, ..., B do 2:

  images. Since the training was done slice by slice, this resulted in a training set of 3500 samples. For the validation set, we used 2 500 × 500 × 500 images not in the training set, resulting in 1000 slice samples. The validation set is used during training to indicate how the model can compute noise on data it has not seen before.For the final evaluation after training, a test set of 2 500 × 500 × 500, different from the validation set is also used.Data augmentation was used to amplify the variance of these samples. The augmentations used were: random holes to simulate texture and squashing the dynamic range of the image. With a probability of 0.1, 20 holes of radius 2 pixels was randomly inserted in the considered slice. Also, the dynamic range of the slices were sometimes squashed to allow the network learn a robust noise model. The original range of the image was between 64 and 192. With a probability of 0.4, the lower limit is shifted to a random position between 64 and 100 and the upper limit is shifted to a random position between 150 and 192.

Figure 4 . 3 :

 43 Figure 4.3: Noise estimation in Kalisphera datasets using a Model based approach and a Neural network based approach

10 )ā

 10 and b represent mean filtered versions of the parameters a k and b k .

(

  grain) phase. Image thresholding is the separation of pixel values into different pixel groups or classes, whose pixel values are sometimes modeled by Gaussian distributions. Most image thresholding algorithms aim to separate the background from the foreground in an image. One of the many preferred methods for achieving such separation is Otsu thresholding. Otsu thresholding maximizes the inter-class variance while minimizing the intra-class variance of pixel values [17]. These classes mostly correspond to the different phases in the grain assembly, i.e., solid/air or solid/water/air.

Figure 6 . 1 : 7 × 7

 6177 Figure 6.1: 7 × 7 image and its partitions into α-connected components [80]

Figure 6 . 2 : 7 × 7 Algorithm 2 Figure 6 . 3 :

 6277263 Figure 6.2: 7 × 7 image and its partitions into α -ω-connected components [80]

Figure 6 . 4 :

 64 Figure 6.4: Oversegmentation of two zeolite grains using regular hierarchical watershed

Figure 6 . 5 :

 65 Figure 6.5: Given structural properties; radius and 3D position, kalisphera generates a 3D image containing 5521 grains as shown in figure 6.5a. A slice of this 3D image is shown in figure 6.5b. Due to near image uniformity, the image is corrupted with additive gaussian noise with 0.1 standard deviation as shown in figure 6.5c. Afteradding noise, the maximum pixel value was thus 255 and the minimum 0. An initial segmentation with 200, 000 segments is obtained using simple linear iterative clustering.[START_REF] Achanta | SLIC superpixels compared to state-of-the-art superpixel methods[END_REF]

Algorithm 3

 3 Contact model based hierarchical segmentation minima = GenerateMinima(image) segments = GenerateInitialMSF(minima, image) while CountMinima(minima) ≥ 1 do attributeRank = RankMinima(minima, segments, attributeFilter) contactRank = RankMinima(minima, segments, contactModel) ranks = attributeRank • ψ contactRank RemoveMinima(minima, ranks)

Figure 6 . 6 :

 66 Figure 6.6: The first three timesteps of the generated kalisphera dataset shown in a, b, and c. The corresponding images with noise added is shown in d, e, and f. The images are thresholded using otsu threshold(g, h, i) and an oversegmentation map is generated for each

Figure 6 . 7 :

 67 Figure 6.7: Slice and 3D image showing Leighton buzzard sand depicted as LB, cropped cross section of Sapphire spheres depicted as SS and cropped Zeolite.

Figure 6 . 9 :

 69 Figure 6.9: BCE-Fragmentation curves showing the performance against watershed variants. The results are averaged over 7 Kalisphera Images

Figure 6 .Figure 6 . 11 :

 6611 Figure 6.10: BCE-Fragmentation curves averaged over 7 Kalisphera images

Figure 6 . 12 :

 612 Figure 6.12: Two sapphire grain segmentation: In sapphire spheres, this shows how a contact aware segmentation can reduce or minimise the oversegmentation of granular materials.

Figure 6 . 13 :

 613 Figure 6.13: Volume attribute function + contact model: Shows the saliency values at the contact of granular materials.

Figure 6 . 14 :

 614 Figure 6.14: Power watershed compared against Contact based hierarchical watershed (Dynamics + contact).

Figure 6 . 15 :

 615 Figure 6.15: Observed contact types in Leighton buzzard sand.

68 Figure 6 . 16 :

 68616 Figure 6.16: Observed segmentation for flat contacts using different values of ψ.

(a) ψ = 0. 1 ( 56 Figure 6 . 17 :

 156617 Figure 6.17: Observed segmentation for contacts due to grain breakage using different values of ψ.

66 Figure 6 . 18 :

 66618 Figure 6.18: Observed segmentation for multiple point contact using different values of ψ.

(a) ψ = 0. 2 (Figure 6 . 19 :

 2619 Figure 6.19: Segmentation results from the application of the contact model to Leighton buzzard sand. c and d shows extracted grains undergoing undersegmentation when contact model is applied

Figure 6 . 20 :

 620 Figure 6.20: Breakage capture in zeolite grains with ψ = 0.38 compared against morphological watershed.

Figure 7 . 1 :

 71 Figure 7.1: More complex breakage operation through the stacking of three breakage modes. Both left and right are further breakaged

Figure 7 . 2 :

 72 Figure 7.2: Localized texture in artificially generated kalisphera grains

Figure 7 . 3 :Figure 7 . 4 :

 7374 Figure 7.3: Normalized correlation coefficient as the breakage centre is varied along the Y-axis

Figure 7 . 5 :

 75 Figure 7.5: Normalized correlation coefficient for broken kalisphera grains, where the number of particles is 2, 3 and 4 respectively.

Figure 7 . 6 :

 76 Figure 7.6: Normalized correlation coefficient for kalisphera grains of different sizes

8. 1

 1 Back-correlation based detectionInstead of computing a forward correlation from time-step 0 to time-step 1 and thresholding the normalized correlation coefficient to detect broken grains, we formulate a backward correlation from time-step 1 to 0. This involves performing a segmentation(labeling) of the grains in time-step 1 and correlating them backwards with the previous time-step. Although this eliminates the dependence on a tuned threshold parameter, it introduces a new dependence on the quality of segmentation.

Figure 8 . 1 : 7 (

 817 Figure 8.1: Slices for 3D image sequence showing Leighton buzzard sand undergoing progressive breakage from time-step 1 to time-step 7

Figure 8 . 2 : 7 (Figure 8 . 3 :

 82783 Figure 8.2: 3D image sequence showing Leighton buzzard sand undergoing progressive breakage from time-step 1 to time-step 7

contains 7

 7 time-series images, and progressive breakage can be observed from the second time-step. The progression of breakage can be observed in sample slices from the 3D volume at increasing time-steps shown in figure 8.1 and figure 8.2.

Figure 8 . 4 :Figure 8 . 5 :

 8485 Figure 8.4: Detecting breakage in Leighton buzzard sand

Figure 8 . 6 :

 86 Figure 8.6: Breakage detection from time-step 1 to time-step 2 using the back correlation algorithm for zeolite grains undergoing compression and progressive breakage

Figure 8 . 7 :

 87 Figure 8.7: Breakage detection from time-step 2 to time-step 3 using the back correlation detection algorithm in Leighton buzzard sand grains undergoing compression and progressive breakage

Figure 8 . 8 :

 88 Figure 8.8: Breakage detection from time-step 2 to time-step 3 using the back correlation algorithm in zeolite grains undergoing compression and progressive breakage

Figure 8 . 9 :

 89 Figure 8.9: Breakage detection from time-step 3 to time-step 4 using the back correlation algorithm in zeolite grains undergoing compression and progressive breakage

Figure 8 . 10 :

 810 Figure 8.10: Breakage detection from time-step 2 to time-step 3 using the back correlation algorithm in zeolite grains undergoing compression and progressive breakage

Figure 8 . 11 :

 811 Figure 8.11: Breakage detection from time-step 3 to time-step 4 using the back correlation detection algorithm in Leighton buzzard sand grains undergoing compression and progressive breakage

Figure 8 . 12 :

 812 Figure 8.12: Breakage detection from time-step 5 to time-step 6 using the back correlation detection algorithm in Leighton buzzard sand grains undergoing compression and progressive breakage

Figure 8 . 13 :

 813 Figure 8.13: Breakage detection from time-step 6 to time-step 7 using the back correlation detection algorithm in Leighton buzzard sand grains undergoing compression and progressive breakage

(a) t = 3 (b) t = 4 Figure 8 . 14 :

 34814 Figure 8.14: Breakage detection from time-step 3 to time-step 4 using the back correlation detection algorithm in Leighton buzzard sand grains undergoing compression and progressive breakage

Figure A. 1 :

 1 Figure A.1: Runtime for noise filtering algorithms for image volumes: The image is a cube volume with defined width. The cube is progressively increased and the runtime is measured

Figure A. 2 :

 2 Figure A.2: Runtime for noise filtering algorithms for image volumes: The image is a cube volume with defined width. The cube is progressively increased and the runtime is measured. Plotting the volume against the measured time

Figure B. 1 :

 1 Figure B.1: Noise regression plots: Noise is progressively to Sapphire spheres, Hostun sand, Caicos ooids, Highly decomposable granite, Leighton buzzzard and Zeolite.The internal noise is the calculated using the noise computation model. The plot is analysed for how well it fits a linear model.

Figure B. 2 :

 2 Figure B.2: Noise regression plots: Noise is progressively to Sapphire spheres, Hostun sand, Caicos ooids, Highly decomposable granite, Leighton buzzard and Zeolite and filtered using Guided filtering, Bilateral filtering, Anisotropic diffusion and total variation filtering. The internal noise is the calculated using the noise computation model. The bar plot representation gives inference into how much noise remains after filtering with a particular filter.

Figure B. 3 :

 3 Figure B.3: Noise regression plots: Noise is progressively to Sapphire spheres, Hostun sand, Caicos ooids, Highly decomposable granite, Leighton buzzard and Zeolite and filtered using Guided filtering, Bilateral filtering, Anisotropic diffusion and total variation filtering. The internal noise is the calculated using the noise computation model. The graph gives inference into how much noise remains after filtering with a particular filter.

Figure C. 2 :

 2 Figure C.2: Different filters are applied to zeolite image: Anisotropic diffusion, bilateral filter, total variation filter, guided filter

Figure C. 3 :

 3 Figure C.3: Different filters are applied to highly decomposable granite image: Anisotropic diffusion, bilateral filter, total variation filter, guided filter

Figure C. 4 :

 4 Figure C.4: Different filters are applied to hostun image: Anisotropic diffusion, bilateral filter, total variation filter, guided filter

Figure C. 5 :

 5 Figure C.5: Blur effects of total variation filtering: The weight parameter in total variation filtering is varied. An increase in the value implies less data fidelity. This is shown to increase the blur in the sapphire spheres image

Table 4 .

 4

2: Architecture details for the encoder. k refers to kernel size for the convolution, s to stride, p to padding.

table 4 .

 4 [START_REF] Rosenblatt | A central limit theorem and a strong mixing condition[END_REF]. One way to curb overfitting would be to introduce artificial grains with more angular shapes. In reality, this might be difficult to achieve. As in DEM, artificial grains with more angularity might be generated by generating clusters of spheres. These clustered spheres should yield complex shapes that aid in overcoming overfitting. However, it would be difficult to generate realistic mechanical configurations and would require significant compute resources.

			Estimated noise
	Method	S.Spheres Zeolite H. sand C. ooids L. Buzzard H.D. Granite
	Model(5000bins)	0.157	0.1625 0.145	0.12	0.043	0.124
	Encoder	0.198	0.174	0.012	0.152	0.01	0.03

The model being trained on spherical kalisphera grains, overfits. It implicitly learns that the grain for which the noise distribution is needed is spherical. As such, it gives realistic noise values for spherical or near-spherical grains. This is can be observed in

Table 4 .

 4 3: Average noise on the natural datasets. All results use the same encoder and the model based approach. S.Spheres refers to Sapphire spheres, H. sand refers to hostun sand, C. ooids refers to Caicos ooids, L. Buzzard to Leighton buzzard,

H.D Granite to Highly decomposable granite.

Table 8 .

 8 1: Change in average grayscale variance computed over grains undergoing progressive breakage is gradually removed. Furthermore, any method to detection of crushed grains hinged on texture based methods such as correlation may become less effective as the granular state approaches the ultimate distribution. Although we did not fully investigate this phenomenon, we show by a simple experiment the gradual reduction in texture, which we approximate as pixel variance. To verify this texture dissipation, we track the average change in variance of 4 grains having similar initial gray value variance, where we observe a decrease in variance as time-step increases. This is shown in table 8.1. Analyzing the next time-step in the Leighton buzzard series (see figure8.11), we observe less false positives. We argue that this is due to its size and inherently more diverse texture component.

Acknowledgements

I would like to start by thanking AGIR, under whose auspices and finance this project was made possible. Next, I thank my supervisors at the GIPSA laboratory Mauro Dalla-Mura and Laurent Bonnaud for their help, time and direction.

6.5 Result and Evaluation

Evaluation of hierarchies on kalisphera

In evaluating the effect of adding a contact model, we determine the optimal ψ value for the attribute combination function. The ψ parameter as iterated is the weighting on the combination function. This weighting determines the contribution of the contact model to the ranking of a minimum. To determine the optimal weighting, we experiment with values from 0 to 1. We observe that optimal values for the combination of dynamics and contact model range between 0.3 and 0.4. But, optimal values for the combination of volume and our model range between 0.4 and 0.5. We show this in figure 6.8a. We should state that ψ values and BCE do not exhibit a linear relationship. This implies that increasing the value of ψ does not guarantee a better segmentation.

Having obtained optimal ψ values that can be used, we analyze the effect of using the contact model. Figure 6.8 and 6.10 shows that the introduction of the contact model allows merging to reach a more optimal segmentation. We characterize optimal hierarchical segmentation as valid image structures persisting in the hierarchy. This implies that the addition of our contact model leads towards the was learned;

• Wider breakage gaps require low threshold values

• When small grains break, they experience a more significant decrease the correlation coefficient than when large grains break.

• Unevenly distributed texture on grains result in a different correlation coefficient depending on which of the particles has the most texture component.

• In a simple splitting breakage event, even sized particles result in higher correlation coefficient, provided that they have uniform texture

In chapter 8 we derive a new method for breakage detection and tracking. In this derived method we map a grain particle (broken) to a previous state. When two labels map to a single grain in the previous time-step, they are characterized as broken. Such backward characterization can be achieved using a minima set or the entire label. We show how our algorithm can be used to rectify oversegmentation in granular materials. Although the algorithm was outlined and shown to work in Leighton Buzzard sand and Zeolite, extensive studies was not made on a wide range of datasets. This is a natural extension to the proposed method.

Appendix A

Algorithmic speed and efficiency of noise removal algorithms

In evaluating the performance of these algorithms, the speeds of guided filtering, bilateral and total variation denoising are evaluated. Anisotropic diffusion is not evaluated due to the fact, that it is iterative and an evaluation of its runtime is inherently very subjective. Classical bilateral filtering is slow, although speed up schemes have been proposed in [START_REF] Paris | Bilateral filtering: Theory and applications[END_REF].

The experiment was run using the Scipy implementation of Total variation denoising, the ITK implementation of Bilateral filtering and an implementation of guided filter. The guided filter was implemented in C++ but coupled to python using Pybind11. The seemed a fair fight as ITK and scipy are C++ developed but have python bindings.

The runtime scales linearly with increase in image volume. This is shown in

Noise response in natural datasets

It is necessary to make a distinction between edge-based methods and model based methods. Although all of the methods are edge preserving, edge-based methods are modelled on the preservation of edges while model based methods generate an output that is described as mathematical transformation of the input and in effect, edge preserving. Total variation models noise as an additive Gaussian while guided filtering models the transformation as a linear transformation.

The addition of two Gaussian, with width σ 1 and σ 2 will yield a Gaussian with σ 1 + σ 2 [START_REF] Valentin V Petrov | Limit theorems of probability theory: sequences of independent random variables[END_REF]. This follows the principle of additivity outlined by linearity.

Thus, if a noise component of an image is an additive Gaussian, it should respond linearly to the addition of noise. In the cases where the noise response can be closely approximated as linear, Highly decomposable granite 3.5 and Leighton buzzard 3.4, Model based algorithms perform better (Guided filtering and Total variation). The linear approximations of its noise response is seen B.1d, B.1e. However even in these cases where the image response to noise is closely depicting that the initial noise value is an additive Gaussian, the variational approach out performs guided filtering at high noise levels. With increase in noise, parameter in guided fails to distinguish between noise and edges. An example is in figure B.3 where guided filtering performs well at low noise but its performance quickly degrades in high noise. It should be noted that bilateral filtering does almost as well in these scenarios, but apart from the fact that it is outperformed, its computational cost makes it insufficient.

In case of non linear response of noise, signifying that the underlying noise is not an additive Gaussian, Bilateral and anisotropic diffusion are consistent in their performance. However, in cases where grain size per pixel is small and images have large sized grains, guided filtering performs better across all noise realisations. This is seen in sapphire spheres B.1a and highly decomposable granite B.1d. It can be said that although noise is high, there is enough statistical information for the filter.

Appendix C

Denoising of dataset images