En 2012, un cancer a été diagnostiqué pour 14 millions de personnes environ, alors qu'approximativement 8,2 millions de personnes meurent chaque année de cette maladie [Stewart and Wild, 2014]. La plus grande cause de mortalité liée au cancer est due à les lésions métastatiques [START_REF] Adams | Capture and enumeration of circulating tumor cells from peripheral blood using microfluids[END_REF]. L'étude de Cristofanilli et al. [START_REF] Cristofanilli | Circulating Tumor Cells, Disease Progression, and Survival in Metastatic Breast Cancer[END_REF] montre le lien entre les cellules tumorales circulantes et les taux de survie de patients atteints de cancers. Seulement 20% de ces patients survivent plus de 18 mois lorsqu'il y a plus que 5 cellules tumorales circulantes pour 7.9 mL de sang. Le taux de survie grimpe jusqu'à 75% pour moins que 5 cellules tumorales circulantes par 7.9 mL de sang. La séparation des cellules tumorales et des fluides corporels est donc une procédure importante dans le cadre du diagnostic clinique et de la thérapie du cancer.

Les développements modernes dans le domaine de la séparation de cellules sont basés sur la chromatographie cellulaire dans laquelle l'interaction spécifique des cellules cibles avec la surface interne des médias est utilisée pour la séparation [START_REF] Adams | Capture and enumeration of circulating tumor cells from peripheral blood using microfluids[END_REF], Dharmasiri et al., 2009, Dharmasiri et al., 2010, Plieva et al., 2011]. En général, les processus chromatographiques sont utilisés pour séparer des éléments composés comme les cellules, les organelles, les gouttelettes ou les particules. Dans le processus de filtrage, une suspension avec une population de cellules est pompée à travers l'espace poreux d'un matériau et un ou plusieurs types de cellules d'une population sont retenus par la surface interne du milieu. L'efficacité de ce processus dépend à la fois de la géométrie du milieu et du revêtement de la surface interne constitué d'une substance qui, activée de manière appropriée, autorise la rétension sélective de cellules, ainsi illustré par la Fig. 1. L'objectif est d'améliorer la séparation des cellules, un aspect important de ce travail résidant dans la conception de filtres dont la capacité de filtrage est suffisamment élevée. Cette thèse est un premier pas vers la conception de filtres optimaux. Dans le but d'optimiser le filtrage chromatographique en termes de taux de déposition de cellules cibles, il est nécessaire de déterminer les caractéristiques géométriques adéquates de l'espace poreux. En particulier, la torsion des trajectoires des particules à travers l'espace poreux a un impact significatif sur le processus chromatographique. v Figure 1: Filtrage chromatographique de cellules rouges et vertes. La taille des pores du milieu est plusieurs fois plus grande que la taille des cellules. La surface est revêtue avec des molécules liantes (représentées par des cheveux violet) qui peuvent capturer les cellules vertes. Les cellules en vert foncé sont actuellement attachées à la surface. Les cellules vertes sont plus durablement immobilisées que les cellules rouges, cet effet étant à l'origine de la séparation.

Défis scientifiques

Les candidats les plus prometteurs en tant que filtres chromatographiques sont les mousses ouvertes (macroscopiquement homogènes) et les systèmes aléatoires de fibres (voir Fig. 2). Pour cette raison, les structures de ces milieux sont étudiées avec soin au regard de leur porosité, de leur surface spécifique, de leur tortuosité, etc.

Mis-à-part ces caractéristiques, la torsion des trajectoires des particules à travers l'espace poreux est une quantité importante pour la caractérisaton des interactions entre cellules et surfaces. Une torsion élevée de l'espace poreux induit une force centrifuge qui déplace les cellules vers les surfaces internes de sorte que les interactions entre cellules et surface augmentent (en d'autres termes, les chemins de particules avec une torsion élevée agissent comme des micro-centrifugeuse).

Il existe diverses approches d'estimation de la torsion locale à partir d'images 3D. Celles-ci diffèrent par leur stratégie de prise en compte des effets de discrérisation. Les approches suggérées dans [START_REF] Kim | Frenet-Serret and the Estimation of Curvature and Torsion[END_REF], [Kehtarnavaz and de Figueiredo, 1988] et [START_REF] Sangalli | Efficient estimation of three-dimensional curves and their derivatives by free-knot regression splines, applied to the analysis of inner carotid artery centrelines[END_REF] sont basées sur des splines de lissage ajustées aux courbes discrètes, où la courbure et la torsion sont estimées à partir des splines à l'aide des formules de géométrie différentielle bien connues. Une approche similaire est utilisée dans [START_REF] Lewinger | Curvature and torsion estimators based on parametric curve fitting[END_REF], Worz et al., 2010], des modèles de courbes continues étant ajustées aux courbes discrètes à l'aide de méthodes des moindres carrés. Cette représentation dans l'espace continu est à nouveau utilisée pour l'évaluation des torsions des courbes. Dans vi 2. Défis scientifiques (a) (b)

Figure 2: Rendu de deux jeux de données tomographiques visualisant des sous-volumes de 512 3 pixels extraits d'un volume de 2048 3 pixels. La longueur des côtés des pixels cubiques est de 1.1 µm. (a) une mousse ouverte et (b) un système de fibres.

[ [START_REF] Pao | Bending and twisting of an in vivo coronary artery at a bifurcation[END_REF] et [START_REF] Crenshaw | Analysis of the three-dimensional trajectories of organisms: estimates of the velocity, curvature and torsion from positional information[END_REF], des splines bi-cubiques sont ajustées respectivement aux positions de repères anatomiques et à celles d'organismes. [Mokhtarian, 1997] utilise une description multi-échelle de la forme des courbes en appliquant à celles-ci un lissage gaussien anisotropique pour atténuer les effets de la discrétisation sur l'estimation des torsions. Il est possible d'accélérer ce lissage dans l'espace de Fourier, ce qui conduit à un calcul plus rapide des torsions [START_REF] Medina | Curvature and torsion estimation for coronary-artery motion analysis[END_REF]. Cette approche est particulièrement intéressante puisque les dérivées requises pour le calcul de la torsion de la courbe 3D peuvent également être calculées dans l'espace de Fourier (simultannément avec le lissage). Le coeur de l'approche algorithmique suggérée par [START_REF] Nguyen | Curvature and Torsion Estimators for 3D Curves[END_REF] et [START_REF] Nguyen | On the Local Properties of Digital Curves[END_REF] réside dans une décomposition des courbes discrètes en des segments dits "bruités", ce qui est fondamentalement équivalement à un lissage.

Pour résumer, excepté la dernière méthode basée sur les "segments bruités", la tendance générale dans la litérature est d'approximer et de lisser les courbes avec des fonctions connues analytiquement et de déterminer leur courbure et torsion à l'aide des formules de géométrie différentielle.

L'ajustement et le lissage dépend usuellement de divers paramètres dont les valeurs doivent souvent être subjectivement fixées. Plus précisément, l'ajustement par des courbes continues dans le but d'atténuer les effets de discrétisation peut uniquement conduire à une estimation appropriée de la torsion lorsque des connaissances a priori (par exemple le type et l'ordre des fonctions) sont disponibles. Le calcul de la torsion est basé sur l'évaluation de dérivées d'ordre élevé fortement influencées par le bruit de discrétisation et étant à l'origine de grandes erreurs d'estimations. Des connaissances a priori exactes sont indisponibles pour la majorité des applications, et en particulier dans le cas de la vii séparation des cellules. Par ailleurs, les valeurs appropriées des paramètres de lissage dépendent de la résolution latérale (c'est-à-dire du taux d'échantillonnage), aussi bien que des changements locaux de la forme des courbes. Enfin, certaines approches sont uniquement adaptées à des courbes formées de points connectés, alors que les trajectoires des cellules suivies dans le temps sont données sous la forme de points discrets la plus part du temps non connectés à cause de la faible résolution temporelle.

Dans le cadre de cette thèse, nous proposons de nouvelles approches pour l'estimation de la courbure et de la torsion de chemins 3D à travers des milieux poreux. Nous avons développé deux approches d'estimation des courbures et des torsions, chacune de ces approches étant dédiée à un type spécifique de données.

La première approche a été developpée pour des données tomographiques 4D, c'està-dire pour des séquences d'images 3D (comme celle montrée dans la Fig. 2) résolue dans le temps. Une suspension avec des cellules est pompée à travers le milieu poreux durant l'acquisition des données. Chaque trajectoire discrète des particules est obtenue en suivant une cellule dans la séquence des images 3D. Ce conduit, pour chaque chemin, à une liste (succession) de points déconnectés avec une distance variable entre les positions. Cet échantillonnage à pas variable des chemins 3D est dû, d'une part, aux trajectoires et vitesses changeantes des particules aux cours des acquisitions et, d'autre parts, à cause de la faible résolution temporelle comparée à la vitesse des particules. Le défi, pour ce type de données, est de trouver une méthode d'estimation de la torsion ne reposant sur aucune hypothèse de modèle de courbe continues comme celles utilisées systématiquement dans la littérature.

Les données 4D décrites dans le paragraphe précédent sont acquises avec un dispositif expérimental sophistiqué proposé dans ce travail. Cette acquisition de données 4D est un challenge, la quantité des données à stocker étant énorme et les temps de traitement très importants. Même si le rôle de ces données 4D est très important dans ce travail, nous avons aussi utilisé des images 3D des milieux poreux sans suspension qui les traversent. Ainsi décrit dans la thèse, des chemins ont été simulés à travers ces volumes. Contrairement aux données 4D, les chemins dans les volumes 3D sont représentés par des séquencs de pixels connectés. Alors que pour les données 4D l'objectif est de trouver une approche purement discrète, le défi pour les données 3D est d'obtenir une solution précise basée sur des modèles de courbes continues sans connaissance a priori quant à leur type et degré, et avec peu de paramètres facilement ajustables.

En résumé, le défi de cette thèse est de trouver deux nouvelles méthodes d'estimation de torsion pour des chemins de pixels connectés et et non connectés. Ces méthodes doivent atténuer l'effet du bruit de discrétisation et du bruit (erreurs) de localisation des particules et permettre un sélection aisée de paramètres. viii 3 Principales contributions Cette section décrit brièvement les quatre principales contributions de la thèse. Le dispositif expérimental d'acquisition des donnéée basé sur la micro-tomographie (➭CT) est d'abord introduit. Nous décrivons ensuite la détermination des chemins dans les données 3D et 4D. Dans les deux dernières sections, la méthode "purement discrète" d'estimation des torsions dans les données 4D et la méthode d'estimation des torsions dans les données statiques 3D sont présentées.

Dispositif expérimental

Un dispositif expérimental a été développé pour obtenir des images3D résolues dans le temps avec le ➭CT de l'Installation Européenne de Rayonnement Synchrotron de Grenoble. Le dispositif expérimental pour l'acquisition des données 4D consiste en un milieu poreux placé dans une cartouche cylindrique qui est connecté à une pompe seringue, ainsi illustré dans la Fig. 3. Durant ces expériences in situ, une suspension avec des particules est pompée dans la colonne avec la cartouche cylindrique et une séquence d'images 3D est acquise. La suspension inclue des billes sphériques en silicate avec un diamètre de 20 ➭m et suspendues dans une solution d'alcool (les particules en silicate agissent en tant que substitut des cellules mammaliennes vivantes). Pour réaliser ces acquisition de ➭CT 4D, la densité du flux des photons, l'augmentation du contraste à l'aide d'une technique de contraste de phase et le matériel utilisé doivent être coordonnés avec soin.

(1)

(2)

(3) (4) Cet aspect du travail a été discuté dans [START_REF] Blankenburg | Torsion estimation of particle paths through porous media observed by in-situ timeresolved microtomography[END_REF].

Extraction des chemins des particules

Dans cette thèse, selon la nature des données (données 3D ou 4D), deux types de chemins de particules sont considérés. Le premier type de trajectoires 3D extraites des données 4D est directement obtenu via l'observation de particules lentes à différents pas temporels (voir Fig. 4). Un suivi complètement automatique des particules est très difficile et hasardeux à cause de la grande taille des images, de la présence d'artéfacts dans les images et de la faible résolution temporelle. Pour cette raison, nous proposons dans ce travail une méthode semi-automatique de suivi de particules.

Le deuxième type de trajectoires, extraites des données 3D, est obtenu avec une version modifiée de la "Méthode de Fast Marching" (MFM) introduite par [START_REF] Malladi | Level set and fast marching methods in image processing and computer vision[END_REF].La MFM calcule la solution discrète de léquation iconale qui est à la base de la détermination du chemin le plus rapide à travers l'espace poreux. La transformée de distances est donc appliquée au milieu poreux, les distances obtenues pouvant être inteterpées en tant que vitesses des particules. Des vitesses élevées sont attribuées à des positions éloignées des surfaces internes. Ainsi illustré dans , l'algorithme est utilisé pour obtenir différents chemins caractérisants l'espace vide du milieu poreux.

Cet aspect de notre travail a été publié dans [Blankenburg et al., 2016a].

Estimation de la torsion discrète à l'aide des quotients de différences

Dans la suite nous présentons une méthode d'estimation des courbure et torsions locales des chemins de particules. Cette méthode est basée sur la discrétisation des formules de géométrie différentielle, où les quotients différentiels apparaissant dans ces équations sont approximés par les quotients de différences correspondants. Cette approche est mise en oeuvre sans devoir recourir au réglage d'une de taille de filtre ou d'un pas d'échantillonnage. Par ailleurs, même si cette méthode a été développée pour un type spécifique de discrétisation (les position des particules mobiles sont échantillonnées sur un ensemble de position discrètes temporelles), elle peut également être appliquée à des chemins de particules qui ont été échantillonnés sur une grille (données 3D).

Le chemin d'une particule peut être représenté par une courbe C = {f (t) : t 0 ≤ t ≤ t 1 } dans l'espace et donnée sous forme de fonction paramétrique rectifiable f : R → R 3 .

x Figure 5: Chemins 3D les plus courts extraits à l'aide de l'adaptation de la "Méthode de Fast Marching" (MFM). Les chemins obtenus sont en rouge, le système poreux des filtres sont en brun et l'espace poreux est en noir.

Soient les dérivées continues ḟ , f et ... f avec ḟ > 0 et ḟ × f > 0, alors la courbure locale κ et la torsion locale τ de f en t existent et on obtient :

κ = ḟ × f ḟ 3 , τ = ( ḟ × f ) • ... f ḟ × f 2
(1) [Spivak, 1979], où × dénote le produit vectoriel.

Ainsi détaillé dans cette thèse, en utilisant les différences ∆f k (t) = f (t + k∆t)f (t + (k -1)∆t) pour k = 1, 2, 3, un pas ∆t > 0 et la moyenne géométrique pour une symmétrisation appropriée, la courbure et la torsion peuvent être ré-écrites respectivement comme suit :

κ = ∆f 1 × ∆f 2 • ∆f 2 × ∆f 3 ∆f 1 • ∆f 2 • ∆f 3 , ( 2
) τ = (∆f 1 × ∆f 2 ) • ∆f 3 ∆f 1 × ∆f 2 • ∆f 2 × ∆f 3 . ( 3
)
xi 3. Principales contributions À partir des expressions situées à droite des égalités dans (2) et (3) on obtient des estimateurs appropriés κ and τ pour respectivement κ and τ aux points f (t + ∆t) de la courbe C.

Cet aspect du travail a été publié dans [START_REF] Blankenburg | Torsion of particle trajectories through pore space and its estimation using information on local pixel configurations[END_REF][START_REF] Blankenburg | Parameter free torsion estimation of curves in 3D images[END_REF].

Estimation de la torsion basée sur l'approximation de Fourier

La méthode d'approximation de Fourier proposée pour l'estimation des courbure et des torsion démarre avec l'ajustement d'une séerie de Fourier aux points d'une courbe discrète C . Pour ce faire, il est supposé que la fonction paramétrique f : R → R 3 a une période a, c'est-à-dire qu'elle correspond à une courbe cyclique C. Si f intégrable au sens de Lebesgue sur (0, a), ces coefficients de Fourier est la m-ième somme partielle de la série de Fourier.

Il est à remarquer que les dérivées ṡm , sm , ... s m et les fonctions ṡm 2 , ṡm × sm , ṡm × sm 2 et ( ṡm × sm ) • ... s m sont aussi des sommes partielles de la série de Fourier. Pour cette raison, la somme partielle est ajustée aux points discrets d'une courbe dans l'espace pour ensuite estimer la courbure et la torsion de la courbe à partir de cette somme partielle.

Soient C les échantillons (points discrets) de la courbe C d'une fonction paramétrique f (t) qui est supposée être périodique sur [0, a). Les échantillons C correspondent à n pixels avec des positions (coordonnées) f 0 , . . . , f n-1 . La m-ième somme partielle s m (t) d'une série de Fourier peut être ajustée aux échantillons C en utilisant l'approximation c ℓ ≈ fℓ , où les f0 , . . . , fn-1 font référence à la transformée de Fourier du vecteur (f 0 , . . . , f n-1 ) vérifiant

fℓ = 1 n n-1 k=0
f k e -2πikℓ n , ℓ = 0, . . . , n -1.

En exploitant une représentation explicite par les moyennes de Cesàro correspondantes, les sommes partielles s m (t) de la courbe discrète C sont maintenant données par : Dans Eq. ( 4), le lissage et le calcul des dérivées sont réalisés simultanément en une seule étape. L'équation (1) est ensuite utilisée pour déterminer la courbure et la torsion.

Cet aspect du travail a été publié dans [Blankenburg et al., 2016a[START_REF] Blankenburg | Parameter free torsion estimation of curves in 3D images[END_REF].

Comparaison des méthodes et analyse des résultats

Les deux méthodes proposées dans cette thèse ont été testées sur diverses données de référence et comparées à la méthode d'approximation par les splines qui fait partie de l'état de l'art. La figure 6 donne un exemple de résultats pour la fonction

f (t) = r    cos t sin t c sin 2 t    , t ∈ [0, 4π), r > 0, (5)
avec r = 32 and c = 1. Cette courbe est discrétisées sur une grille de la manière suivante. D'abord la courbe est soumise à des rotations distribuées uniformément dans l'espace. Puis, des offsets sur f (et répartis de manière uniforme sur les mailles de la grille) servent de points départ sur les courbes. Finalement, les courbes sont échantillonnées sur la grille selon la discrétisation extérieure de Jordan. La figure 6 donne la valeur moyenne des torsions et l'intervalle de confiance des torsions pour chaque position t. L'approximation de Fourier utilise un paramètre de lissage valant m = 4 et la méthode des splines utilise le réglage p = 4 • 10 -4 . Il est rappelé que la méthode basée sur les quotients de différences n'implique aucun réglage de paramètre. Il est visible sur Fig. 6 que l'approche par l'approximation de Fourier a le plus petit intervalle de confiance et conduit donc à la plus petite erreur d'estimation. Toutefois, le graphe des valeurs moyennes de la méthode des quotients de différences conduit au meilleur ajustement avec la vraie torsion, tout en ayant un intervalle de confiance plus large que l'approche d'approximation de Fourier. L'approche par la méthode des splines conduit à une intervalle de confiance légèrement plus large que celui des quotients des différences et présente une valeur moyenne de torsions s'écartant le plus de la vérité terrain. Les résultats avec la méthode des splines peuvent être améliorés en adaptant le paramétre de lissage au type de la courbe, tout en conduisant à une augmentation des erreurs d'estimation des courbures et torsions lorsque la forme des courbes (qui dépend de r et c) change. En pratique, la forme des courbes est inconnue et il est donc impossible d'adapter la méthode des splines à un jeu de données avec une grande variété de courbes. Au contraire, l'approche des quotients des différences peut être utilisée sans adaptation de paramètres et offre des résultats fiables à travers l'ensemble des valeurs de l'espace des paramètres r et c. Ainsi montré dans cette thèse, à cause de son indépendence vis-à-vis des paramètres, l'approche des quotients des différences est la mieux appropriée pour l'estimation des torsions positions discrètes non connectées de particules extraites des séquences vidéos d'images 3D. Ainsi xiii également montré dans ce travail, l'approche par approximation de Fourier, est la plus adaptée pour ces pixels discret connectés obtenu à partir des espace poreux des volumes statiques. Globalement, les deux méthodes proposées ont un meilleur intervalle de confiance tout en étant moins paramètre dépendant que la méthode des splines.

Figure 6: Comparaison des résultats d'estimation de la torsion basées sur (a) l'approche populaire des splines, (b) la nouvelle approche des quotients de différences et (c) l'approche par l'approximation de Fourier avec (d) les vraies valeurs de torsion. Les lignes correspondent aux valeurs moyennes des torsions, alors que les surfaces colorées représentent les intervalles de confiance comprenant 90% des valeurs de torsions estimées.

Figure 7: Distribution des torsions des chemins des particules à travers l'espace poreux d'une mousse ouverte (en vert) et d'un système de fibres (en rouge). Les valeurs des torsions sont calculées avec la méthode d'approximation de Fourier.

La figure 7 montre la distribution des torsions des chemins de particules à travers l'espace poreux d'une mousse et d'un système de fibres. Ces histogrammes normalisés des torsions ont été calculés par l'approche d'approximation de Fourier appliquée aux chemins les plus courts extraits des données 3D. Il est à noter que la variation des torsions des trajectoires à travers le système de fibres est nettement plus grande que celle des mousses. Cette différence est due au fait que les systèmes de fibres ont une porosité beaucoup plus élevée que les mousses.

Dans la Fig. 7, il est possible d'observer un décalage significatif de l'histogramme vers les torsions négatives, tandis que l'histogramme des torsions du système de fibre est xiv

ッ N セ ウ @ , lO 111:\:o - ., • O.IS セ H [ @ " • • li セ @ :.: ' -<11 ' •o.:s • , • ' • , ' 1 •(l.IM N i I c セ @ " Q Q L A Q セ @ (l Z N セ @ .. b•,.-• 5.
Liste des publications personnelles presque symétrique et centré sur zéro. Cet effet provient du procécédé de production du système de fibre par soufflage à l'état fondu, processus dans lequel le réseau des filières des fibres est en rotation. Même si les fibres sont fortement remuées durant leur dépôt, Le système des fibres hérite de la direction de rotation du réseau de la filière de fibres.

Les conséquences de cet héritage sont double. 1) La torsion est une quantité signée et la distribution des torsion d'un chemin typique à travers l'espace des pores est une quantité caractérisant l'orientation de l'espace des pores, c'est-à-dire cela indique que le flot de haut en bas à travers l'espace des pores est orienté soit vers la droite, soit vers la gauche.

2) La distribution des torsions des chemins à travers l'espace d'un milieu poreux peut être influencé par son processus de production. À notre connaissance, notre méthode d'estimation des torsions de chemins de pixels connectés est la seule qui mesure cet effet.

Cet aspect du travail a été publié dans [Blankenburg et al., 2016a[START_REF] Blankenburg | Parameter free torsion estimation of curves in 3D images[END_REF].

Liste des publications personnelles

Journaux internationaux avec comité de lecture Blankenburg, C., Daul, C., and Ohser, J. (2016a). Estimating torsion of digital curves using 3D image analysis. Image Analysis and Stereology, 35(2):81-91. [START_REF] Blankenburg | Torsion estimation of particle paths through porous media observed by in-situ timeresolved microtomography[END_REF]. Torsion estimation of particle paths through porous media observed by in-situ time-resolved microtomography. Journal of Microscopy, 266(2):141-152.

Conférences internationales avec comité de lecture [START_REF] Blankenburg | Torsion of particle trajectories through pore space and its estimation using information on local pixel configurations[END_REF]. Torsion of particle trajectories through pore space and its estimation using information on local pixel configurations. In Pirard, E., editor, Proceedings of 14th International Congress for Stereology and Image Analysis (ICSIA), page WA04, Liège. ISS. [START_REF] Blankenburg | Parameter free torsion estimation of curves in 3D images[END_REF]. Parameter free torsion estimation of curves in 3D images. In Proceedings 23rd IEEE International Conference on Image Processing (ICIP), pages 1081-1085, Phoenix, Arizona USA.

xv

Introduction

This thesis was written at Université de Lorraine in the department of Santé-Biologie-Signal from the CRAN laboratory (Centre de Recherche en Automatique de Nancy, UMR 7039 UL/CNRS). This work resulted from the close collaboration with the Hochschule Darmstadt University of Applied Sciences (Schöfferstraße 3, 64295 Darmstadt, Germany) and was funded by the German Federal Ministry of Education and Research (BMBF) under grant MNT/05M13RCA/AniS. In this AniS project, the companies Merck KGaA, BASF SE, MANN+HUMMEL GmbH, HARTING KGaA are cooperating with the universities Darmstadt, Heidelberg, Kaiserslautern, Ulm and the Fraunhofer Institute for Industrial Mathematics. The topic of the collaborative project was the analysis of low dimensional structures in three-dimensional (3D) images. One task of this project was to investigate porous media with percolating pore space regarding their ability to separate cells out of a cell suspension. The results of these investigations are part of the presented thesis. The aim of this work is to improve the cell separation and, thus, a crucial point is to design filters so that the corresponding filtering efficiency is sufficiently high. In order to observe the microstructure of the porous media as well as the fluid-solid flow through the pore space, experiments at the European Synchrotron Radiation Facility in Grenoble were designed.

The separation of cancer cells from body fluids is an important procedure applied in clinical diagnostics and therapy for cancer detection and treatment [START_REF] Adams | Capture and enumeration of circulating tumor cells from peripheral blood using microfluids[END_REF]. For this purpose, different methods like fluorescence activated cell sorting, magnetic bead assisted cell sorting and micro-fluid separation can be applied. These laboratory methods are only feasible for small amounts of cell suspensions. Therefore, the most common separation technique is centrifugation, but centrifugation has some disadvantages regarding cell separation. On the one hand, the investment costs and maintenance time can be very high, because the uncontrolled proliferation of cells and the contamination with extraneous cells must be avoided. On the other hand, the cells can be damaged by the centrifugal process, which can cause unintentional mutations. For these reasons, alternative approaches for cell separation like cell chromatography have to be investigated.

In the classical cell chromatography, a suspension, the so-called mobile phase, is flown over a thin layer the so-called stationary phase. During the process, the different cells are temporally bounded for different periods of time on the stationary phase like shown in Fig. 1. This time difference between immobilized and non-immobilized cells Figure 1: One-dimensional chromatographic process: the ligands on the stationary phase (purple hair) are interacting with the receptors of the red cells, which causes a longer immobilization of the red cells. The resulting propagation time difference between immobilized and non-immobilized cells causes the separation. causes the separation. Modern cell chromatographs use 3D porous media as stationary phase instead of thin layers, which lead to new challenges. For such filters, it is no longer possible to investigate the suspension flow, and in particular, the cell motion with simple microscopy. Moreover, the geometric characteristics of porous media have a considerable influence on the chromatographic process. One possibility to examine 3D particle motion in open porous media is micro-computed tomography (➭CT) and in particular timeresolved ➭CT. However, the X-ray radiation involved by ➭CT, the flow conditions, the choice of the suspension and stationary phase of the chromatographic filter need all to be carefully adapted for the investigation of chromatographic processes. Moreover, the image reconstruction method, the acquisition rate, and the lateral resolution have to be chosen accordingly to the used filter material, flow rate and the particle size, respectively. Therefore, one aim of this thesis is to design an experimental setup to observe 3D motions of at least slow particles passing through porous media. Despite recent progress in the field of 3D image acquisition, it cannot be expected that the temporal resolution of the time-resolved ➭CT is sufficient to observe fast particles. An alternative is the extraction of the particle paths from a simulation of a solid/liquid two-phase flow. Such simulations are usually based on the numerical solution of discretized systems of partial differential equations, which results in very high computational effort especially when applying them to large 3D images. Therefore, another aim of this thesis is to find and adapt a fast simulation method of particle paths through porous media based on image processing techniques.

This thesis focuses on the characterization of filters using image processing techniques. Indeed, it is necessary to determine geometric characteristics of the pore space in order to optimize chromatographic filtering regarding a maximum deposition rate of the target cells. Particularly, the curvature and torsion of particle paths through the pore space have a significant impact on the chromatographic process.

In the case of discrete data, there are various approaches to estimate local curvature and torsion from 3D images, which differ in their strategies for handling the discretization effects. The general technique applied in literature is to approximate and smooth 3D curves with analytically known functions and to determine curvature and torsion using the differential geometric formulas. Fitting and smoothing usually depend on various parameters which might be chosen subjectively. More precisely, fitting continuous curves in order to attenuate the effects of discretization can only lead to appropriate curvature and torsion estimation when prior knowledge (required type of function, the order of the functions, etc.) is available. One reason for large estimation errors lies in the fact that the computation of torsion is based on higher order derivatives which are very sensitive towards discretization noise and curve approximation errors. Curvature and torsion estimation is particularly challenging in the cell separation task since no prior information is available for this application. Furthermore, the choice of the appropriate smoothing parameter depends also on the lateral resolution (i.e. the image sampling rate) as well as on the local change of the curve. One aim of this thesis is to improve the curvature and torsion estimation in the presence of discretization noise. This estimation should be widely independent of the curve shape and lateral resolution. To reach these goals, the thesis is structured as follows:

In Chapter 1, the medical context of cell chromatography and cell chromatography itself are explained. The quantities having the most significant influence on the filter efficiency are geometric characteristics of the filter medium, the chemical surface activation, and the flow through the media. Therefore, each of the quantities is discussed in a separate section. One interesting possibility to investigate chromatographic processes is the ➭CT. With an appropriately chosen distance between the sample and the camera as well as with the use of a phase retrieval method, it is possible to increase the contrast between the solid matter, the fluid and the particles. This acquisition technique is based on a setup in which projections of the specimen are taken from different angles. The inverse Radon transform can be applied to these projections to reconstruct the 3D information of the sample. At the end of the chapter, the objectives of the thesis are summarized.

In Chapter 2, a relationship between physical properties of a filter and characteristics obtained from 3D images of the filter will be established. The physical properties like volume, inner surface area, cell size and cell number have their mathematical equivalence in the intrinsic volumes. Therefore, the foundations of integral and stochastic geometry are presented. Besides, the intrinsic volumes there are further characteristics, like the contact distribution function or the percolation probability, which are influencing the chromatographic filter efficiency. However, the most important characteristics under in-vestigation in this thesis are the curvature and torsion of particle paths. A strong torsion of particle paths through the pore space induces a centrifugal force that moves particles to the surface. This effect increases the particle-surface interaction, i.e. particle paths with high torsion serve as micro-centrifuges. The estimation of curvature and torsion of discrete particle paths is a difficult and open task, which leads to the scientific challenge to define accurate estimation approaches.

In Chapter 3, two new estimation methods for curvature and torsion from discrete curves are presented. These methods are closely related to two different discretization schemes (leading to two different data sets) corresponding each to a particular acquisition setup proposed in this thesis. In the first scheme, fast particle paths are simulated by an adapted version of the geodesic dilatation, which is applied to 3D images of porous media. The corresponding curvature and torsion estimation is based on a novel Fourier approximation method. The second scheme is closely related to the discretized paths observed in sequences of 3D images. The corresponding estimation method is based on the discretization of the differential-geometric formulas. Notably, the estimation of the third derivative of a real function is studied, where the estimation methods introduced in this chapter are compared with a reference method commonly used in literature.

In Chapter 4, the estimation methods of curvature and torsion are evaluated on examples of space curves with known curvature and torsion. Special attention is paid to the boundary handling and the influence of the 3D path discretization on the estimation methods. Finally, the mean estimation errors of the presented methods are investigated for two different discretization schemes and different space curves.

In Chapter 5, the experiments for acquiring 3D images of porous media and 3D image sequences of moving particles through porous media are described. Special attention is paid to the description of time-resolved ➭CT setup applied for observing liquid-solid flow through the pore space of chromatographic filters. The acquired 3D images are used to evaluate two porous media regarding their suitability for the chromatographic filtering. Basic geometric characteristics of the chromatographic filters are determined and the new estimation methods are applied for estimating curvature and torsion of particle paths.

Finally, a general conclusion summarizes the main contributions of this thesis and gives perspectives for improvements and further work.

Chapter 1

Chromatographic Cell Separation

In the present chapter, the application, the medical context and the objectives of the thesis are described. The main goal of investigating the geometry of chromatographic filters is to improve the chromatographic efficiency of cell separation, which is very important, notably in cancer diagnostic and treatment. The fundamentals of chromatographic filtering are explained in Section 1.1. The quantities having the most significant influence on the filter efficiency are the geometric characteristics of the filter medium, the chemical surface activation and the flow through the media, respectively, described in Section 1.2, 1.3 and 1.4. The micro-computed tomography (➭CT) is used to investigate the geometry of chromatographic filters. The ➭CT is introduced in Section 1.5. Phase contrast images have to be used to increase the low contrast between filter medium and cells, where phase contrast images can be seen as Laplacian filtered absorption contrast images. In Section 1.5.2 it is explained in detail how the virtual (i.e. not existing) absorption contrast images can be reconstructed with a phase retrieval process. Finally, as described in Section 1.5.3, the inverse Radon transform is used to calculate the 3D information from a sequence of 2D images. Phase retrieval and the inverse Radon transform, considered as two separate image processing steps, are linear transforms and, therefore, they are exchangeable in the image processing chain. Moreover, both transforms can be simultaneously performed in the Fourier domain, which can reduce the computation time.

Medical Context of Chromatography

In 2012, about 14 million people were diagnosed with cancer and each year about 8.2 million people die from the effects of cancer. This means cancer is responsible for 13% of all deaths worldwide [Stewart and Wild, 2014]. The most cancer-related mortalities are caused by metastatic disease [START_REF] Adams | Capture and enumeration of circulating tumor cells from peripheral blood using microfluids[END_REF]. The study of [START_REF] Cristofanilli | Circulating Tumor Cells, Disease Progression, and Survival in Metastatic Breast Cancer[END_REF] shows the relationship between circulating tumor cells and survival rates of the cancer patients. Only 20% of the patients survive more than 18 months if there are more than 5 circulating tumor cells per 7.9 mL blood. With less than 5 tumor cells per 7.9 mL blood the survival rate increases to 75%. Therefore, the isolation of cancer cells from body fluids is an important procedure applied in clinical diagnostics and therapy for cancer detection and treatment. The isolation is performed using techniques like fluorescence activated cell sorting, magnetic bead assisted cell sorting and micro-fluid separation [START_REF] Adams | Capture and enumeration of circulating tumor cells from peripheral blood using microfluids[END_REF], Dharmasiri et al., 2009, Dharmasiri et al., 2010]. These laboratory methods are only applicable for small amounts of cells and sample volumes [Willasch, 2010]. The technique mostly used for separation of cells from human blood circulatory systems, e.g. leukemia cells from healthy ones, or for harvesting cells from industrial incubators is centrifugation. However, industrial centrifuges require high expenditures of investment and maintenance and cause centrifugal damage, which alters the cell surface properties and interior structures up to the DNA [START_REF] Peterson | Bacterial Cell Surface Damage Due to Centrifugal Compaction[END_REF]. Therefore, modern development of technical cell separation is based on alternative approaches such as cell chromatography, where the specific interaction of the target cells with the inner surface is used for the separation. For details see the monographs [START_REF] Plieva | Preparation of macroporous crystructurated gel monoliths, their characterization and main applications[END_REF], Kumar et al., 2007].

In general, chromatographic processes are used to separate composed components like cells, organelles, droplets or particles. In these processes, a suspension of a cell population is pumped through the pore space of a porous material and one or more cell types of the population are bind at the inner surface of the medium. The efficiency of this process depends on the geometry of the medium as well as on the coating of the inner surface with an appropriately activated substance allowing selective cell binding.

In many medical, biotechnical, environmental and synthetic processes, composed components are moved by the convective flow of a surrounding medium. Cells expose a large number of surface-bound molecules to the medium, which are diversified regarding their chemical structure. Classes of cells are characterized by patterns of exposed molecules, e.g. proteins, carbohydrates or lipids. It is often desired to bind one or more classes of cells out of cell mixture. A suitable principle for such a pattern-selective binding can be a combination of long and short reaching interactions between the activated surface of chromatographic filters and the surfaces of target cells . In order to establish a selective surface interaction to the filter, the cells should percolate through the pore space of a micro-structured material driven by the homogeneous flow of the surrounding liquid. A probability of contact events between cells and activated surface exists, depending on structural properties of the porous material and on the long-range interactions. The interaction between cells and a porous medium is sketched in Fig. 1.1.

The contact itself with the unspecialized surface is not sufficient to a selective binding regarding the surface properties of the target specimen. This missing selectivity of surface binding can be encountered by short-range interaction (e.g. between antibodies and cells) with specific binder molecules at the surface of the porous material. The binder molecules provide a high binding constant for the characteristic surface properties of the target component. They can be generated by classical immunology as polyclonal antibodies, by hybridoma cell lines as monoclonal antibodies or by a semisynthetic design of libraries of camelid or shark antibody partial sequences (variable new antigen receptor) Figure 1.1: Chromatographic filtering of green and red cells. The pore size of the medium is multiple times larger than the cell size. The surface is coated with binder molecules (purple hair) that can immobilize green cells. Dark green cells are currently bounded to the surface. The green cells are longer immobilized than the red ones. This effect causes the cell separation.

etc. The conjugation of the selected binder with an anchor molecule is described in common literature [START_REF] Stephanopoulos | [END_REF]Francis, 2011, Novick andRubinstein, 2012]. However, even if such a description seems to resemble a chromatographic separation, important differences to the use of porous media have to be pinpointed. In contrast to chromatography, the mass transfer between the mobile phase and the stationary phase is driven by convection instead of diffusion. The ratio of the kinematic energy of the cells and the free reaction enthalpy is much higher than for molecules that are surface bound in chromatography. This requires a strict control of the attraction and repulsion forces in order to keep the equilibrium. This control can only be maintained when the impact of the three major properties, namely the geometry of the porous phase, the chemical composition of the inner surface and the hydrodynamic flow within the pore space, are well understood and adjusted. For chromatographic filters, the three properties are identical in every volume element of the pore space. The high degree of homogeneity and adjustability of the chromatographic filter gives extensive control on its properties and makes it suitable for "design of experiment" approaches [Kuehl, 1999, Siebertz et al., 2010]. The theoretical approach to identify controlling and noise factors is given by [START_REF] Taguchi | Introduction to Off-line Quality Control[END_REF].

Chromatographic Filters

Most promising candidates for chromatographic filters applied in cell chromatography are macroscopically homogeneous open foams and random fiber systems which are both widely used as particle filters for instance. Therefore, structures of these media are carefully investigated concerning their porosity, specific surface area, tortuosity, finite percolation probability [Ohser et al., 2009] etc. It has been proven that these characteristics significantly influence filtration properties such as flow rate, filtration efficiency and particle holding capability.

A macroscopic homogeneous medium has a probability distribution of its characteristics that is invariant with respect to translations. The macroscopically homogeneity of the pore space is transferred onto the hydrodynamic flow of the cell suspensions driven through the pores. The resulting flow leads to macroscopic homogeneous deposition rates of the cells at the inner surface of the porous medium, which allows for an optimal utilization of each volume element, which is important for high filter efficiency.

Usually, a macroscopic homogeneous foam is produced by a foaming process, where the decomposition of a blowing agent forms bubbles (also said to be the cells of the foam). In an open foam, the walls between neighboring bubbles (the facets of the cells) are removed such that the pore space of the foam is completely percolating. However, improvements of rapid prototyping via 3D printing [START_REF] Fee | 3D printed porous media columns with fine control of column packing morphology[END_REF] open up a new range of possibilities to design chromatographic filters. With 3D printing, the geometry of the inner surface can be easily adjusted in order to allow the selective binding of specific target components. However, efficient filter capabilities also imply accurate geometric models.

Geometric modeling of open foams by computer simulation is usually based on Laguerre tessellations with respect to random hard spheres packing's [Lautensack, 2007a], where the model parameter is the diameter distribution of the spheres. The so-called random sequential addition packing of spheres of gamma distributed diameters have been proven very useful for modeling of polyurethane foams [Ohser et al., 2009, p. 273]. This means that the random sequential addition packing depend on only the two parameters of the gamma distribution. Two further model parameters are necessary in order to describe the shape of the struts' cross sections. For instance, in [START_REF] Kanaun | Representative volume element and effective elastic properties of open cell foam materials with random microstructures[END_REF] the cross section is described by a parametric function along the strut axis.

This geometric model can be fitted to the real structure of partially open foam by varying the model parameters and comparing the geometric characteristics of the model with the estimates from a tomographic image of the foam.

The foam model can be generated in such a way that it is periodic with respect to the bounded periodic unit, e.g. shown in Fig. 1.2a. The periodic patterns are notably at the basis of the simulation of fluid flow through the pore space of the model as well as of the construction of large filters by periodic extension of the bounded periodic unit, seen in Fig. 1.2b. Such a designed foam is not macroscopically homogeneous anymore, as the consequence of the periodicity, there exists straight particle path through the pore space.

Notice that an appropriate geometric foam model is the Weaire-Phelan foam, shown in Fig. 1.2. It is the smallest and thus most simple bounded periodic unit [Weaire and R.Phelan, 1994]. Its adaption to real foams is described in [START_REF] Redenbach | A random weairephelan foam. 8th International Conference on Stereology and Image Analysis in Material Science[END_REF] and [Ohser et al., 2009, p. 266].

Geometric modeling of fiber systems is based on the Poisson's straight line distribution, e.g. Poission straight lines were developed in [Matheron, 1975] and Poission thick lines (cylinders with any cross-section) by [Jeulin, 1991]. These models describe overlapping fibers with infinite length. Important characteristics like specific surface and contact distribution can be explicitly calculated for convex fiber cross sections [Hoffmann, 2008, Spiess and[START_REF] Spiess | [END_REF]. A method to model compact, curved and nonoverlapping fiber systems is introduced in [START_REF] Altendorf | Random-walk-based stochastic modeling of three-dimensional fiber systems[END_REF].

Chemical Activation of the Inner Surface

Synthetic polymers and artificial resins such as polyurethane possess an inherent low surface energy that limits the application for most bio-coatings. Thus, giving a surface a cell interacting character is challenging because the direct application of hydrophilic surface coatings is only accessible to a limited extent. In most cases, the poor polarizable surface needs to be pretreated with low pressure plasma activation to make the surface attractive for subsequent ligand functionalization [START_REF] Sanchis | Surface modification of a polyurethane film by low pressure glow discharge oxygen plasma treatment[END_REF]. This can be performed by physical and chemical adsorption of ligands that promote cell attachment. One example is the application of synthetic ligands such as Poly-L-Lysine or Poly-D-Lysine, commonly used to coat tissue cell cultures by introducing positive charges to the surface. Another example is the application of proteins containing specific sequence motifs such as the arginine-glycine-aspartate (known as RGD-sequence) [Bellis, 2011] that interact with membrane-spanning receptor proteins called integrins and extracellular matrix components such as collagen and fibrin [D' Souza et al., 1991]. Both described ways are suitable to create a bio-functionalized surface that is attractive for mostly all types of cells.

Organosilanes are known for their function as adhesive agents that promote the bonding of organic layers to inorganic layers under formation of surface monolayers. To enable the deposition of organic layers, most surfaces that possess an inherent low surface energy such as polyurethane or a variant thereof, need to be pretreated with low plasma activation or piranha solution to create surface exposed-OH groups. The created hydroxyl groups can further undergo condensation reactions with the silanol groups of the silanes under formation of the surface monolayers. The silanes used for that application contain terminal functional groups that can be seen as anchors for biomolecule conjugation (organosilanes). In order to create bio-functional surfaces that feature a high degree of homogeneity, the cell-binding molecules can be covalently coupled through a variety of chemical strategies in a very defined way. In this context, the term cell-binding molecule means any protein featuring specific properties that interacts with the targeted cell type.

A second option is a composed activated monolayer based on amphiphilic molecules covering the hydrophobic surface of the structure. This stable oriented monolayer exposes the hydrophilic residues of its molecular constituents to the liquid phase, while the hydrophobic residues are in contact with the hydrophobic polymer. Second, the amphiphilic monolayer serves as anchoring phase for the ligands which are selectively binding the suspended target cells. Complex biogenic ligands can be fixed to the covering monolayer using mild conditions. An important feature of the monolayer bound binders is their ability for lateral diffusion. The mobility of selective protein binders and of charge carrying binders allows for the accumulation of both binder types into the geometric interface area of molecular interaction in a given time of direct interaction. This yields into a maximum number of selective molecular interactions as well as balanced ionic interactions between target cells and activated surface. The described activated monolayer surface does not require demanding covalent coupling chemistry. Complex and sensitive binder molecules can be fixed to the surface in an even distribution. Finally, the simple process of binder fixation allows for serial variation of binder surface concentration and composition. Using these properties, the optimal design of a surface activation can be achieved by systematic variation. In this context, the term hydrophobic or amphiphilic anchor molecule comprises any lipid such as fatty acid, triglyceride, steroid alcohols (cholesterol), phospholipid, sphingolipid, glycerolipid, lipopolysaccharide or synthetic derivatives, thereof 1.4. Hydrodynamic Flow through Pore Space that features a functional group suitable for bioconjugation. Within the use of these anchor molecules the need of rendering the surface hydrophilic through the creation of surface exposed hydroxyl groups is not necessary. The cell binder molecule can be coupled directly in a one-step reaction to the anchor molecule which assembles itself chemically by principles of molecular recognition under formation of a self-assembled monolayer. This will yield in a high degree of topological homogeneity needed for the herein described chromatographic approach.

Hydrodynamic Flow through Pore Space

The transport of cells through a porous phase and to the activated surface is a convective process. A hydrodynamic flow path equipped with torsion allows to hit the concave surfaces for the floating cells. The convective flow through the pore space induces the boundary conditions for the selective initial adhesion of target cells to the activated surface. For an initial adhesion the cells have got in direct contact with the activated surface. The floating cell might drift out a zone of high volume flow velocity into a volume element with low volume flow nearby the activated surface. Thereby, it might be attracted into geometric contact by mechanic impact of cell mass, by micro-turbulent flow and/or electrostatic attraction. The cells in geometric contact might have an average residence time at the selective surface which is required to establish the specific molecular interaction between the binder molecules of the activated surface and the socalled reporter groups which are exposed at the surface of the target cells. Considering the complex pattern of convective transport of cells and the multistep binding process, it is very difficult to predict the conditions for most selective or most productive binding of target cells. Thus, the design of experiments approach is required for the development of an efficient cell separation device. The prerequisite for a successful design of experiment approach is the uniformity of microscopic process conditions in all volume elements of the stationary phase which are flown through by the suspended target cells. As described above, the most suitable geometry of the stationary phase is a defined isotropic structure as a multiplication of 3D unit elements, e.g. a face removed Weaire-Phelan foam (see Fig. 1.2) or an open foam with a narrow random distribution of pore diameters. The macroscopic flow of cells through the stationary phase has to be invariant for every 3D unit element. The flow can be generated either by moving a fluid medium through a closed device which contains the stationary phase or by moving the stationary phase in an open device through a fluid medium. In both cases the flow conditions of the entrance of the stationary phase have to be very similar to the flow conditions during percolation and at the outlet. This setup cannot be realized by a chromatographic column which has an inlet and an outlet fluid connection with diameters far smaller than that of the column bed or stationary phase. Large scale chromatography columns are often equipped with flow distributing spare parts at the entrance. These functional elements are efficient for solutions of molecules in which diffusion controls the mass transfer. In contrast, the movement of suspended cells floating through a stationary phase is convective. An en-trance flow which passes narrow inlets or even flow distributors results in an irregular movement of cells characterized by jet streams of high cell density and volumes in between with low cell densities. This uneven and not controllable flow scheme renders any tuning of the cell binding process impossible. As a simple alternative, a set of parallel capillaries can be used to generate an even and controllable flow. But capillary structures have typically very low property of cells to surface contact event and therefore a low binding capacity for cells.

Investigation of Chromatographic Filters by Tomography

The classical investigation of chromatographic filters is focused on the cell adhesion realized by the chemical activation of the surface and the composition of the mobile phase. The investigation setup is normally simplified to one-dimensional flow over thinlayered stationary phase, which can be easily observed by microscopes [START_REF] Oh | Analytical cell adhesion chromatography reveals impaired persistence of metastatic cell rolling adhesion to P-selectin[END_REF]. However, in this thesis, the influence of the geometry of the pore space and the 3D flow through it is the matter of investigation. Therefore, microscopes are no longer possible and the more sufficient method of computed tomography is used for the investigation.

Tomography is the process to create 3D images by using projections of a rotating 3D object. When discrete calculations are mandatory for the tomography process, it is called computed tomography (CT). In this work, the 3D images have at least a resolution more accurate than 10 ➭m. So, to distinguish the data under interest from CT-scans in medical application with a resolution of millimeter, the process is specified as ➭CT.

Tomography can be classified by the shape of the beam. There are ➭CT with cone beam and collimated beam. Cone beam ➭CT is used in laboratory CT scanners with a lens less magnification. The magnification factor depends only on the distance between the sample and the radiation detector. Collimated beams are generated by a synchrotron. A synchrotron produces monochromatic collimated coherent radiation with a high photon flux, which avoids hardening artifacts and provides a high signal to noise ratio images. The synchrotron radiation enables the selection of different wave lengths and the use of phase contrast images. The increasing brilliance of X-ray sources available nowadays gives via the tomographic approach access to the fourth dimension: the time [START_REF] Rack | On the possibilities of hard X-ray imaging with high spatio-temporal resolution using polychromatic synchrotron radiation[END_REF], Mokso et al., 2013, Maire and Withers, 2014, Di Michiel et al., 2005, Zabler et al., 2012]. This enables the application of time-resolved ➭CT for scanning processes.

Synchrotron radiation is exclusively used in this work.

Synchrotron Light Source

The 3D images used in this thesis are produced at the European Synchrotron Radiation Facility (ESRF) on beam line ID19. A synchrotron consists of three parts shown in Fig. 1.3. At first, the electrons will be emitted and accelerated to near-light speed in the linear accelerator. Afterwards, the electron energy will be further increased to 6 GeV in the booster synchrotron. In the third step, the electrons are passed to the storage ring [Facility (ESRF), 2009]. Due to the high speed of electrons in the storage ring, relativistic effects take place and the low speed dipole radiation of electrons becomes to a narrow cone radiation pointing to the traveling direction [Schleede, 2013] and causing a high collimation.

Phase Contrast and Phase Retrieval

The absorption contrast is typically used in X-ray tomography. In the case investigated in this thesis, the cell simulators and the porous media are too small to weaken the radiation significantly to get interpretable data. For this reason phase contrast imaging is used. The phase shift is caused by a higher phase velocity in a medium represented by a complex refractive index n = 1δiβ.

(1.1)

The real part 1δ stands for the elastic Rayleigh scattering caused by a forced oscillation of bounded electrons, which are reemitting the forcing wave [START_REF] Pedrotti | Optik für Ingenieure, Grundlagen[END_REF].

The inelastic Compton scattering accounts for the imaginary part -β. It is caused by electrons, which absorb a part of the primary radiation and scatter a lower powered radiation. The refractive index is wave length dependent, whereas the material has an absolute refractive index n > 1 for visible light, the factor decreases to values smaller one for high energetic radiation. This causes a defocusing of a beam passing a convex body. Phase contrast is only available for spatial coherent radiation. Such radiations are produced by synchrotrons. Let us consider a single ray parallel to the x 3 -axis and pass-1.5. Investigation of Chromatographic Filters by Tomography ing through the medium with the time and spatial depending energy E(x 3 , t). In such condition the wave equation

E(x 3 , t) = Êe i(ωt-nk 0 x 3 ) = Êe i(ωt-k 0 x 3 +δk 0 x 3 +iβk 0 x 3 ) (1.
2) [Schleede, 2013] holds. Eq. (1.2) factories itself in the energy E in = Êe i(ωt-k 0 x 3 ) of the vacuum wave, the phase shift e ϕ = e δk 0 x 3 and the amplitude decay e -βk 0 x 3 , where Ê is the amplitude, k 0 = 2π λ is the wave number, ω is the angular frequency, t is the time and x 3 is the position. Because, on the one hand, the phase shift and amplitude decay are time independent and, on the other hand, in the following only the spatial changing of the wave is considered, the time depending factor of the vacuum wave is omitted. The x 3 axis is chosen in such a way, that the right boundary of the object is on position x 3 = 0. Thus, as sketched in Fig. 1.4, the position of the left object boundary corresponds to the negative thickness of the object. The energy of the wave directly after the sample is denoted by E 0 . The attenuation I 0 /I in of the material can be obtained from

I 0 (x 3 ) I in = |E 0 (x 3 )| 2 |E in (x 3 )| 2 = e -2k 0 βx 3 , (1.3)
where I 0 (x 3 ) is the intensity of the attenuated wave and I in is the initial wave intensity. I in is independent of x 3 because the absolute value of a planar wave is calculated. Now, an X-ray beam is considered, which is parallel to the x 3 axis with a vector r = x 1 x 2 lying in the plane perpendicular to the optical axis x 3 and indicates the position of one point in the plane. Then, for a given object thickness d(r), the intensity I 0 is (1.4) which corresponds to the Beer-Lambert law I(r, 0) = I in e -µd(r) , (1.5)

I 0 = e -2k 0 βd(r) I in ,
ϕ(r, 0) = -k 0 δd(r), (1.6)
where I in is the input intensity, µ is the attenuation coefficient of the material insight the object and δ is the real part of the deviation of the material's refractive index from one [START_REF] Wilkins | Phase-contrast imaging using polychromatic hard X-rays[END_REF]. Than d(r) is proportional to the projection of the object onto the plane x 3 = 0. A comparison of formula (1.5) with (1.4) gives µ = 2k 0 β.

The phase shift cannot be measured directly, but, there are different methods [Banhart, 2008] to modulate the intensity by quantities related to the phase shift. The most common and useful method for synchrotron radiation is inline phase contrast. Inline phase contrast means that the phase shift of the sample is encoded as intensity interference pattern appearing at a particular distance. To calculate the intensity pattern the propagation formula is used

∇ I 0 (r)∇(ϕ(r) = -k 0 ∂I 0 (r, x 3 ) ∂x 3 , ∇ = ∂ ∂x 1 ∂ ∂x 2
(1.7)

1.5. Investigation of Chromatographic Filters by Tomography [ [START_REF] Rytov | Wave Propagation Through Random Media[END_REF], where ∇ is the two dimensional nabla operator. Substituting Eq.

x 3 ≤ -d x 3 = 0 (x 3,0 ) 0 < x 3 ≪ d 2 λ (x 3,1 ) x 3 ≈ d 2 λ x 3 ≫ d 2
(1.4) and then Eq. (1.6) into Eq. (1.7) yields [START_REF] Paganin | Simultaneous phase and amplitude extraction from a single defocused image of a homogeneous object[END_REF], where ∆ is the Laplacian operator. Notice, that the identity

1 -2k 0 β I in ∆ 1 d(r) e -2k 0 βd(r) ϕ(r) = -k 0 ∂I(r, x 3 ) ∂x 3 , δ 2β I in ∆e -2k 0 βd(r) = -k 0 ∂I(r, x 3 ) ∂x 3 [
δ∇(e -2k 0 βd(r) ∇d(r)) = - δ 2k 0 β ∆e -2k 0 βd(r)
is used during substitution. The derivation with respect to x 3 at the positions x 3,0 = 0 and x 3,1 can be approximated with

∂I 0 (r, x 3 ) ∂x 3 ≈ I 1 (r) -e -2k 0 βd(r) I in x 3,1 -x 3,0
and it yields [ [START_REF] Paganin | Simultaneous phase and amplitude extraction from a single defocused image of a homogeneous object[END_REF]. The result is a Laplacian filtered absorption contrast image at the distance x 3,1 . Notice, that this observation is valid for the near field regime, only. Fig. 1.4 shows the diffraction pattern on different distances x 3,1 . The absorption contrast images can be acquired directly after the sample at position x 3,0 . The wave interference shown in Fig. 1.5 is caused by the interference of the phase shifted wave inducted by the sample and the unrefracted planar wave. This results in a very simple setup for phase contrast images. The X-ray has only to propagate a distance in free space before it is focused to a CCD chip. If the distance is chosen properly, phase contrast images enhance the object edges and homogeneous regions are suppressed. Such images can be seen as the second derivative of the absorption contrast images, which is mathematically described by the Laplacian operator in Eq. (1.8). However, the segmentation of edges enhanced images is difficult, as seen in Fig. 1.6a. It is preferable to retrieve the phase for an easy binarization, as seen in Fig. 1.6b.

-δx 3,1 2βk 0 ∆ + 1 e -2k 0 βd(r) = I 1 (r) I in (1.
The wave propagation formula has to be considered to understand the phase retrieval process and, thus, the relation between absorption contrast and phase contrast images (see Fig. 1.4). The absorption contrast image can be acquired on position x 3 = 0. To retrieve the phase the projected thickness of the sample has to be calculated from the near field contrast image. Therefore, Paganin developed the following formula [ [START_REF] Paganin | Simultaneous phase and amplitude extraction from a single defocused image of a homogeneous object[END_REF], where F ξ denotes the two-dimensional Fourier transform respect to the vector r, and Fξ is the corresponding co-transform with respect to the variable ξ. The division by ξ 2 in Fourier domain corresponds to an inverse Laplace filtering in spatial domain, which retrievals the complete phase information. Interestingly, only one image I(r, x 3 ), the input radiation I in and the refractive index n are needed to calculate the projected thickness of a homogeneous sample. One may think of the algorithm as an aberration-free virtual lens, which brings into focus the necessarily defocused images [START_REF] Paganin | Simultaneous phase and amplitude extraction from a single defocused image of a homogeneous object[END_REF].

d(r) = - 1 µ ln Fξ µ (F ξ I(r, x 3,1 )) (ξ) I in x 3,1 δ ξ 2 + µ (r) (1.9) (a) (b)

The Inverse Radon Transform

The Radon transform is an integral transform, where the integrals of parallel intersection lines with a 3D object f (x) are projected onto a plane. These projection images Rf done for multiple angles ϑ are created by CT. The inverse Radon transform is used to get back the 3D information from a sequence of 2D projections.

The back projection B is given by

BRf (x) = π 0 Rf (x 1 cos ϑ + x 2 sin ϑ, ϑ) dϑ, x = x 1 x 2 ∈ R 2 ,
(1.10) [Toft, 1996], where ϑ is the angle of the normal of the straight line. Notice, that in general the backprojection does not retrieve the original signal f = BRf . In order to explain tomographic reconstruction, the Fourier slice theorem is applied, which says

Rf (s, ϑ) = 1 2π ∞ -∞ F s Rf (-σ sin ϑ, σ cos ϑ)e -2πisσ dσ, s ∈ R, ϑ ∈ [0, π),
where -σ sin ϑ and σ cos ϑ are the polar coordinates and F s is the one-dimensional Fourier transform with respect to the radial coordinate s. Using this, the Radon transform Rf is firstly considered in frequency domain,

F s Rf (σ, ϑ) = ∞ -∞ Rf (s, ϑ) e 2πisσ ds (1.11) = f (σ cos ϑ, σ sin ϑ), σ ∈ R, ϑ ∈ [0, π).
Now, the 2D inverse Fourier transform can be applied to get f = F f ,

f (x) = 1 (2π) 2 R 2 f (ω) e -2πixω dω = 1 (2π) 2 ∞ -∞ ∞ -∞ f (ω 1 , ω 2 ) e -2πi x 1 x 2 ω 1 ω 2 dω 2 dω 1 = 1 (2π) 2 2π 0 ∞ 0 f (σ cos ϑ, σ sin ϑ) e -2πi x 1 x 2 σ cos ϑ σ sin ϑ σ dσdϑ = 1 2π π 0 1 2π ∞ -∞
f (σ cos ϑ, σ sin ϑ) e -2πi(x 1 cos ϑ+x 2 sin ϑ)σ |σ| dσ dϑ.

(1.12)

Equation (1.12) can be reformulated with the back projection formula (1.10) and the one dimensional Fourier transform (1.11),

f (x) = 1 2π B Fσ (abs • F s Rf ),
where abs σ = |σ| is the absolute value of σ. The multiplication of function abs in frequency domain corresponds to a convolution with the one-dimensional function ( F abs)(x) in spatial domain

f (x) = 1 2π B (Rf ) * F (abs) (x).
(1. 13) In image processing, convolution is equated with filtering. That is the reason why this inversion is called filtered backprojection [Deans, 2007]. Of course, since the function abs 1.6. General Aims of the Thesis is not measurable, its Fourier transform does not exist, but for the moment it is assumed that it exists. This is made only for demonstration reason. In the discrete case, where periodicity is implied, the periodic extension of the function abs is applied which makes it irrelevant that it is not measurable. Therefore, in image processing the nonexistent function F (abs) can be replaced by the discrete function

g k =        1 4∆s , k = 0 0, k = 2, 4, 6, . . . , m 2 -2 -1 π 2 k 2 ∆s , k = 1, 3, 5, . . . , m 2 -1 (1.14)
where ∆s is the step width and m ≥ n is a multiple of four. The transfer function g(σ) of a one-dimensional filter mask (g 0 , . . . , g m-1 )) is

g(σ) = 1 4 - 2 π 2 m/4 k=0 cos 2π(2k + 1)σ (2k + 1) 2 , σ ∈ R, (1.15)
which is the periodically extension of abs σ. Now, switching again to the continuous case, it is possible to calculate the inverse Radon transform in the spatial domain, without the need of the Fourier transform

f (x) ≈ 1 2π B (Rf ) * g (x), (1.16)
where g(x) denotes the (not existing) inverse Fourier transform of abs. This means, the Fourier transform is helpful to derive the formulas for tomographic reconstruction, but it is not needed in its algorithmic implementation. Nevertheless, since the Paganin filtering is done in the frequency domain, it is useful to calculate the inverse Radon transform and the phase retrieval simultaneously in the frequency domain, as suggested in [START_REF] Burvall | Phase retrieval in X-ray phase-contrast imaging suitable for tomography[END_REF].

General Aims of the Thesis

Chromatographic cell filtering is a technique appearing as a proper alternative to fluorescence activated cell sorting or centrifugation. Appropriate filters for cell separation are the widely used open foams or fiber fleeces. However, the underlying question which emerges with this technology is:

"How a filter has to be designed to maximize the filter efficiency?"

In this thesis, one step is made to design efficient filters. It is necessary to determine geometric characteristics of the pore space in order to optimize chromatographic filtering regarding a maximum deposition rate of the target cells. Characteristics widely used in industrial quality control of porous media are the porosity and the specific surface area. A further characteristic of foams widely used from foam producers and foam consumers is the mean pore size expressed in terms of the mean number of cells per inch (pores per inch, PPI). Important characteristics of fleece filters are the specific fiber length and the fiber direction distribution. There are further quantities suggested in the literature, e.g. the density of the integral of mean curvature, the specific Euler number [Ohser and Schladitz, 2009], the spherical contact distribution [START_REF] Muche | Contact and chord length distributions of the Poisson-Voronoi tessellation[END_REF], the spherical granulometry [START_REF] Chiu | Stochastic Geometry and Its Applications, 3rd Edition[END_REF], the autocorrelation function (or the spectral density), the percolation distribution [START_REF] Ohser | Estimation of the probability of finite percolation in porous microstructures from tomographic images[END_REF], the tortuosity [START_REF]A Novel Curvature Based Algorithm for Automatic Grading of Retinal Blood Vessel Tortuosity[END_REF] etc. Many of these characteristics depend on each other. For example, the density of the integral of the mean curvature is up to a constant factor the same as the specific fiber length. A unique characterization of the pore space would be possible with the Choquet capacity [Wasserman and Kadane, 1990, mat, ], but a fully description of the pore space implies the intersection with the set of all compact bodies, which is practically not feasible.

In the present thesis, the particle paths of a two-phase flow through the pore space of porous media are studied, where the solid phase of the flow consists of approximately spherical particles. Up to translation and rotation, each path is uniquely characterized by its curvature and torsion [Pressley, 2010] and, thus, the distribution of the curvature and torsion of the particle paths are important characteristics of the two-phase flow. In particular, the torsion distribution of the particle paths depend on the orientation of porous media, i.e. it depends on whether the turbulent flow induced by the pore space is either left rotating or right rotating.

There are various techniques to investigate chromatography. An obvious and very simple one is the comparison of the cells' densities of the in-and output flow. However, such experiments cannot give deeper insights into the processes at the microscale. A more promising approach is the observation of particles flowing through the media. An appropriate way to achieve this goal is high speed ➭CT. Unfortunately, the required time resolution to detect full speed particles cannot be reached by state of the art CTscanners. Consequently, only slow particles near the surface can be observed. Two-phase flow simulations offer the opportunity to investigate faster particles. However, a flow simulation requires deep knowledge about the composition and viscosity of the suspension and knowledge about the material, surface layer, stiffness, chemical activation, used antibodies etc. Moreover, these properties depend on temperature, flow rates and pressure and the calculation effort of a 3D flow simulation on 3D images of about 16 GB each is immense. Therefore, this thesis is focused on the calculation of particle paths directly from the geometric characteristics of the porous media. This allows for evaluating porous media independently of chemical activation, flow conditions and used material.

The aim of this thesis is to develop methods to characterize porous filters regarding their suitability for cell chromatography that only depend on the geometric characteristics of the filters and the particle size. The first general aim is to establish the link between the geometric structure and the characteristic of porous media. The second general aim lies in the design of an experimental set-up for the investigation of the capturing process of target particles at the inner surface of porous media by three-dimensional image processing. Finally, two different kinds of porous media should be exemplary investigated and compared with the developed methods regarding their suitability for cell chromatography.

Chapter 2

Characteristics of Porous Media

Classical methods to characterize open porous media by geometric features are based on counting nodes, struts, pores or fibers via microscopy. The volumes and surfaces of samples of such media are typically determined by the Archimedean density principal and BET (Brunauer, Emmett and Teller) method [START_REF] Harald | Elastischer Schaumstoff auf Basis eines Melamin/Formaldehyd-Kondensationsproduktes und seine Verwendung[END_REF], Stork et al., 2005, Klank, 2004]. A wide range of equipment and methods are used for these classical characteristic determinations. The counting of struts, nodes and pores is very errorprone because it is usually done by humans and the three-dimensionality of the foam is insufficiently considered. Features like curvature and torsion of particle paths cannot be directly measured with such approaches.

New tools of 3D image processing and analysis are developed in order to characterize porous media with respect to their application in chromatography. In this chapter, a relation between physical measurements and characteristics obtained from 3D images will be established. The physical characteristics like porosity, inner surface area, cell size and cell number have their mathematical equivalence in the intrinsic volumes. Therefore, the foundations of integral and stochastic geometry are presented in Section 2.1. Besides, the intrinsic volumes there are further characteristics, which are influencing the chromatographic filter efficiency like the contact distribution function or the percolation probability. Both characteristics are presented in Section 2.2.1. In this thesis, two kinds of porous media are investigated, namely open porous foams and random fiber fleeces. The intrinsic volumes can have different meanings for open foams and fiber fleeces and there are characteristics that only apply for one of the both. These differences are explained in Sections 2.2.2 and 2.2.3. However, the most important characteristics detailed in Section 2.3 are curvature and torsion of particle paths through the media. A strong torsion of particle paths through the pore space induces a centrifugal force that moves particles to the surface so that the particle-surface interaction increases, i.e. particle paths with high torsion serve as micro-centrifuges. However, it is very challenging to estimate the torsion of particle paths, which leads to the scientific challenge in Section 2.4.

Firstly, some terms of integral geometry and stochastic geometry which will be frequently used throughout this thesis are introduced. For details on integral geometry, the reader may refer to [Schneider, 1993], while a sound introduction to stochastic geometry can be found in [START_REF] Schneider | Stochastic and Integral Geometry[END_REF], Chiu et al., 2013, Matheron, 1975]. The aim of this section is to introduce some fundamental characteristics (features) describing geometric objects (particles, pores, etc.) or constituents of microstructures.

Characteristics of Objects

Let X be a convex body (in the language of image processing an object or a connected component), which is a compact and convex subset of the 3-dimensional Euclidean space, X ⊂ R 3 . Its volume V (X) is simply the Lebesgue measure [Kestelman, 1960] of X. It is well known from integral geometry that the volume is a characteristic belonging to the class of so-called intrinsic volumes. Up to multiplicative constants, these are the just introduced volume V , the surface area S, the mean width b and the Euler number χ. The definitions of S, b and χ trace back to that of V . Let X ⊕ B r be a parallel set of X, where ⊕ denotes the Minkowski addition and B r is the ball of radius r centered at the origin. Clearly, X ⊕ B r is the dilation of X with the structuring element B r since the center of the ball B r is symmetric. Then the Steiner formula [Gray, 2004] states that the volume of the parallel set is a polynomial of 3rd degree in r, where the coefficients are even the four intrinsic volumes (up to multiplicative constants).

The volume seems to be the simplest characteristic which can easily be estimated from 3D images and intuitively be interpreted in applications. This is due to the fact that the volume is additive, i.e.

V (X 1 ∪ X 2 ) = V (X 1 ) + V (X 2 ) -V (X 1 ∩ X 2 )
for all bodies X 1 and X 2 . Furthermore, the volume is invariant with respect to translations, isotropy, 3-homogeneous and continuous. By means of the Steiner formula these properties are inherited to the other intrinsic volumes. This makes the intrinsic volumes interesting for application in image analysis. In particular, the additivity allows to extend the intrinsic volumes on non-convex bodies. From the inclusion-exclusion principle it follows, that the intrinsic volumes are well defined for all finite unions of convex bodies -the class of finite unions of convex bodies, also said to be the class of polyconvex sets.

The four intrinsic volumes are unique in the sense of Hadwiger's characterization theorem, which states that every additive, translation-invariant, isotropic and continuous characteristic can be represented as a linear combination of the intrinsic volumes [Hadwiger, 1957]. In other words, the intrinsic volumes form a basis of the space of features. This fact is of high importance also for application since any further characteristic having the above properties and is estimated by image analysis would be redundant.

The mean width is also known as the mean caliper diameter or the mean Ferret diameter. Firstly note, that for convex polyhedral the mean width can be computed by means of Santaló's formula. Let P be a polyhedron consisting of m edges with length ℓ k , where k is the edge number in [1, m] and with the angle γ k between the outer normals of the faces hitting the k-th edge. Then it holds [Santaló, 1976, p. 226]. This formula can be extended onto non-convex polyhedral, which allows to compute the mean width from a polygonalization obtained by surface rendering.

b(P ) = 1 2π m k=1 ℓ k γ k ,
Moreover, for compact sets with a sufficiently smooth surface, one can give a relationship between the mean width b and the integral of the mean curvature M (Germain's curvature). More precisely, if radius r is strictly positive so that the convex body is morphologically open with respect to B r , the two principal curvatures exist for every surface elements of the body, and M is the mean of the principal curvatures integrated over all surface elements. For example, the surface area of parallel set X ⊕ B r is smooth everywhere. This parallel set X ⊕ B r converges to X as r → 0, i.e. dist (X ⊕ B r , X) → 0 as r → 0, where dist is the Hausdorff distance [Hausdorff, 1949]. One gets

2π b(X) = lim r→0 M (X ⊕ B r )
for all convex bodies X ⊂ R 3 . The last formula can be generalized for polyconvex sets with non-empty interior. By X

• B r = (X ⊕ B r ) ⊖ B r and X • B r = (X ⊖ B r ) ⊕ B r
the morphological closing, respectively, the morphological opening of X with B r are denoted, where ⊖ is the Minkowski subtraction. Let X be polyconvex, than the sets (X • B r ) • B r and (X • B r ) • B r have smooth surfaces and converge to X for vanishing r. It follows

2π b(X) = lim r→0 M (X • B r ) • B r = lim r→0 M (X • B r ) • B r
for all polyconvex sets X ⊂ R 3 with non-empty interior. Thus, throughout this thesis the setting

M (X) = 2π b(X) is used.
Similar relationships exist between the Euler number χ and the integral of the total curvature K (Gaussian's curvature), which is defined as the product of the principal curvatures integrated over all surface elements. For all polyconvex sets X ⊂ R 3 with non-empty interior it yields

4πχ(X) = lim r→0 K (X • B r ) • B r = lim r→0 K (X • B r ) • B r .
Analogously, K(X) = 4πχ(X) is used for all polyconvex X with non-empty interior.

Finally, the shape factors

̺ 1 = 6 √ π V (X) S 3 (X) , ̺ 2 = 48π 2 V (X) M 3 (X) , ̺ 3 = 4π S(X) M 2 (X)
of polyconvex sets X are introduced. The shape factors are normalized in such a way that

̺ 1 = ̺ 2 = ̺ 1 = 1 if X is a ball.
The object shape deviates from a ball shape when the shape factors deviate from 1 [Ohser et al., 2009, p. 163]. From the isoperimetric inequalities given in [Schneider, 1993, pp. 318, 325] it follows that 0 < ̺ 1 ≤ 1 is valid for all polyconvex sets and 0 < ̺ 2 , ̺ 3 ≤ 1 is valid for all convex sets.

The Estimations of V , S, M and K from binary images B of the set X are based on the additivity, the translation invariance and the isotropy of these quantities. As a consequence, one can estimate them from numbers of local pixel configurations. More precisely, let L 3 = aZ 3 be a homogeneous point lattice with the spacing a > 0. The sampling X ∩ L 3 is the set of the foreground pixels of the binary image, while X c ∩ L 3 forms the background. By h = (h l ) the vector of the numbers h l of 2 × 2 × 2-pixel configurations, l = 0, . . . , 255 is denoted. For a unique indexing of the pixel configurations a 3D filter mask

M = 32 16 128 64 , 2 1 8 4
can be introduced, where the gray square marks the position of the origin. The filtering (i.e. the discrete convolution) G = B * M of the binary image B with the mask M yields to an 8-bit gray-tone image whose pixels are the indices of the pixel configurations of B. As a consequence, the vector h is the so-called gray-tone histogram of G, where h l is the total number of pixels of G with the gray-tone l. Finally, pictograms are introduced in order to depict pixel configurations. For example, the pixel configuration ξ 135 with the index l = 135 is depicted as . Using the above mentioned notations and the convention ∞ • 0 = 0, each of the quantities V , S, M and K can be estimated from a k hw with k = 3 for the volume V , k = 2 for the surface area S, k = 1 for the integral of the mean curvature M and k = 0 for the integral of the total curvature K (or equivalents of estimating the Euler number χ). Here w = (w l ) is the vector weights w l , l = 0, . . . , 255, given in Table 2.1 for representatives of the 22 pixel classes of configurations which are equivalent up to rotations and reflections. This method of estimating the Euler number traces back to that of Jean Serra, see e.g. [Serra, 1982]. The weights for S and M are computed based on discretization of the corresponding Crofton formulas [START_REF] Schladitz | Measurement of intrinsic volumes of sets observed on lattices[END_REF]Schladitz, 2009]. The estimation of the integral of the total curvature is based on the application of the inclusion-exclusion principle [START_REF] Nagel | An integral-geometric approach for the Euler-Poincaré characteristic of spatial images[END_REF], Ohser et al., 2002, Ohser et al., 2003, Ohser and Schladitz, 2009]. Table 2.1 presents the weights of χ for the 6-neighborhood and the 26-neighborhood, respectively. Notice that estimates of the Euler number strongly depend on the neighborhood of the pixels, where Table 2.1: The weights w l for estimating the volume V , the surface area S and the fiber length M/π resp. the Euler number χ = K/4π. Different surface area weights from a) [START_REF] Schladitz | Measurement of intrinsic volumes of sets observed on lattices[END_REF], b) [Lindblad, 2005] and c) [START_REF] Ziegel | Estimation of the surface area and surface area measure of three-dimensional sets from digitizations[END_REF] are presented. The weights in [START_REF] Ziegel | Estimation of the surface area and surface area measure of three-dimensional sets from digitizations[END_REF] depend on the parameters r and s which can be chosen arbitrarily in the intervals 0 ≤ r ≤ 1/2 and 1.663 ≤ s ≤ 1.745.

l Config. V S M χ a) b) c) 6er-N. 26er-N. 0 ξ 0 0 0 0 0 0 0 0 1 ξ 1 0.125 0.
pixel neighborhoods on 3D images are uniquely characterized based on the concept of adjacency systems developed in [START_REF] Nagel | An integral-geometric approach for the Euler-Poincaré characteristic of spatial images[END_REF], Ohser et al., 2003]. This concept is the key for a systematic study of pixel adjacencies and their complementarity in higher dimensions. It also allows for a clear description of the 14.1-and 14.2-adjacencies for 3D images which are the only self-complementary adjacency systems. However, the 14.1and 14.2-adjacencies are anisotropic and, thus, the corresponding weights for the Euler number cannot be presented for representatives of the 2×2×2-pixel configurations, only.

There are a lot of further approaches of estimating the surface area. In Table 2.1, the weights for the surface area are compared with those from [Lindblad, 2005], which are optimal in the sense that they minimize the mean square error of the surface area estimation, see also [START_REF] Legland | Computation of Minkowski measures on 2D and 3D binary images[END_REF]. The work in [START_REF] Ziegel | Estimation of the surface area and surface area measure of three-dimensional sets from digitizations[END_REF] presents a parametric estimation of the surface area which minimizes the maximum error for vanishing pixel size, i.e. for increasing lateral resolution. As pointed out in [START_REF] Ohser | Improved estimation of the fiber length from 3-dimensional images[END_REF], with a sufficiently large lateral resolution, the relative error of estimating the lengths M/π can be considerably reduced by investigating the 3 × 3 × 3-pixel configurations of a binary image instead of the 2 × 2 × 2-pixel configurations.

Further approaches of estimating intrinsic volumes are based on surface rendering [START_REF] Mullikin | Surface estimation of digital planes[END_REF], Lohmann, 1998, Klette et al., 1998, Lindblad and Nyström, 2002, Lindblad, 2003, Windreich et al., 2003, Klette and Rosenfeld, 2004] and the extended versions of the Steiner formula [START_REF] Schmidt | Joint estimators for the specific intrinsic volumes of stationary random sets[END_REF], Klenk et al., 2006, Guderlei et al., 2007, Mrkvička and Rataj, 2008].

Basic Characteristics of Random Structures

The solid constituent of partially open foams or fiber fleeces can be modeled as macroscopic homogeneous random sets. More precisely, the solid constituent of a porous media is a realization of a macroscopically homogeneous random set. For this reason random sets are briefly introduced. Let C and F denote the systems of all compact and closed subsets of R 3 , respectively. Furthermore, let B be the smallest σ-algebra of F containing all sets

F C = {F ∈ F : F ∩ C = ∅}, C ∈ C.
The σ-algebra B forms a Borel σ-algebra with respect to a suitable topology on F [START_REF] Schneider | Stochastic and Integral Geometry[END_REF], Ohser and Schladitz, 2009, Chiu et al., 2013]. Now, a random closed set Ξ can be defined as a measurable mapping from a probability space [Ω,

A , P] into [F, B], which generates a distribution P on [F, B] with P (B) = P(Ξ ∈ B) for every B ∈ B.
Two random closed sets Ξ and Ψ having the same distribution are said to be stochastically equivalent (Ξ ∼ Ψ) and the random closed sets Ξ and Ψ are said to be stochastically independent if

P(Ξ ∈ B 1 , Ψ ∈ B 2 ) = P(Ξ ∈ B 1 ) • P(Ψ ∈ B 2 ), B 1 , B 2 ∈ B.
Furthermore, a random closed set Ξ is called macroscopically homogeneous if its distribution is invariant under translations in R 3 , i.e. all translations Ξ + x of Ξ have the same distribution,

Ξ + x ∼ Ξ, x ∈ R 3 .
Notice that in the literature "macroscopic homogeneity" is also known as "stationarity" or "spatial stationarity". Finally, a macroscopically homogeneous random closed set Ξ is said to be isotropic if the distribution of Ξ is invariant under all rotations in R 3 . Let SO(R 3 ) denote the rotation group (proper orthogonal group) of R 3 , then isotropy of Ξ means that θΞ ∼ Ξ for all θ ∈ SO(R 3 ).

From the Choquet theorem, it follows that a random set Ξ is uniquely characterized by the Choquet capacity T of infinite order, [Schneider and Weil, 2008, mat, ]. For example, macroscopic homogeneity and isotropy of Ξ can respectively be expressed as

T (C) = P(Ξ ∩ C = ∅), C ∈ C,
T (C -x) = T (C), x ∈ R 3 and T (θ -1 C) = T (C), θ ∈ SO(R 3 ),
for all C ∈ C. However, from the image analysis point of view, the capacity functional T is an extremely complicate quantity which cannot be experimentally determined from realizations of Ξ. For this reason, simpler quantities are considered, which, of course, can only characterize different aspects of the distribution of Ξ. The most simple but obvious quantities for characterizing a macroscopically homogeneous random set Ξ are the volume fraction V V , the specific surface area S V , the density of the integral of mean curvature M V , and the density of the integral of the total curvature K V . These quantities are up to multiplicative constants the densities of the intrinsic volumes.

Let W ⊂ R 3 be a compact and convex window with non-empty interior, then Ξ ∩ W is an observation of Ξ in W and the volume fraction of Ξ can be defined as

V V = EV (Ξ ∩ W ) V (W ) ,
where E denotes the expectation with respect to the probability P of Ξ. The densities S V , M V and K V can be defined analogously, but edge effects occurring from the intersection Ξ ∩ W cannot be neglected. For this reason, instead of the window W the window aW with V (aW ) = a 3 V (W ) for the scaling factor a → ∞ is considered. This means that the densities S V , M V and K V are defined as

S V = lim a→∞ ES(Ξ ∩ aW ) V (aW ) , M V = lim a→∞ EM (Ξ ∩ aW ) V (aW ) , K V = lim a→∞ EK(Ξ ∩ aW ) V (aW ) ,
respectively [START_REF] Schneider | Stochastic and Integral Geometry[END_REF]. Finally, it should be noted that the specific Euler number χ V can be defined in an analogous way, but it might be sufficient to make use of the relationship K V = 4πχ V . The solid matter of the open, the partially open foams and fiber fleeces can be modeled by macroscopically homogeneous random sets.

Characteristics of Porous Media

Besides the intrinsic volume densities introduced in Section 2.1.2, there are lots of further characteristics of macroscopically homogeneous random sets known in stochastic geometry and its applications. In the following, characteristics, which are useful to describe the geometry of solid matter or the pore space of porous media are introduced. Some of them are closely related to the intrinsic volume densities.

General Characteristics

The spherical contact distribution function F (r) of a macroscopically homogeneous random set Ξ with volume fraction V V (Ξ) < 1 is defined as

F (r) = 1 - V V ((Ξ c ⊖ B r ) V V (Ξ c ) , r ≥ 0 (2.1)
[ [START_REF] Chiu | Stochastic Geometry and Its Applications, 3rd Edition[END_REF]. The name "contact" distribution comes from the fact that F (r) is the distribution of the radius of the smallest sphere hitting the topological closing Ξ c of Ξ c under the condition that the origin belongs to Ξ. Notice that in the literature the spherical contact distribution of the shortest distance dist (x, Ξ) assumes that x belongs to Ξ c , see e.g. [START_REF] Stoyan | Stochastic Geometry and Its Applications[END_REF].

The spherical granulometry distribution function G(r) of a macroscopically homogeneous random set Ξ is defined as

G(r) = 1 - V V (Ξ • B r ) V V (Ξ) , r ≥ 0, (2.
2) [Matheron, 1975]. Let Ξ be a macroscopically homogeneous system of non-overlapping spheres, then the granulometry distribution G(r) is the volume weighted distribution of sphere radii. The interpretations of the spherical contact distribution and the spherical granulometry with respect to their application in the characterization of the pore space of porous media are very similar. Both distributions can tell something about the "size" of the pore space. The main differences are as follows: If the solid matter can be modeled as a Boolean scheme, then there exist explicit formulas for contact distribution functions [START_REF] Schneider | Stochastic and Integral Geometry[END_REF]. This is not the case for granulometry distribution functions.

On the other hand, neither the solid matter of the partially open foam nor the fleece investigated in this thesis can be modeled as Boolean schemes. The use of the spherical granulometry distribution has the clear advantage compared to the spherical contact distribution that granulometry can be seen as a "generalized concept of size" [START_REF] Haas | Morphologie mathématique et granulométries en place[END_REF]. Insofar, the spherical granulometry distribution is very useful to characterize the "size distribution" of the pore space in partially open foams and fleeces.

The percolation probability describes the connectivity of pores of foams or the connectivity of cavities of fiber fleeces. In order to explain what is meant by percolation probability, a thick slice of a porous medium is considered, i.e. the part of the porous medium between two parallel planes f distance d > 0 and normal direction ϑ. In the following ρ(d, r) is the volume fraction of pores through that a particle of radius r can percolate from the first to the second plane, where the normal direction is fixed as the z-direction (i.e. the principal flow direction). The quantity ρ(d, r) is defined as the ratio

ρ(d, r) = V V (Ξ c p ) V V (Ξ c ) , (2.3)
where Ξ c p are the pores of Ξ c ⊖ B r hitting the first plane as well as the second one. For details of boundary treatment see [START_REF] Ohser | Estimation of the probability of finite percolation in porous microstructures from tomographic images[END_REF].

Instead of the percolation probability, one could also consider the permeability of the porous media which could be computed directly from 3D images [START_REF] Willot | The Permeability of Boolean Sets of Cylinders[END_REF]. The permeability would be much closer to the physics of random flow. However, the focus of the present thesis is on geometric characterization and, therefore, here the percolation probability is the preferred quantity. An open foam as shown in Fig. 2.1 is a macroscopically homogeneous random structure with the volume fraction V V of the solid matter and the porosity 1 -V V . Let B r + x be a ball of radius r > 0 with the center at the point x ∈ R 3 . Furthermore, let S be a finite set of balls, then the Laguerre cell C i of the ball B r i + x i of S is

Characteristics of Open Foams

C i = {y ∈ R 3 : x i -y -r i ≤ x -y -r, B(r) + x ∈ S}.
(2.4)

The Laguerre diagram T(S) is the set of all Laguerre cells corresponding to the balls belonging to S. A random Laguerre tessellation is a Laguerre diagram with respect to a random system of balls, see [Aurenhammer, 1987] for further details. Finally, open foams can be modeled based on the edge system of a random Laguerre tessellation [START_REF] Lautensack ; Lautensack | Random Laguerre Tessellations[END_REF], Redenbach, 2009]. More precisely, a parallel set of the edge system of a realization of a macroscopically homogeneous random Laguerre tessellation is a proper model of a partially open foam, where the dilated edges are said to be the struts of the foam.

For chromatography, the surface area or the surface density S V is the most important characteristic. A larger surface area offers more binding capabilities and is required for an efficient separation process. The other densities M V and K V have a specific meaning for open foams. For example,

M V π(1 -V V ) (2.5)
is an appropriate estimator of the specific strut length L V . The accuracy of the estimation increases with a low boundary length of the strut cross-sections. Furthermore, from the network formula of [START_REF] Mecke | The specific connectivity number of random networks[END_REF], one obtains

χ V = K V 4π .
This is an estimate of the mean number of nodes per unit volume. That can be assumed from the fact, that in open foams four struts meet in a node. For details see [Lautensack, 2007a].

An open foam is produced by removing the faces (or parts of the faces) of the cells of a closed foam. By techniques of image processing, it is possible to reconstruct these faces. The reconstruction result is equivalent to the segmentation of the cells from the strut system [Ohser et al., 2009]. Typical characteristics of the cells are the mean cell volume V , the cell surface area S, the cell width b and the cell number. For foams the mean chord length l plays an important role. It is given by

l = V V N L ,
where N L is the mean number of chords per unit length. For a macroscopically homogeneous system of spherical particles with constant diameter d, one obtains In Fig. 2.2a the relationship between the cell width 2r, the chord length l and the PPI value is shown. If the cell volume increases then the circle radius r and the chord length l will increase, too. The principal curvature 1/r will decrease. With increasing chord length, there are fewer pores per inch, so the PPI-value is decreasing.

l = 2 3 d. (2.6) In an isotropic case N L = S V /4 is given in m -1 . However, if N L is given in inch -1 , then N L is
Even if these characteristics can be calculated with the intrinsic volumes as explained in Section 2.1.1, their determination yields to underestimation because small inclusions or burrs caused by the manufacturing process have a significant impact on the result. These disturbances can be excluded with a segmentation process. The segmentation process consists of three steps. In the first step as segmentation techniques watershed is used to separate cells. In the second step the characteristics for each cell are calculated and in the last step, the statistic values over all cells are computed.

The mean cell width is the integration of the shortest distance b(X) between two parallel tangential planes, over all solid angles. With the cell volume, surface and width the roundness ̺ of the cells can be defined. The roundness ̺ is equivalent to the shape factors, which are already introduced in Section 2.1.1. In opposite to [Leichtweiß, 1980], there are only three independent definitions for the roundness ̺ 1 . . . ̺ 3 for 3D objects. The visualized sub-volume consists of 512 3 out of 2048 3 pixels with a size of 0.65 ➭m. The data were acquired by means of phase-contrast ➭CT using synchrotron radiation.

Characteristics of Random Fiber Fleeces

As illustrated in Fig. 2.3, fiber systems are non-woven fabrics made from continuous fibers or of fibers with limited length.

Let Φ be a random system of rectifiable space curves as defined in [Mecke andNagel, 1980, Nagel, 1983]. In this setting, a random fiber fleece is the parallel set Ξ = Φ ⊕ B r , where r > 0 is the fiber radius, and Φ forms the random system of fiber cores. In order to ensure non-overlapping and smoothness of the fibers, it is supposed that there is a ε > 0 such that the fiber system Ξ is morphological closed with respect to a ball B ε of radius ε, [START_REF] Wirjadi | Estimating fibre direction distributions of reinforced composites from tomographic images[END_REF]. Realizations of macroscopic homogeneous random fiber fleece Ξ can be seen as a geometric model of fiber fleeces with circular cross sections of the fibers and a constant fiber diameter 2r.

(Ξ • B ε ) = Ξ,
The distribution of the fiber fleece is invariant with respect to rotations around the z-axis. Besides porosity 1 -V V and surface density S V , the specific fiber length and the direction distribution are important characteristics. The specific fiber length can be estimated with the integral of mean curvature formulated in Eq. (2.5) divided by π. Interestingly, this also holds for curved fibers and for fibers with an arbitrary cross section.

There are numerous approaches to estimate the direction distribution of random fiber fleeces. One of the first estimation approaches deals with the cross section of the fibers produced by planar sections through the material [START_REF] Zhu | Determination of Non-Symmetric 3-D Fiber-Orientation Distribution and Average Fiber Length in Short-Fiber Composites[END_REF]. Obviously, the intersection of a plane with a cylinder produces an ellipse apart from boundary effects. The aspect ratio of the ellipse depends on the angle between the cylinder and the intersection plane. However, this aspect ratio does not uniquely identify one angle. In the paper [START_REF] Wirjadi | Applications of anisotropic image filters for computing 2D and 3D-fiber orientations[END_REF] a solution based on anisotropic and prolate Gaussian filtering is presented. The main idea is that the filter responds to a convolution with anisotropic Gaussian filters gets maximal if the direction of the longest principal axis of the filter is aligned to the fiber direction. The drawback of this approach is the dependency between accuracy and number of considered directions, which is limited by the runtime. A very elegant and faster way to estimate the direction distribution is based on the directed distance transform. The directed distance transform calculates the distance of each foreground pixel to the background regarding all spatial directions given by the maximal pixel adjacency. This distance information is used to get the chord length for each direction in the fiber, which is assigned to each pixel of the chord [START_REF] Altendorf | 3D directional mathematical morphology for analysis of fiber orientations[END_REF]. Now, the inertia moments and axes of the chords can be determined and the eigenvalues of the inertia matrices can be calculated. The eigenvector having the lowest eigenvalue is indicating the fiber direction [START_REF] Altendorf | 3D directional mathematical morphology for analysis of fiber orientations[END_REF]. Nevertheless, the computation effort for 3D images remains high. Shorter computation time can be reached when estimating the local fiber direction as the eigenvector belonging to the longest eigenvalue of the local Hessian matrix of the smoothed image f * g σ [START_REF] Tankyevych | Direction-adaptive grey-level morphology. application to 3D vascular brain imaging[END_REF], Wirjadi et al., 2016], where g s σ is a Gaussian kernel of width σ.

A further method of estimating the density function of the fiber direction distribution is the inverse cosine transform of the translative (inner) integral of the Crofton formula either for the integral of the mean curvature [Ohser and Schladitz, 2009] or from counts of fiber sections in planar sections through the 3D image [Kinderlen andPfrang, 2005, Riplinger andSpiess, 2012]. Finally, it has to be remarked that also parametric methods of estimating densities of the fiber direction distribution exist, see e.g. [Ohser and Schladitz, 2009].

Besides the more or less volume based characteristics, the geometry of particle paths plays an important role for the efficiency of chromatographic processes.

Particle Paths

Local curvature and torsion are basic characteristics of cell paths through porous media. These properties are very interesting from the viewpoint of cell chromatography. Cell paths with a given curvature and torsion are acting as microcentrifuges driven by the fluid-solid flow. High curvature and torsion induced by the shape of the pore space are closely related to a high deposition rate of the cells at the inner surface. Thus, investigating the distributions of curvature and torsion of cell paths in a fluid-solid flow through a macroscopically homogeneous porous medium can help to evaluate the applicability of the porous medium as a chromatographic filter.

Curvature and torsion of space curves are investigated in numerous application fields [Blankenburg et al., 2016a]. For example, in materials characterization it is well known for a long time that screw and edge dislocations in crystals fundamentally distinguish themselves by their torsion [Rhines, 1977, Zhao et al., 2012]. Particle paths in turbulent two-phase flows are considered in [START_REF] Braun | Geometry of particle paths in turbulent flows[END_REF], Xu, 2007], where the curvature distribution is estimated from Lagrangian paths in fully developed turbulent flows. These investigations are extended on the common distribution of curvature and torsion of particle paths in [Scagliarini, 2011]. It is shown that small scale turbulence features of the flow, as well as the velocity magnitude of the particles, can have a considerable impact on path shapes. Similar studies are made in biology, where, for instance, path torsion distributions are used to quantify the biological activity of freely swimming organisms [START_REF] Crenshaw | Analysis of the three-dimensional trajectories of organisms: estimates of the velocity, curvature and torsion from positional information[END_REF]. Continuous space curves have to be introduced firstly to explain the meaning of curvature and torsion in differential geometry.

Characteristics of Continuous Space Curves

A space curve C may be given by

C = {f (t) : t 0 ≤ t ≤ t 1 }
where f (t) is a rectifiable parametric function, f : R → R 3 . For the sake of simplicity, the parameter t is omitted, if it is not necessary. Assuming that the derivatives ḟ and f are continuous, the local arc length ς(t) of the curve C exist

ς(t) = t t 0 ḟ dt. (2.7)
Then L (C) = ς(t 1 ) is the curve length and ς(t) = ḟ (t) is the length density. Moreover, the curvature κ and the torsion τ can be introduced

κ = ḟ × f ḟ 3 , (2.8) τ = ( ḟ × f ) • ... f ḟ × f 2 ,
(2.9) [Kehtarnavaz and de Figueiredo, 1988], where × denotes the cross product and where ( ḟ × f )• ... f is the triple product of the first, second and third derivatives. The torsion only exists for κ > 0 and when the third derivative is continuous. Notice that, in contrast to the curvature, the torsion is a signed quantity. For a curve C, the integral of curvature K and the integral of torsion T can be expressed as

K (C) = t 1 t 0 κ(t) ḟ (t) dt and T (C) = t 1 t 0 τ (t) ḟ (t) dt.
The curve length L (C), the integral of curvature K (C) and the integral of torsion T (C) are basic characteristics of C.

The unit of κ and τ is m -1 , and thus, their integrals K resp. T are unscaled, i. e.

K (aC) = K (C), T (aC) = T (C)
for all non-zero a ∈ R. Furthermore, K and T are invariant with respect to translations and rotations of C. Finally, it can be noticed that the reflection of C on a plane changes the sign of T (C).

In many papers the curve is expressed by arc length parameterization [Kehtarnavaz and de Figueiredo, 1988, Kim et al., 2013, Lewinger et al., 2005]. The arc length parameterization is useful to illustrate the meaning of the derivations and the relation of curvature and torsion. It is assumed that the inverse ς -1 of the local arc length u = ς(t) exists ς -1 (u) = t(u), then one get the length parameterized function g(u)

g(u) = f t(u) , u ∈ [0, L (C)].
It follows g ′ (u) = 1. Now it is possible to introduce the tangent unit vector T , the normal unit vector N and the binormal unit vector B, given by

T = g ′ (u), (2.10
)

N = T ′ T ′ , (2.11) B = T × N, (2.12)
where T is perpendicular to N , and B is the normal to the osculating plane containing T and N . These formulas can be summarized to the system of Frenet-Serret equations

T ′ = κN, N ′ = -κT + τ B, B ′ = -τ N,      37
2.3. Particle Paths [START_REF] Frenet | Sur les courbes à double courbure[END_REF] and which becomes in matrix notation

   T ′ N ′ B ′    =    0 κ 0 -κ ′ 0 -τ 0 τ 0       T N B    ,
where

κ = g ′′ = g ′ × g ′′ , (2.13) τ = (g ′ × g ′′ )g ′′′ κ 2 , (2.14)
are the curvature and torsion of the arc length parameterized space curve [Spivak, 1979].

Since the length of g ′ is one, the length of vector g ′′ and the area g ′ ×g ′′ of a parallelogram spanned by g ′ and g ′′ have the same value. From Eq. (2.14) it is obvious that the torsion is only well defined, for κ > 0. With Eq. (2.13) the curvature can be seen as the absolute rate of the direction change g ′′ (u) of the unit tangent vector g ′ (u). The torsion depends on the first derivation of the binormal vector, so it can be thought of as measuring the deviation of g(u) from planarity [START_REF] Klette | Digital Geometry[END_REF]. Notice that the curvature and torsion of g on position u is equal to the curvature and torsion of f on position t(u).

Rewriting the derivatives of g(u), one can see that Eqs. (2.8) and (2.13) are equivalent. The same is true for Eqs. (2.9) and (2.14). It follows from the rule for the derivative of an inverse function that t ′ (u) = 1/ ς t(u) and the first derivation yields

g ′ (u) = ḟ t(u) t ′ (u) = ḟ t(u) ς(t(u)) = ḟ t(u) ḟ (t(u)) .
Now, by means of the substitutions one gets the first, second and third derivatives of function g(u),

g ′ (u) = ḟ ḟ , (2.15) g ′′ (u) = f -g ′ (u)ς ḟ 2 = f ḟ -ḟ ς ḟ 3 , (2.16) g ′′′ (u) = ... f -3g ′′ (u)ς ς -g ′ (u) ... ς ḟ 3 . (2.17)
Since ḟ × ḟ is zero, one immediately gets (2.8) from the cross product of the right-hand sides of (2.15) and (2.16). Analogously, (2.9) follows from (2.17).

In image processing, there exists the concept of tortuosity. The terms "torsion" and "tortuosity" are often used synonymously. However, the tortuosity fundamentally differs from torsion as introduced above. The tortuosity is defined as the ratio of the Euclidean distance f (t 1 )f (t 0 ) between the end points of a curve and the curve length L (C)

L (C) -f (t 1 ) -f (t 0 ) L (C) , (2.18)
where the curve length is equivalent to the geodesic distance. The geodesic distance is the length of the shortest path through the pore space connecting the points f (t 0 ) and f (t 1 ) belonging to the pore space [Soille, 1998].

The paper of [START_REF] Patasius | Model based investigation of retinal vessel tortuosity as a function of blood pressure: preliminary results[END_REF] gives a survey on various variants of the definition of tortuosity known from the literature. In particular, the tortuosity depends on how the geodesic distance is measured [START_REF] Grisan | A novel method for the automatic evaluation of retinal vessel tortuosity[END_REF]. This problem is addressed in detail in [START_REF] Chen-Wiegart | Tortuosity characterization of 3D microstructure at nano-scale for energy storage and conversion materials[END_REF], where the context of application is porous LiCoO 2 cathode electrodes of lithium ion batteries. A similar approach is used in [START_REF] Gommes | Practical Methods for Measuring the Tortuosity of Porous Materials from Binary or Gray-Tone Tomographic Reconstructions[END_REF] for the characterization of porous ethylene vinyl acetate. In [START_REF] Peyrega | Estimation of tortuosity and reconstruction of geodesic paths in 3D[END_REF] the tortuosity is estimated from ➭CT images of Thermisorel structure (a fiber material applied for acoustic absorption), where special attention is paid to the reconstruction of geodesic paths. It should be also referred to the article [START_REF] Gaiselmann | Extraction of curved fibers from 3D data[END_REF] where the tortuosity of the fibers of a non-woven material is computed.

One difference between torsion and tortuosity is the value range. The tortuosity takes values between 0 and 1, whereas the torsion can take nearly any value. Moreover, the torsion can remain zero even if the tortuosity reaches the maximum. However, with the curvature and torsion information for each curve point, it is possible to calculate the tortuosity, but the curvature and torsion cannot be calculated by given tortuosity.

Discretization

Two kinds of discretization are considered having their origins in different ways of observing particle paths. The first discretization is induced by sampling the pore space on a cubic lattice. However, the sampling of a space curve on this lattice would almost certainly be empty. For this reason, a discretization scheme of space curves must be introduced that leads to a sequence of connected pixels [START_REF] Klette | Digital Geometry[END_REF]. This representation fundamentally depends on the adjacency of the connected pixels.

The second discretization is related to observing particle paths in sequences of 3D images obtained by time-resolved ➭CT. In such a sequence, the discrete curve forms a set of discrete particle positions acquired at different times, where the discrete particle positions are acquired by estimating the center of gravity of the particles. These discrete positions depend on the particle speed, on the lateral resolution and on the temporal resolution. The main difference in the two discretizations lies in the fact that for the temporal sampling the discrete points of the 3D curves are disconnected and the distance between two positions are varying, while for the spatial discretization the pixels of the discrete curve are connected and equidistant.

The probably most popular and widespread techniques to compute curvature and torsion of discrete curves are based on splines [Kehtarnavaz and de Figueiredo, 1988, Kim et al., 2013, Sangalli et al., 2009]. These splines allow for a simultaneous smoothing effect of the data and recovering a continuous representation of a discrete curve. Smoothing splines offer the possibility to use the continuous formulas for curvature (2.8) and torsion (2.9).

Spline Approximation

A spline is a piecewise defined polynomial g : R → R of degree w and satisfying specific smoothness conditions at the points. More precisely, the w-1 derivatives of g(x) must be continuous at the connecting points. Let (x i , y i ), i = 0, . . . , n -1 be a set of points with x 0 < . . . < x n-1 and let s i (x) be polynomials of degree w > 0, than g(x) can be written as g(x) = s i (x) for x i ≤ x < x i+1 . The polynomials s i (x) have to be at least of third degree to allow the estimation of the third derivation ... g . The third derivation of these cubic splines may not be continuous at the given points, but higher order polynomials loose the shape preserving property. Furthermore, higher order polynomial often leads to oscillations which can alter the sign, the monotony and the convexity of the data [Späth, 1990].

The cubic smoothing splines are firstly presented by [Reinsch, 1967]. They minimize the functional

p n-1 i=0 y i -g(x i ) δy i 2 + x n-1 x 0 g ′′ (t) 2 dt, (2.19)
where p is a smoothing parameter and δy i > 0 is the standard deviation in y i [ de Boor, 1978]. Parameter p controls the smoothness of the solution. For p → ∞ the spline g(x) purely interpolates the data points without smoothing. For 0 < p < ∞ one obtains a smoothed approximation of the curve defined by the data points and for p = 0 function g becomes constant. For polynomials of third degree

s i (x) = a i + b i (x -x i ) + c i (x -x i ) 2 + d i (x -x i ) 3 , x i ≤ x < x i+1 , (2.20)
with the coefficients a i , b i , c i , d i ∈ R and specific boundary conditions to be formulated for s 0 and s n-1 , the minimization of the functional (2.19) yields to an approximating cubic spline g(x) [Pollock, 1999]. The boundary conditions can be chosen arbitrarily and often there is no motivation for the choice. Thus, the so called natural boundary conditions, i.e. s0 (x 0 ) = 0 and sn-2 (x n-1 ) = 0 are chosen. Nevertheless, even for natural boundary conditions, the s i (x), k = 1, . . . , n -3 are influenced by specifying these conditions. Thus, in order to be free of boundary effects one often assumes periodicity, i.e. for y i = y i mod n for i ∈ Z. In this case a periodic approximating spline is adapted to the points (x i , y i ), i = 0, . . . , n -1.

Cubic splines have some disadvantages. The third derivation is piecewise constant and higher order splines lose their shape preserving property. Furthermore, in practical applications, the choice of an appropriate smoothing parameter is very challenging. Moreover, the degree of freedom increases with the choice of boundary handling, the choice of parameterization for higher dimensional curves and with the use of none constant spatial dependent smoothing parameters. All these parameters have a significant impact on the approximation and have to be adjusted accordantly regarding the function and data type, the resolution and the expected displacement error.

Previous Work for Torsion Estimation

There are various approaches to estimate local torsion from 3D images, which differ in their strategies for handling the discretization effects. The approaches suggested in [Kehtarnavaz and de Figueiredo, 1988, Kim et al., 2013, Sangalli et al., 2009] are based on smoothing splines fitted to discrete curves, where the torsion is estimated from the splines. In [START_REF] Pao | Bending and twisting of an in vivo coronary artery at a bifurcation[END_REF] and [START_REF] Crenshaw | Analysis of the three-dimensional trajectories of organisms: estimates of the velocity, curvature and torsion from positional information[END_REF], cubic splines are fitted to the positions of anatomic landmarks and organisms, respectively. As mentioned in the previous section, the used cubic splines have the disadvantage of a piecewise constant third derivation. Furthermore, the boundary condition and the smoothing parameter must be accordantly adapted to the kind of curve and the lateral resolution.

A concept similar to splines is used in [START_REF] Lewinger | Curvature and torsion estimators based on parametric curve fitting[END_REF], Worz et al., 2010], where models of continuous curves are fitted to discrete ones by least-square methods, and the continuous curve representations are again used for torsion assessment. Fitting the points with continuous curves to attenuate the effects of discretization is only efficient with precise prior knowledge (required type of function, degree of the functions, etc.). This knowledge is not available for many applications and in particular for the task of chromatographic filter evaluation.

[ Mokhtarian, 1997] follows a multiscale description of curve shape and applies anisotropic Gaussian smoothing to attenuate discretization effects on torsion estimations. It is obvious to speed up smoothing by Fourier techniques, which leads to fast torsion estimation [START_REF] Medina | Curvature and torsion estimation for coronary-artery motion analysis[END_REF]. This approach is very skillful since the derivatives of the space curve necessary for the calculation of the torsion can also be computed via the inverse space (simultaneously with the smoothing). However, the choice of the appropriate smoothing parameter depends on the resolution in relation to the size of the underlying continuous curve. In the case of the geometric description of arteries [START_REF] Medina | Curvature and torsion estimation for coronary-artery motion analysis[END_REF] assumptions of the kind of curve are made, but in the case of curvature and torsion estimation of particle paths through chromatographic filters investigated in this thesis, such assumptions are not possible.

The core of the algorithmic approach suggested by [Nguyen andDebled-Rennesson, 2008, Nguyen andDebled-Bennesson, 2009] is a decomposition of discrete curves into socalled "blurred" segments, which is fundamentally the same as smoothing. This approach 2.4. Scientific Challenge has the advantage that multidimensional curves can be decomposed into 2D curves, which makes it easily extendable. Nevertheless, the choice of the right step width, which correspond to a smoothing parameter is more or less intuitively and the suggested approach has a relatively high complexity of O(n log 2 n).

To sum up, all approaches need a smoothing parameter for an accurate torsion estimation. This smoothing parameter needs to be accordingly adapted to the curve, which is only possible with precise prior knowledge about the curve shape and resolution of the data. In the case of curvature and torsion estimation of particle paths through chromatographic filters, such knowledge is not available. Therefore, objective criteria are needed to select the smoothing parameter or the estimation approach needs to be widely independent of the smoothing parameter. Finally, some of the approaches only work with curves of connected points, whereas traced cells motion tracked from sequences of 3D images offer only discrete points for low time resolution.

Scientific Challenge

Most of the quantities mentioned above can be estimated from 3D images with moderate statistical and systematic errors as long as the image is of sufficient size and lateral resolution. However, the estimation of torsion is more challenging, since it involves the third derivation of the parametric function. Discretization, as well as image noise, destroy the information on torsion, which must be retrieved with smoothing. However, torsion of particle paths is an important characteristic of porous media regarding their suitability for cell chromatography. In the same way, curvature is an important characteristic of particle paths, although the torsion is likely to have a greater influence on the chromatic separation. In addition, the torsion calculation includes the first three derivatives, whereas the calculation of the curvature requires only the first two. Therefore, it can be expected that an increase in the accuracy of torsion estimation also improves the curvature estimation. That is the reason why in some cases in this thesis only the torsion estimation will be considered.

The aim of this work is to estimate local curvature and torsion as accurate as possible. In this thesis, two new accurate torsion estimation approaches are developed for two discretization schemes of continuous particle paths. Both approaches are widely independent of a smoothing parameter.

Chapter 3

Curvature and Torsion of Particle Paths

In this chapter, two new approaches for torsion estimation based on discrete curves are introduced, where each approach is closely related to a discretization scheme for the curves. The first scheme is based on the extraction of paths of fast particles from 3D

images obtained by ➭CT. The result is a discrete curve of connected pixels which can be seen as an outer Jordan discretization (digitization) [START_REF] Klette | Digital Geometry[END_REF] of a space curve. The second one is directly obtained from the observation of the positions of slow particles at discrete time steps, where moving particles are observed in sequences of 3D images. It will be pointed out that for the first discretization scheme the Fourier approximation of the discrete curves appears as an appropriate method for estimating local curvature and torsion. The second discretization scheme is based on the approximation of the derivatives occurring in the differential geometric formulas by difference quotients. In principle, it is possible to apply the latter approach to the outer Jordan discretization, too, but this would require an additional effort to distinguish discretization effects and image noise from high torsion values of strongly twisting paths.

Conversely, applying the Fourier approximation on sequences of particle positions would require a generalization of the discrete Fourier transform to inhomogeneous lattices (see non-uniform discrete Fourier transform [START_REF] Bagchi | The Nonuniform Discrete Fourier Transform[END_REF]). However, this would have incalculable consequences for the approximation. In particular, it is not yet clear how the Fejér's theorem could be extended to non-uniform sampling. In Section 3.1 the two discretization schemes are described in detail. The extraction of the paths from the 3D images is presented in Section 3.2. The Fourier approximation of discrete particle paths is presented in Section 3.3, whereas discretizations of the differential geometric formulas are given in Section 3.4. Finally, a simulation study is presented in Section 3.5 to give a glimpse on how blurring can affect the computations of derivatives. In particular, the Fourier approximation and the discretization of the differential geometric formulas are compared with spline approximation, which is widely used for torsion estimation [Kehtarnavaz and de Figueiredo, 1988, Kim et al., 2013, Sangalli et al., 2009].

The description of the Fourier approximation on curvature and torsion estimation of particle paths was published in [Blankenburg et al., 2016a]. The discretizations of the differential geometric formulas were firstly presented in [START_REF] Blankenburg | Torsion of particle trajectories through pore space and its estimation using information on local pixel configurations[END_REF].

Discretization Schemes

The discretization of particle paths is induced by the kind of observing particle motion. In the first case, particle motion is simulated based on 3D images of the porous media. Then, the discretization of particle paths leads to a discrete curve, which can be seen as an outer Jordan discretization of a space curve. Under certain conditions, in particular, if the curve is smooth enough, then its Jordan discretization forms a skeleton line which is connected with respect to the 6-adjacency of the pixels (also known as the pixel neighborhood). As far as known, there is no discretization scheme for space curves that leads to unique skeleton lines with respect to other adjacencies, e.g. analogously to the extension of Breshams line on the 3D case. In general, it will be impossible to link the positions of the pixel sequence with moving particle positions at points in time.

In the second case, particle motion is observed in a sequence of 3D images. Assuming that the particle position in each of these images is estimated as the pixel position closest to the center of gravity of the particles, the corresponding sequence of pixels does not form a connected discrete curve with respect to a chosen adjacency of pixels. On the other hand, these pixels represent estimates of positions of moving particles observed at discrete points in time with constant step width on the time scale.

Outer Jordan Discretization

Let Z 3 be the homogeneous 3D lattice of the lattice spacing 1 (which is identical with the pixel size), where Z is the set of integers and [0, 1) 3 is the half-open unit cell. The sampling C of C on the lattice Z 3 is the set of all pixels x ∈ Z 3 with the property that C hits the shifted cell [0, 1

) 3 + x, C = x ∈ Z 3 : C ∩ [0, 1) 3 + x = ∅ .
The union x∈C [0, 1) 3 + x of the lattice cells hitting C is well known as the outer Jordan discretization of C [Jordan, 1892]. Under certain assumptions for the shape of C, the sampling C is a discrete curve with respect to the 6-adjacency of the pixels. More precisely, it is supposed that C may

m T > 0 T < 0 x z y Figure 3
.1: Outer Jordan discretization with m + 3 pixels of a curve piece with a tangent direction (0, 0, 1). m is an integer representing the length of the configurations. As represented in this figure, this length is bounded by the centers of the two pixels with the coordinates (1, 0, 0) and (m, 0, 0). The left configuration has a negative integral of torsion and the right pixel configuration a positive integral.

have a shape such that C is a skeleton line with respect to the 6-adjacency; i. e. C is the skeleton of itself. A necessary condition for that is that the curvature of C must be bounded.

Examples of outer Jordan discretization of curve pieces are shown in Fig. 3.1. Using the right-hand rule one can check whether the integral of torsion T is negative, zero or positive. Such a decision is based on the assumption that the curve pieces behind the discretization are sufficiently smooth with almost everywhere existing local curvature κ.

Nonetheless, fitting curves to discrete ones could be a challenging problem. This is demonstrated in Fig. 3.2, where the outer Jordan discretization of line segments with increasing lengths and direction (1, -1, 1) are shown. Obviously, the straight line has no torsion, but summing up the torsion obtained from the local pixel configuration as shown in Fig. 3.1 would yield values differing from zero. This means one cannot distinguish whether the underlying curve forms a line segment or the path of a twisting particle.

Sequence of Particle Positions

The second discretization scheme covers datasets of 3D image sequences. Although 3D images allow for a torsion assessment which is interesting from the practical point of view, which means a more simple acquisition setup and smaller datasets, the proposed algorithms have also to be evaluated on particle paths extracted from sequences of 3D images. In these images, the positions of moving particles are sampled on a set of discrete temporal positions, i.e. f (t) is discretized with respect to the step ∆t > 0. In this case f (t) can be seen as a four-dimensional parametric function including the time component. In opposite to the first scheme, it is not necessarily an outer Jordan discretization, because the temporal resolution is normally insufficient to get connected positions. As detailed in next paragraphs, there is a localization or displacement error ε k of the segmented particles due to their position representation on the lattice with rounded numbers. Therefore, the noisy position can be expressed with x k = f (k∆t) + ε k , where the random displacement ε k are identically distributed position errors. It is clear that the variance σ 2 of the position error decreases with increasing lateral resolution of 3D images, however, an important question has to be answered: how is σ related to the pixel size α of the image?

To establish a relationship between the standard deviation σ of the pixel position and the pixel size α, it should be assumed that the curve of f (t) and the cubic lattice αZ 3 of the pixels are chosen independently of each other, where Z is the set of integers. More precisely, f (t) may have a random offset uniformly distributed on the unit cell [0, α) 3 . Then the blurred particle positions α⌊f (t)/α⌋ form lattice points of αZ 3 (i.e. pixels of the 3D image). Here, ⌊x⌋ denotes the floor function, i.e. ⌊x⌋ is that point on the unit lattice Z 3 whose components are the integers of the components of x ∈ R 3 . The random displacements ε k = α⌊f (k∆t)/α⌋f (k∆t) are uniformly distributed on the unit cell [0, α) 3 of αZ 3 . As a consequence, each component of the ε k has a standard deviation of σ = α/ √ 12 (pixel accuracy). Alternatively, if the positions x k are estimated as the centers of gravity of the particles, one can reach subpixel accuracy, σ ≪ α/ √ 12, where σ decreases with increasing particle size, see e.g. [van Assen et al., 2002] and references therein.

It should be noted that for the outer Jordan discretization an increased lateral resolution (i.e. reducing the pixel size α) is equivalent with the scaling of C. More precisely, the sampling of C on a lattice αZ 3 with the pixel size α > 0 is the same as the sampling of 1 α f on Z 3 . Similarly, an increase of the temporal resolution of sequences of 3D images would allow for the observation of faster particles.

Simulation of Particle Paths

The simulation of particle paths through porous media is closely related to the simulation of a solid/liquid two-phase flow. In this thesis, the suspension is blood consisting of plasma and blood cells representing the solid phase. Blood is a highly concentrated suspension of diverse blood cells in an aqueous liquid (the blood plasma) and its biological, as well as mechanical behavior, is well known [START_REF] Thiriet | Biology and Mechanics of Blood Flows: Part I: Biology[END_REF], Thiriet et al., 2010]. However, modeling and simulation of the flow dynamics of blood in applications of blood-wash (apheresis) or cardiovascular simulations [START_REF] Formaggia | Cardiovascular Mathematics: Modeling and simulation of the circulatory system[END_REF] are usually restricted to the description of blood flow as an effective, single-phase, non-Newtonian suspension current and its interaction with the deformable wires and effective deposition models e.g. with respect to arteriosclerosis. The formulation of blood rheology depending on the local concentration of the blood cells and its change induced by the selective deposition yields to a multi-phase system according to the multi-phase mixture theory [START_REF] Ehlers | Porous Media -Theory, Experiments and Numerical Applications[END_REF]]. The resulting system of partial differential equations of coupled advectiondiffusion equations with nonlinear Navier-Stokes equations for the blood plasma can be discretized and solved with the usual numerical methods. On the other hand, there are various papers which consider the dynamics of single viscoelastically deformable particles (describing the red blood cells) in a flow. Specific computation methods for the fluid-structure interaction of several particles are e.g. immersed boundary methods or embedded lattice Boltzmann methods, where a current overview in particular on red blood cells is given in [START_REF] Ju | A review of numerical methods for red blood cell flow simulation[END_REF]. These methods are often computationally expensive, so that not yet efficient scale-connecting approaches have been established, although such approaches are currently discussed in literature [Fedosov, 2010, Pan et al., 2011]. The modeling of deposition and separation of cells at inner surfaces of porous media, which is of high importance for the cell chromatography, is only inadequately treated in these approaches. However, the reactive aggregation models for the description of cells clumping according to a Ligand-receptor kinetics [START_REF] Bagchi | Computational fluid dynamic simulation of aggregation of deformable cells in a shear flow[END_REF], Fedosov, 2010] can be adapted for the selective deposition. Significant parameters are the effective minimum distance between the cells and the inner surface because the specific electrostatic potential significantly controls this interaction. This shows that the local curvature of the inner surface influences the local deposition rate. Furthermore, current results on classical particle filter simulations are important for modeling and simulation of cell chromatography, since analytical, as well as numerical multi-scale methods, are used in industrial applications (in addition to various model approaches at different scales). For a review of these methods see [START_REF] Iliev | Numerics[END_REF].

The deposition at the inner surface is described by a fine-scale model, i.e. based on Lagrangian particle physics, which describes the motion of single particles in the flow field and influenced by electrostatic attraction as well as individual interaction with the surface [START_REF] Rief | Computer simulation of filtration including electric surface changes in 3-dimensional fibrous microstructures[END_REF].

Nevertheless, the simulation of blood flowing through porous media remains an open problem. In particular, the formulation of the boundary conditions requires a comprehensive investigation of the interaction of the different cells with the chemical activated inner surface. Solving these simulation problems with the necessary investigation are beyond the scope of this thesis. Moreover, the numerical solution of discretized systems of partial differential equations results in very high computational effort especially when applied to large 3D images.

In the present thesis, a simulation based on techniques of image processing is used. The implementation is usually much faster than the numerical solutions and the formulation of the boundary conditions do not need to be as accurate as in the solution of partial differential equations.

Path of Fast Particles

In order to introduce a simulation of particle motion through the pore space X the geodesic distance gdist X is considered as the shortest path in X between two points x, y ∈ X. Let C X be a curve with the continuous parametric function f (t) satisfying f (0) = x, f (1) = y and f (t) ∈ X for t ∈ [0, 1]. Furthermore, let ζ X be the set of all curves C X , then the geodesic distance gdist X (x, y) is the shortest curve in X connecting x and y, i.e.

gdist X (x, y) = inf{L (C X ) : C X ∈ ζ X },
where L (C X ) is the length of the curve introduced in Section 2.3.1. Now, the geodesic dilation is defined as the set {x ∈ X : gdist X (x, Y ) ≤ t} with Y X and a threshold t ≥ 0 [Soille, 1994]. In order to generalize the geodesic distance, the pore space X is replaced by a continuous integrable function g : R 3 → R with the support of X. Moreover, tdist f is the geodesic time, defined as

tdist f (x, y) = inf      C X g f (t) || ḟ (t)||dt : C X ∈ ζ X      .
(3.1)

The corresponding dilation for the function g

(x) is x ∈ X : tdist X (x, Y ) ≤ t, Y X. (3.2)
If the set Y is the intersection of X with the xy-plane, then the principal motion direction of the particles is the z-direction, which is also the principal flow direction of the suspension. However, it is possible to choose the set Y of starting points arbitrary, for Y ⊂ X.

For the simulation of the paths of fast particles the function -g(x) is the Euclidean distance of x to the complementary set X C of the pore space X,

g(x) = -inf{||x -y|| : y ∈ X C }, x ∈ R 3 .
Unfortunately, the direct implementation of the geodesic dilation is extremely timeconsuming. In order to speed up the dilation, alternative techniques are proposed in some papers [Malladi andSethian, 1996, Baerentzen, 2001]. One of them is based on a numerical solution of the eikonal equation, called Fast Marching Method (FMM) [Sethian, 1999]. Let t(x) be the time at position x ∈ X and v(x) is a scalar speed, v : R 3 → R, with v(x) ≥ 0 and x ∈ X. Then the geodesic dilation is the solution of the boundary value problem

||∇t(x)||v(x) = 1, x ∈ X, (3.3) t(x) = q(x), x ∈ ∂X, ( 3.4)
(eikonal boundary value problem), where q(x) is an exit-cost penalty. The function is chosen as q(x) = ∞ for x ∈ ∂X \ Y and q(x) = 0 for x ∈ Y . Finally, the speed v is set to the Euclidean distance,

v(x) = inf{||x -y|| : y ∈ X C }, x ∈ R 3 .
The speed improvement of the FMM is caused by focusing on neighboring pixels and avoiding recomputations. Each pixel can be labeled with three states: unknown, socalled narrow band and frozen [START_REF] Valero-Gomez | The Path to Efficiency: Fast Marching Method for Safer, More Efficient Mobile Robot Trajectories[END_REF]. At the beginning, all pixels except the start points are set to unknown. The starting points are set to frozen. Now, the geodesic distances of the neighbors are computed and assigned to the corresponding pixels. These pixels are marked as the narrow band. The narrow band pixels with the shortest distance will set to frozen and added to the shortest path. Then the distances of neighbors of the new frozen pixels are calculated and then added to the narrow band. This process is repeated until the end pixels are reached or all pixels are frozen [Baerentzen, 2001].

Local Pixel Configurations

The path obtained by the FMM through the pore space can be seen as a discrete curve, whose pixels are connected with respect to a chosen adjacency system. The curve length can be estimated from local pixel configurations using the technique described in [START_REF] Ohser | Estimation of the probability of finite percolation in porous microstructures from tomographic images[END_REF]. This technique is an alternative to those published in [START_REF] Coeurjolly | Discrete curvature based on osculating circles estimation[END_REF]Svensson, 2003] for instance. It is based on a discretization of the Crofton intersection formula for the integral of the mean curvature [Schneider, 1993]. The backbone of the algorithm to estimate the curve length is simply the determination of the number of 2 × 2 × 2-pixel configurations of a binary image and the computation of a weighted sum of these numbers. The simplicity of this method is due to the fact that the integral of the mean curvature is an additive functional.

The idea arises if this method for estimating the curve length from local pixel configuration is also applicable for the curvature and torsion estimation. For the sake of simplicity, the following considerations are focused on the torsion, but similar examinations can be made for the curvature, too.

Up to a translation and a rotation, there are only two classes of pixel configurations contributing to torsion. These configurations yield from the fact that torsion τ = 0 indicates a deviation from a plane, see Section 2.3.1. The possible pixel classes are illustrated in Fig. 3.1. These configurations can have an arbitrary number of m pixels and therefore, it is not possible to estimate the torsion with a local 2 × 2 × 2 or 3 × 3 × 3-pixel configuration. Moreover, the integral of torsion is not additive and, therefore, it cannot be expected that it can be estimated from local knowledge of a binary 3D image. In other words, the integral of torsion cannot be estimated using the sum of local torsions of curve parts representing different connected pieces of the curve. This can be seen from examples of outer Jordan discretizations shown in Fig. 3.2. In these examples, the pixel sequences correspond to a line segment with tangent direction (1, -1, 1), which has zero curvature and undefined torsion, but lead to an erroneous curvature and torsion estimation due to discretization effects. It is important to keep in mind that the robustness of the curvature and torsion estimation depends on the handling of such discretization effects. One possibility to suppress these discretization effects is to distinguish between pixel configurations inducted by discretization and those of significant torsion. A second possibility is to use a low pass filter in frequencies domain realized by the Fourier approximation.

Fourier Approximation

The proposed Fourier approximation for torsion and curvature estimation starts with the fit of a Fourier series to the points of a discrete curve C . To do so, it is assumed that parametric function f : R → R 3 is periodic with the period length a, i.e. the corresponding curve C is cyclic. If f is Lebesgue integrable on (0, a), its Fourier coefficients

c ℓ = 1 a a 0 f (t) e -2πiℓt a dt, ℓ ∈ Z
exist, where c ℓ = (c ℓ,1 , c ℓ,2 , c ℓ,3 ) ∈ C 3 and C is the space of complex numbers. Furthermore,

s m (t) = m ℓ=-m c ℓ e 2πiℓt a , m = 0, 1, . . .
is the mth partial sum of the Fourier series. It can be noticed that the derivatives ṡm , sm , ... s m and the functions ṡm 2 , ṡm × sm , ṡm × sm 2 and ( ṡm × sm ) • ... s m are partial sums of Fourier series, too. This motivates to fit a partial sum to a sampling of a space curve and to estimate the curvature and torsion of the curve from the partial sum.

In order to discuss the problems that can arise in the Fourier approximation, the helix with the radius r and the slope h should be considered. Its nonzero Fourier coefficients are

c -1,1 = c 1,1 = r/2, c -1,2 = ri/2, c 1,2 = -ri/2, c 0,3 = rhπ and c k,3 = rhi/k for k = 0.
The Fourier series is convergent, i.e. s m (t) → f (t) as m → ∞. The Parseval identity yields

1 2π 2π 0 f (t) 2 dt = ∞ k=-∞ c k 2 = r 2 1 + 4π 2 3 h 2 .
However, the periodicity assumption implies a discontinuity of f (t) at t = 0 and, therefore, the first derivative of the Fourier series is divergent. In particular, one gets

lim m→∞ 1 2π 2π 0 ṡm (t) 2 dt = lim m→∞ m k=-m ikc k 2 = ∞.
As a consequence, the length of the space curve given by the Fourier series of the helix with h = 0 does not exist. Analogously, the L 2 -norms of the higher-order derivatives of s m are divergent, and thus, the integral of curvature K and the integral of torsion T of s m diverge, too. From the previous considerations it can be seen, that the higher-order derivatives of s m do not necessarily converge to the higher-order derivatives of f even if s m (t) → f (t) as m → ∞. To ensure the convergence, the Fourier series must be limited either by partial sums or another low pass filter. Thus, instead of the partial sum s m the corresponding Cesàro mean σ m is considered, which is simply the arithmetic mean of the first partial sums,

σ m (t) = 1 m + 1 m ℓ=0 s ℓ (t).
Fejér's theorem states that if the left and right limits f (t + 0 -) and f (t + 0 + ) exist, then

lim m→∞ σ m (t) = f (t + 0 -) + f (t -0 + ) 2
for all t ∈ [0, a) [Zygmund, 1988]. For this reason, the Fejér summability of a Fourier series is even more important than its convergence.

Let now C be a sampling of the curve C of a parametric function f (t) which is assumed to be periodic on [0, a). The sampling C may consist of n pixels with the positions f 0 , . . . , f n-1 . The mth partial sum s m (t) of a Fourier series can be fitted to the sampling C using the approximation c ℓ ≈ fℓ , where ( f0 , . . . , fn-1 ) denotes the discrete Fourier transform of the vector (f 0 , . . . , f n-1 ) with the setting

fℓ = 1 n n-1 k=0 f k e -2πikℓ n , ℓ = 0, . . . , n -1.
With regard to an explicit presentation of the corresponding Cesàro mean, the partial sums s m (t) of the discrete curve C are now given by

s m (t) = f0 + m ℓ=1 fℓ e 2πiℓt a + n-1 ℓ=n-m fℓ e -2πiℓt a (3.5)
for m ≤ n/2 (Fourier approximation). Notice that the pixel positions f k are points on continuous and infinitely differentiable parametric function, i. e. pieces of the curve of s n/2 (t) are interpolating between the f k (Fourier interpolation). Now Cesàro's mean σ m (t) of the s ℓ (t), ℓ = 0, . . . , m, can be written as

σ m (t) = f0 + m ℓ=1 m + 1 -ℓ m + 1 fℓ e 2πiℓt a + n-1 ℓ=n-m ℓ + m -n + 1 m + 1 fℓ e -2πiℓt a .
From Fejér's theorem it follows, that in the case of a continuous parametric function f the σ m (t) converges uniformly to f (t) on [0, a) for increasing lateral resolution (i. e. for decreasing pixel size) and as m → ∞. The function σ m can be seen as a reconstruction of f from C . Furthermore, Cesàro's mean involves a smoothing that suppresses discretization noise, where the power of smoothing can be controlled by the choice of the parameter m. As a consequence, estimates of the local curvature and torsion from σ m (t) usually have much smaller errors than those from the partial sums s m (t). Moreover, even if the curvature of s m (t) does not converge to the true curvature κ(t), the curvature of σ m (t) can pointwisely converge to κ(t) for increasing lateral resolution and as m → ∞.

Analogously, one can find examples of discrete curves with a divergent torsion of s m (t) but with a convergent torsion of σ m (t). The k-th derivatives σ m,k (t) of σ m (t) can be expressed as are the Fourier coefficients of the k-th derivatives s m,k (t) of s m (t). In other words, in Eq. (3.6) the smoothing and the computation of the derivatives is made simultaneously in one step.

σ m,k (t) = m ℓ=1 m + 1 -ℓ m + 1 2πiℓ a k fℓ e 2πiℓt a + n-1 ℓ=n-m ℓ + m -n + 1 m + 1 - 2πiℓ a k fℓ e -
A similar regularization can be obtained by filtering the discrete signal (f 0 , ..., f n-1 ) with a binomial filter of order ν > 1. For even ν, the transfer function of the binomial filter is cos ν (πℓ/n) and, therefore, a smoothing of f can be rewritten by the binomial Fourier approximation

b ν (t) = n-1 ℓ=0 cos ν πℓ n fℓ e 2πiℓt a . (3.7)
Notice, the meaning of parameter ν differs from that of the smoothing parameter m. However, in order to be able to compare the low-pass filters later, the parameter ν is substituted by ν = nm for 0 < m < n.

The transfer functions can be seen in Fig. 3.3. The use of Eqs. (3.5), (3.6) and (3.7), respectively, is equivalent with the application of a low pass filter.

The application of the Fourier approximation to a non-periodic function f (t) sampled on an interval [0, 2π] with a step width ∆t and the sample point number n usually induces a discontinuity of f (t) on the interval bounds by the periodic continuation of f (t) (boundary effects). As a consequence, the spectrum of periodically continued function contains high frequencies with high amplitudes, and it would be hard or impossible to separate these frequencies from that of the noise. To attenuate or avoid these boundary effects one can proceed as follows. The function f (t) can be reflected at the origin, which doubles the number of sample points. Therefore, the smoothing parameter m adjusting the low pass filters must be doubled, too. An alternative boundary treatment is the subtraction of a polynomial. Let g(x) be a polynomial of third degree satisfying

g(0) = 0, g(2π) = f n -f 0 , ġ(0) = f 1 -f 0 ∆t and ġ(2π) = fn-f n-1 ∆t .
Then the Fourier approximation is computed for g k = f k -g(k∆t) and the derivatives of f (t) are estimated by the derivatives of the Fourier approximation plus the corresponding derivatives of g(t). Clearly, if the constrains ġ(0) = f 1 -f 0 ∆t and ġ(2π) = fn-f n-1 ∆t are dropped, then a first degree polynomial (straight line) can be used and the error prone calculation of the first derivative is obsolete.

Finally, the derivatives of f can be inserted in (2.8) and (2.9), to estimate the curvature and torsion, respectively.

The Fourier approximation is well suited for the simulated path with equidistant pixel positions. However, positions obtained from a sequence of 3D images are not equidistant. Therefore, another curvature and torsion estimation approach is needed.

Discretization of Differential Geometric Formulas

The basic idea of the difference method used as the numerical solution of the differential equation is to replace differential quotients by difference quotients in order to solve the discrete equation obtained in this way. In the following, an estimation method for local curvature and torsion of particle paths based on a discretization of the differentialgeometric formulas (2.8) and (2.9) is presented, where the differential quotients occurring in these formulas are approximated by corresponding difference quotients. Furthermore, the estimation method is well adapted to the specific kind of discretization, where either particle paths are sampled on a point lattice or the positions of moving particles are sampled on a set of discrete temporal positions.

Using the forward differences

∆ + f k (t) = f (t + k∆t) -f (t + (k -1)∆t) (3.8)
for k = 1, 2, 3 and a step ∆t > 0, the derivatives of f can be written as

ḟ (t) = lim ∆t→0 f (t + ∆t) -f (t) ∆t = lim ∆t→0 ∆ + f 1 ∆t (3.9) f (t) = lim ∆t→0 f (t + 2 ∆t) -2 f (t + ∆t) + f (t) ∆t 2 = lim ∆t→0 ∆ + f 2 -∆ + f 1 ∆t 2 (3.10) ... f (t) = lim ∆t→0 f (t + 3 ∆t) -3 f (t + 2 ∆t) + 3 f (t + ∆t) -f (t) ∆t 3 = lim ∆t→0 ∆ + f 3 -2∆ + f 2 + ∆ + f 1 ∆t 3 . (3.11)
Combining (3.9) and (3.10) with the continuous curvature formula (2.8) yields to

κ = lim ∆t→0 ∆ + f 1 × (∆ + f 2 -∆ + f 1 ) ∆ + f 1 3 = lim ∆t→0 ∆ + f 1 × ∆ + f 2 ∆ + f 1 3 (3.12)
and with the third derivation (3.11) and the torsion formula (2.9) it follows

τ = lim ∆t→0 (∆ + f 1 × ∆ + f 2 ) • (∆ + f 3 -2∆ + f 2 + ∆ + f 1 ) ∆ + f 1 × (∆ + f 2 -∆ + f 1 ) 2 = lim ∆t→0 (∆ + f 1 × ∆ + f 2 ) • ∆ + f 3 ∆ + f 1 × ∆ + f 2 2 . (3.13)
From the right hand sides of (3.12) and (3.13) one obtains appropriate estimators κ and τ of κ and τ at a curve point. For further improvement the estimators given by

κ = ∆ + f 1 × ∆ + f 2 • ∆ + f 2 × ∆ + f 3 ∆ + f 1 • ∆ + f 2 • ∆ + f 3 , (3.14) τ = (∆ + f 1 × ∆ + f 2 ) • ∆ + f 3 ∆ + f 1 × ∆ + f 2 • ∆ + f 2 × ∆ + f 3 , (3.15)
are introduced, respectively. It has been proven that κ and τ behave better than that of the right-hand sides of (3.14) and (3.15), respectively, i.e. their application usually leads to smaller estimation errors. A justification for that is that κ and τ can be represented as geometric means of estimators. For example, analogously to (3.15) one can show that

τ = lim ∆t→0 (∆ + f 2 × ∆ + f 3 ) • ∆ + f 1 ∆ + f 2 × ∆ + f 3 2 , (3.16)
and hence, τ is the geometric mean of the right hand sides of (3.13) and (3.16). But what is the position x on the curve at which κ and τ are "best" estimators? For symmetric reasons, x = f t + t 0 ∆t is chosen, where t 0 = t + 3 2 ∆t is the center of the discrete values t + k∆t, k = 0, . . . , 3.

An alternative for a symmetric estimation of curvature and torsion is based on symmetrized difference quotients. For this reason the backward differences are introduced,

∆ -f k (t) = f (t -(k -1)∆t) -f (t -k∆t),
(3.17)

for k = 1, 2, 3 and ∆t > 0. The symmetrized derivatives can now be written as

ḟ (t) = lim ∆t→0 f (t + ∆t) -f (t -∆t) 2∆t = lim ∆t→0 ∆ + f 1 + ∆ -f 1 2∆t (3.18) f (t) = lim ∆t→0 f (t + ∆t) -2 f (t) + f (t -∆t) ∆t 2 = lim ∆t→0 ∆ + f 1 -∆ -f 1 ∆t 2 (3.19) ... f (t) = lim ∆t→0 f (t + 2∆t) -3f (t + ∆t) + 3f (t) -f (t -∆t) 2∆t 3 + f (t + ∆t) -3f (t) + 3(f (t -∆t) -f (t -2∆t) 2∆t 3 = lim ∆t→0 f (t + 2 ∆t) -2 f (t + ∆t) + 2 f (t -∆t) -f (t -2 ∆t) 2∆t 3 = lim ∆t→0 ∆ + f 2 + ∆ -f 2 -(∆ + f 1 + ∆ -f 1 ) 2∆t 3 . (3.20)
Notice, that the shifted second derivation f (t + ∆t) in Eq. (3.19) is equal to the second derivation f (t) in Eq. (3.10), but the expression of first derivatives in Eq. (3.18) and (3.9) differ. A similar observation can be made for the third derivatives in (3.20) and (3.11). The above formulas can also be applied for the curvature and torsion estimation. It holds

ḟ (t) × f (t) = lim ∆t→0 (∆ + f 1 + ∆ -f 1 ) × (∆ + f 1 -∆ -f 1 ) 2∆t 3 = lim ∆t→0 ∆ + f 1 × ∆ -f 1 ∆t 3 (3.21) κ(t) = ∆ + f 1 ×∆ -f 1 ∆t 3 ∆ + f 1 +∆ -f 1 2∆t 3 = 8 ∆ + f 1 × ∆ -f 1 ∆ + f 1 + ∆ -f 1 3 (3.22) τ (t) = (∆ + f 1 × ∆ -f 1 ) • (∆ + f 2 + ∆ -f 2 ) -(∆ + f 1 + ∆ -f 1 ) ∆ + f 1 × ∆ -f 1 2 ∆t 3 2 2∆t 3 ∆t 3 = (∆ + f 1 × ∆ -f 1 ) • (∆ + f 2 + ∆ -f 2 ) 2 ∆ + f 1 × ∆ -f 1 2 . (3.23)
Notice that in Eqs. (3.22) and (3.23) an additional point is needed for the torsion estimation.

One of the main issues of the difference method is the choice of an appropriate discrete operator. The minimum requirement is to meet the consistency, where a differential operator is consistently discretized when the obtained discrete operator converges to the continuous one as the step width ∆t → 0. In image processing, this necessary variability of the step width is usually not taken into consideration, but the initial pixel size is used as the finest step width ∆t. The determination of the step width is the real reason why the above concept of consistency is not the best choice for evaluating discrete operators.

This consideration was the motivation of Hanno Scharr [Scharr, 2000] to extend the consistency concept. In this respect the common first order derivation operator written on terms of image processing as the filter mask

1 2∆t
( 1, 0, -1 ), (3.24)

is replaced with the mask

D 1 = 1 12∆t ( -1, 8, 0, -8, 1 ) (3.25)
representing higher order of consistency. Analogously, the discrete Laplace filter

1 ∆t 2 ( 1, -2, 1 ) (3.26)
of length three can be replaced with D 2 = 1 12∆t 2 ( -1, 16, -30, 16, -1 ),

(3.27) [Scharr, 2000].

In the terminology of image processing the third order derivation operator is written as

D 3 = 1 2∆t 3 ( 1, -3, 3, -1, 0 ) + ( 0, 1, -3, 3, -1 ) = 1 2∆t 3 ( 1, -2, 0, 2, -1 ), (3.28)
which is simply the convolution of (3.24) and (3.26). However, is D 3 really the best choice? In fact, following Scharr's formalism, this discrete third order derivation operator D 3 of length five is the solution of the linear equation

       1 1 1 1 1 -2 -1 0 1 2 4 1 0 1 4 -8 -1 0 1 8 16 1 0 1 16        D * 3 =        0 0 0 3! ∆t 3 0        , (3.29) -2 0 2 4 -3 -2 -1 0 1 2 3 Intensity ξ I C2 C3
Figure 3.4: Imaginary part of the transfer functions for a step width ∆t = 1 of the Fourier transform of the continuous third order derivation operator (I), of the filter mask for the third order derivation filter with second order consistency (C2), of the filter mask for the third order derivation filter with third order consistency (C3).

where D * 3 is the reflection of D 3 at the origin. This means, D 3 is a discrete deviation operator of second order consistency.

Clearly, one can consider higher-order, consistency. For example, third order consistency would yield to the third order derivation operator 1 8∆t 3 ( -1, 8, -13, 0, 13, -8, 1 )

(3.30) Fig. 3.4 shows the imaginary parts of transfer functions of the filter masks for the third order derivation filters with second respectively third order consistency compared to the imaginary part of the Fourier transform -ξ 3 of the continuous third order derivation operator. For ∆t = 1 the filter mask D 3 given by (3.28) has the transfer function -i(2 sin ξsin 2ξ).

The corresponding filter mask of third order consistency has the transfer function -i(13 sin ξ -8 sin 2ξ + sin 3ξ) 4 .

The consistency specifies only consistent behavior between an ideal continuous derivation operator and the corresponding discrete operator at the zero point in the Fourier domain (see Fig. 3.4). For an decreasing step width ∆t → 0 the restriction of the condition to the origin does not provide any problems. However, in image processing, the pixel size is the minimum step width. Therefore, increasing the order of the filter consistency decrease the error between the discrete derivation operator and the ideal one until an optimum is reached. Further increase does not make sense.

A possible regularization in the use of discrete derivation operators is the smoothing with an appropriately chosen smoothing filter. For example M = 1 m (1, . . . , 1) be the mask of the mean value filter of odd length m ≥ 3,

D 3 * M 3 = 1 2m∆t 3 ( 1, -1, -1, 0, . . . , 0, m-2 times 1, 1, -1 ),
which can reduce the estimation error of the third derivative, where m serves as the smoothing parameter. Another possibility is the choice of a multiple step width m∆t. For instance, one gets for the double step width 1 4∆t 3 ( 1, 0, -2, 0, 0, 0, 2, 0, -1 ), etc.

Finally, in view of the next chapter f (t) is turned to the length parameterization form, where the constant step ∆t on the parameter space is replaced with a constant step ∆ς along the curve. More precisely, assume a constant distance ∆ς between the curve points x k and x k+1 measured along the curve, k = 0, 1, 2, where ς(t) = t 0 ḟ (t) dt is the arc length of the curve C. Then the corresponding step widths ∆t k fulfilling x k+1 = f (t+∆t k ) are obtained by solving the (usually nonlinear) equation ς(t) + k∆ς = ς(t + ∆t k ) for k = 0, 1, 2.

Estimation of the Third Derivatives of Real Functions -a Case Study

The accuracy of curvature and torsion estimation of the curve of a parametric function depends on the accuracy of the estimated derivatives occurring in the corresponding differential geometric formulas. For this reason, in this section, it will be examined how the choice of an estimation method of the third derivative of a real function affects the accuracy of the estimate. Usually, the errors of estimating derivatives of functions increase with increasing order of the derivatives. Thus, one can expect that the error of estimating torsion is dominated by the errors occurring in the third derivative. Let now f be a real function, f : R → R, with an everywhere existing third derivative. Assume, that f is sampled with a uniform step width ∆t > 0. Now, in order to study the influence of noise (blurring of the function) on the estimates, independent and identically distributed ǫ k are added to the sampled function f (k∆t),

f k = f (k∆t) + ǫ k , k ∈ Z,
where the ǫ k are independent and identically distributed random numbers with expectation zero.

- An example is shown in Fig. 3.5, where f (t) is the cosine function and the ǫ k are Gaussian distributed with zero expectation and a variance σ 2 . More precisely, let f k given by f k = cos ∆t + ǫ k with ∆t = 2π/n, where n is the sampling rate, for n > 1. The random numbers ǫ 0 , . . . , ǫ n-1 are pairwise independent and ǫ k = ǫ k+n such that f k is periodic, with f k = f k+n for k ∈ Z. The error of estimating the third derivative depends on the step width ∆t, the variance σ 2 and the choice of the estimation method. If an approximating cyclic cubic spline with the smoothing parameter p is applied to f 0 , . . . , f n-1 then the estimated third derivative is a piecewise constant monotone function, see Fig. 3.5. The Fourier approximation is a continuous function, independent of the choice of the smoothing parameter m > 1. Clearly, for m < n 2 the Fourier approximation coincide with f (t) since the values of the discrete Fourier transform of f (k∆t) are zero for k = 2, . . . , n -2. Notice that using the discrete differential-geometric formulas, the local curvature and torsion can only be estimated at discrete positions shown as crosses in Fig. 3.5. This is a significant difference to the use of the Fourier and spline approximations.

Notice that the Fourier and the spline approximation involve a smoothing, whereas the discretization of the differential-geometric formulas is used without smoothing, which explains the rough graph of the discretization of the differential-geometric formulas. The Fourier and spline approximation require smoothing parameters which must be specifically adjusted to step size ∆t and the expected displacement error ǫ with variance σ 2 . Nevertheless, the Fourier approximation is even more accurate (closer to the ground truth) than the spline approximation. However, in practice, where the ground truth is unknown it is quite difficult to find the best parameter. This is especially true for the spline smoothing parameter, which can only be found approximately because it is a floating point number. In contrast, the Fourier smoothing parameter is a positive integer number, which in practice rarely exceeds 20. In the following, the influence of smoothing on the estimation error is studied in more detail. Let

δ = 1 n n k=1 ... f (k∆t) - ... f (k∆t) 2 (3.31)
be the mean error of estimation. An overview on the estimation error for different smoothing parameters is given in Figs. 3.6 and 3.7.

Fig. 3.6 shows the impact of different smoothing parameters on the accuracy of the Fourier approximation. Each point in the images of Fig. 3.6 gives the mean error of a cosine function calculated for a specific step width ∆t and for a specific variance σ 2 of the Gaussian distributed position error ε. It is difficult to give a relative error since the integral of the cosine function on [0, 2π] is zero.

Clearly, for the Fourier approximation, as well as for the spline approximation, the mean error increases with increasing σ, see Figs. 3.6 and 3.7, respectively. Furthermore, for the Fourier approximation, the mean error increases for increasing step width ∆t, as well as for decreasing smoothing (increasing m as long as m > 1). The latter is not so clear for the spline approximation. A moderate step width ∆t can contribute to smoothing and, thus, it reduces the error of estimating the third derivative. Analogously, a moderate smoothing parameter p (p ≈ 14) yields the best results, while the spline approximation uniformly tends to zero as p → ∞. From

1 2π 2π 0 cos 2 t dt = 1 2 (3.32)
it follows δ → 1 √ 2 as p → ∞ and σ → 0 for arbitrary step width ∆t. Notice that the spectrum of the test function consists of only one frequency and so the Fourier approximation can be well adapted. Such a test situation is not representative of real data in which the signals correspond to a band of frequencies.

The influence of the Fourier smoothing parameter depends strongly on the used low pass filter. Moreover, the filter significantly changes the ratio between low and high error Notice that the point number is changing with different step width. Therefore, the best results for the binomial filter can be reached with a smoothing parameter depending on the point number n. However, the highest accuracy over the whole parameter space can be reached with the partial sum, so this filter will be used for the further comparisons.

The Cesàro mean has the interesting property that for large and small step widths the influence of the displacement error is suppressed. The diagram of the binomial filter has the same shape like the spline approximation, but the mean error is significantly higher in most cases.

In principle, the discretization of the differential-geometric formulas does not involve any smoothing, but a smoothing (e.g. by applying a linear smoothing filter) can be involved as a pre-or postprocessing step. The effect of different masks on the estimation error is shown in Fig. 3.9. The third-order derivative operator becomes numerically unstable for small step widths because the third derivation is normalized by the division with the term ∆t 3 . However, the accuracy for high step widths can be significantly improved with a higher consistency, but it also increases the mean error for σ > 10 -1 . Notice that the σ-axis is in logarithmic scale. The use of a smoothing filter has the opposite effect. Unfortunately, a higher consistency or smoothing implies the use of at least two additional positions, which is inappropriate in the presented application, because the number of torsion estimates would become very small. Therefore, smoothing is avoided.

The best results are achieved by the Fourier approximation. It has a small error nearly independent of the step width. This can be explained by the fact that the whole function is taken into account during calculation rather than only local information. In opposite, the spline approximation gets only accurate results for high-resolution data and small displacement errors. This last result is a first indication that the spline approximation is effectively not adapted to sequences of 3D images in which the resolution is low and the sampling is done with a nonconstant step.

To sum up this preliminary results on the third derivative indicate that the Fourier approximation is more accurate than the spline approximation when prior knowledge is available and the smoothing parameters can be roughly adjusted to the data, the first results, see Figs. 3.6 and 3.7. Moreover, without prior knowledge or with data including many different curve types or data with only a few number of points, the discretization of the differential-geometric formulas can effectively represent a potential alternative solution.

Conclusion

In this chapter, a particle path simulation based on the geodesic dilation was proposed. In particular, the geodesic dilation is applied on the Euclidean distance transform of the pore space to extract the fastest particle path. It was shown, that the geodesic dilation can be seen as a solution to the eikonal boundary value problem. Therefore, the FMM can be used to increase the computation speed of the geodesic dilation. In this frame, the distance to the inner surface is seen as a scalar speed function v(x) and the flow direction can be specified with the choice of appropriated exit-cost penalty q(x).

Furthermore, the curvature and torsion estimation methods based on the Fourier approximation and the discretization of the differential-geometric formulas were introduced. These approaches were published in [Blankenburg et al., 2016a] and [START_REF] Blankenburg | Torsion estimation of particle paths through porous media observed by in-situ timeresolved microtomography[END_REF], respectively. The Fourier approximation leads to a continuous function, which can be derived several times and provides smaller errors than the spline approximation. In particular, for a weak as well as strong noise, the estimation error of the Fourier approximation is usually smaller than that of the spline approximation. The reason for this is that, by uses the Fourier approximation, the high frequencies belonging to noise can be separated from the low frequencies of a sufficiently smooth function by simple thresholding. However, the Fourier approximation requires periodicity and equidistant discretization steps.

The conditions are not fulfilled by time-resolved particle positions. Therefore, a new parameter free approach based on the discretization of the torsion and curvature formulas was presented. The derivative approximation for higher step widths looks promising if considered that no smoothing is done. However, the new discrete formulation of curvature and torsion are not tested, because the curvature and torsion are only defined for 2D and 3D curves, respectively. The Fourier approximation can be easily extended to higher dimensional space, which is the basis for the curvature and torsion estimation.

Chapter 4

Evaluation of the Estimation Methods

In this part of the thesis, the methods of estimating curvature and torsion presented in the previous chapter are evaluated based on examples of space curves with known curvature and torsion. For the Fourier approximation, the estimation errors result from the periodicity assumption involved into the discrete Fourier transform, as well as from the outer Jordan discretization. These effects will be separately studied in Section 4.2.1 and 4.2.2, respectively, where parts of the results presented in these sections are published in [Blankenburg et al., 2016a]. The estimation method based on replacing the differential quotients by differences quotients in the differential geometric formulas is investigated in Section 4.3 which contains results presented in [START_REF] Blankenburg | Parameter free torsion estimation of curves in 3D images[END_REF]. Finally, in Section 4.4 both estimation methods are compared with those known from the literature.

Space Curves Examples

Two classes of space curve examples with analytically known curvature and torsion are considered. Let

f 1 (t) = r    cos t sin t ct    , t ∈ [0, 2π), r > 0 (4.1)
and

f 2 (t) = r    cos t sin t c sin 2 t    , t ∈ [0, 2π), r > 0, (4.2)
be these two classes of parametric functions with the additional parameters r and c.

Parameter r is a scaling factor, whereas the local torsion can be adapted with c. These classes of space curves are so rich that probably most curves occurring in applications can be piecewise fitted by class members. In particular, for f 2 , the parameters r and c can be chosen such that the local change of the torsion strongly varies. The advantage of f 1 in the frame of the algorithm evaluation is that the local curvature and torsion are constant.

Eq. ( 4.1) describes the class of helices of pitch 2πrc and length L (C) = 2πr √ 1 + c 2 , see Fig. 4.1a for a representation of this class of curves. One can easily verify that the local curvature κ(t) and torsion τ (t) are independent of t. It holds

κ = 1 r(1 + c 2 ) (4.3) and τ = c r(1 + c 2 ) . (4.4)
It is noticeable that both the integral of curvature

K (C) = 2π √ 1 + c 2 (4.5)
and the integral of torsion

T (C) = 2πc √ 1 + c 2 (4.6)
are independent of parameter r.

The curve given by Eq. ( 4.2) and shown in Fig. 4.1b is cyclic since the components of the parametric function are periodic with a period length 2π. The curve length can be determined using the derivative of f 2 given by

ḟ2 (t) = r √ 2 2 + c 2 (1 -cos 4t). (4.7)
The curve length is obtained using Eq. (2.7) in the interval [0, 2π]

L (C) = 2π 0 ḟ2 (t) dt = 2r E(-c 2 ) + 1 + c 2 E c 2 1 + c 2 , (4.8)
where E is the (incomplete) Legendre elliptic integral of the 2nd kind given by which yields the complete elliptic integral of the 2nd kind when φ = π/2:

E(φ, x) = φ 0 1 -x sin 2 ϕ dϕ, (4.9) (a) (b)
E(x) = E( π 2 , x) = π/2 0 1 -x sin 2 ϕ dϕ, x ≤ 1.
The local curvature can be expressed as

κ(t) = 2 r 2 + c 2 (5 + 3 cos 4t) 2 + c 2 (1 -cos 4t) 3 (4.10)
and the integral of curvature is given by the expression (4.11) which can be computed numerically, only. In a similar way, the local torsion can be calculated as

K (C) = √ 2 2π 0 2 + c 2 (5 + 3 cos 4t) 2 + c 2 (1 -cos 4t) dt,
τ (t) = - 6c sin 2t r 2 + c 2 (5 + 3 cos 4t) , (4.12)
and its integration yields (4.13) which converges to -π 2 as c → ∞.

T (C) = π/2 0 τ (t)dt = arctan c -2 arctan 2c,
As detailed in Section 3.3, both the curvature and the torsion can be estimated by inserting the Fourier partial sum s m,k (t), see Eq. (3.5), into the differential geometric formulas of κ and τ (see Eqs. (2.8) and (2.9)) given in Section 2.3.1. In particular, the real curve is approximated by a curve of the m-th first Fourier coefficients computed with the discrete points of C and where k stands for the order of the derivative. In a similar way, the curvature and the torsion can also be estimated by inserting Cesàro's mean σ m,k (t) or the binomial filter b m,k (t) into the differential geometric formulas. This yields the estimates

κ s = ṡ m × s m ṡ m 3 , τ s = ( ṡ m × s m ) • ... s m ṡ m × s m 2 , (4.14)
and (4.15) respectively. In the same way, the curvature κ b and torsion τ b with binomial filter can be calculated.

κ σ = σ m × σ m σ m 3 , τ σ = ( σ m × σ m ) • ... σ m σ m × σ m 2 ,
Both examples of space curves are discretized regarding the outer Jordan discretization introduced in Section 3.1. The discrete curves C are assumed to be samplings of independent and isotropic uniformly random (IUR) rotations of C [START_REF] Coleman | Random rotations in simulation with computer 3D reconstruction[END_REF]. A random offset is added which is uniformly distributed in the lattice cell [0, 1) 3 . As a consequence, curvature and torsion estimates depend only on the function type, the scaling defined by parameters r and c and smoothing parameter m.

Boundary Effects

The Fourier series decompose functions into a sum of cosine and sine functions. Naturally, only periodic continuous functions could be represented with a finite number of summands and without boundary effects. The periodic continuation of an originally nonperiodic function usually induces a sharp step at the bounds of the periodicity interval. This causes high coefficient values at high frequencies in the Fourier series. Consequently, first and higher order derivatives of the partial sums of the Fourier series do not converge to the derivatives of the function. The worst situation in terms of curvature and torsion estimation accuracy is when this boundary effect is associated with noise affecting the discrete point positions. Therefore, the aim of this section is to investigate the impact of open curves on the proposed Fourier approximation.

An extreme example is the periodically continued helix f 1 (t) having a strong discontinuity at t = 0. The results of curvature and torsion estimation are shown in 2650 sampling points which results in the total sample size of about n = 2.7 • 10 6 . Notice, that the continuous t parameter for each sample point is stored and, therefore, the curvature and torsion estimates can be linked with the corresponding t values. Moreover, the mean curvature κ (t) and torsion τ (t) for a particular t-value can be calculated. The mean error

δm = 1 n n i=1 τ m (t i ) -τ (t i ) 2 (4.16)
between torsion estimates τ m (t i ) and the true torsion τ (t i ) quantifies the mean error of estimation, where i = 1, . . . , n is the index of the ith sampling point. Analogously, the mean error can be introduced for the curvature estimation. Notice that the curvature and the torsion have the unit 1/pixel size, but the pixel size in the simulated data has no physical unit. Therefore, the unit is omitted.

The frequencies resulting from the discontinuity are superposed with those caused by the sampling of the curve on the lattice. Clearly, the deviation shown in Fig. 4.2 are dramatically large, but one can also see that in general the deviation of the local curvature and torsion estimated form Cesàro's mean are much smaller than those from the corresponding partial sum. The corresponding mean errors are -3 for the torsion estimation with the partial sum, δ4 = 6.9 • 10 -3 for the torsion estimation with Cesàro's mean.

As shown later, it turns out that the error caused by the sampling is much lower than the one shown here. This means that the estimation errors mainly result from the discontinuity of f at t = 0.

Reflection of f at t = 0 removes the discontinuity. An appropriate reflection condition for the proposed Fourier approximation can be reached by replacing the sampling

(f 1,0 , . . . , f 1,n-1 ) of f 1 (t) with (f 1,0 , . . . , f n-1 , f 1,n-1 , . . . , f 1,0 )
, where now the pixel number is 2n. The smoothing parameter m has to be adapted to the increased pixel number; instead of m the smoothing parameter 2m is used. The effect is shown in Fig. 4.3. The reduction of the mean estimation errors is considerable when boundary values are omitted. The maximum of mean curvature and torsion drop down from 0.015 and 0.016 (see Fig. 4.2) to 0.001 and 0.0015 (see Fig. 4.3), respectively. The errors are diminished by about one order of magnitude when the reflection is used: An alternative method to remove the discontinuity of f 1 (t) is to decompose f into the sum

f (t) = f * 1 (t) + f 1,n-1 -f 1,0 n -1 t, (4.17)
where n is still the point number and f * 1 (t) is continuous, see Section 3.3. Now, the calculation of the derivatives via the Fourier approximation is applied only on the first summand f * 1 (t). One gets

ḟ1 (t) = ḟ * 1 (t) + f 1,n-1 -f 1,0 n -1 , (4.18)
for the first derivation, whereas the higher order derivations of f 1 (t) and f * 1 (t) coincide. The effect of restricting the Fourier approximation on f * 1 (t) is shown in Fig. 4.4 for function f 1 (t) which is not reflected. This restriction leads to a significant reduction of the deviation compared to the handling of the boundary effects by reflecting of f 1 (t) at t = 0. Furthermore, when restricting the Fourier approximation on f * (t), the use of Cesàro's mean causes a systematic overestimation of curvature and torsion. However, the standard deviation is very small, which leads to a smaller error than by using the partial sum. This overestimation can be softened with the selection of a higher smoothing parameter, for instant m = 8 (as seen in Fig. 4.4), which decreases the mean error from δ4 = 0.66 • 10 -3 to δ8 = 0.34 • 10 -3 . However, this solution also results in a higher standard deviation of the curvature and torsion estimation. Unfortunately, the mean value of the Cesàro's mean asymptotically approaches the true torsion value for m → ∞, so the overestimation cannot be avoided completely.

It is possible to use higher order polynomials instead of f * 1 (t) and add further constraints regarding the continuity of higher order derivations at the boundaries. However, the calculation of derivatives at the boundary is very error prone resulting in polynomials with strong slope, which has a negative effect on the derivative estimation. Moreover, in some cases, the calculation of the coefficients of the polynomial tends to get numerically unstable, which is causing an unstable derivative estimation, too. The values in Table 4.2 show that the lowest error for the torsion estimation can be reached with Cesàro's mean. The mean error of the torsion estimation is more than three times smaller than for the partial sum and a magnitude smaller than applying the binomial filter. The curvature error is typically lower than the torsion error, because the curvature calculation needs not more than the 2nd derivatives, whereas the torsion calculation needs the 3rd one. Therefore, it is appropriate only to consider the torsion error to evaluate different low pass filters and estimation methods. The binomial filter has the lowest curvature estimation error and interestingly, the mean error of the torsion estimation with reflection is smaller than the error with restriction on f * (t). The latter 1.17 0.34 4.33 can be explained with the omission of the boundary pixels. However, the binomial filter will not be considered in the following sections anymore, because the high errors of torsion estimation. The higher error of the partial sum is mainly caused by the high standard deviation, whereas the mean torsion value is much closer to the true torsion value compared to Cesàro's mean. If another helix with the parameters r = 32 and c = 1 is used, the overestimation of Cesàro's mean results in a higher mean error δ8 = 6.21 than applying the partial sum δ4 = 1.62. Furthermore, the error of the binomial filter increases to δ8 = 6.71.

The Influence of Changes of Local Curvature and Torsion along the Curve

Contrary to the parametric function of helix f 1 , 3D point data extracted from real scenes (for instance 3D cell paths in open porous media) are characterized by changing curvature and torsion along the underlying curves. The Fourier approximation simultaneously smooths the curves and allows for the computation of the derivatives required for curvature and torsion estimation. It is clear that a constant smoothing parameter m will have a different impact in terms of smoothing and, consequently, in terms of accuracy according to the local curvature and torsion values. On the one hand, smoothing along the curve can reduce discretization noise and, on the other hand, this smoothing can also destroy the information on local curvature and torsion. The aim of this section is to study the latter effect on the estimation accuracy based on the cyclic function f 2 (t) given in Fig. Notice, that since the local torsion of f 2 (t) is π periodic, it is sufficient to plot the simulation results for 0 ≤ t ≤ π. The mean errors δm (see Eq. 4.16) of torsion estimates from the Fourier approximation are 6.39 • 10 -3 for the use of the partial sum with m = 4 7.25 • 10 -3 for the use of Cesàro's mean with m = 4 12.9 • 10 -3 for the use of the partial sum with m = 8 4.09 • 10 -3 for the use of Cesàro's mean with m = 8

These errors must be normalized by the maximal true torsion value 37.5 • 10 -3 known for f 2 (t). The smallest error for the partial sum is about 17 % with m = 4 and 11 % for the Cesàro mean with m = 8. The low error of Cesàro's mean is mainly caused by a small variance of torsion estimates τ 8 , which is not shown in Fig. 4.5. Further results are presented in Table 4.3.

In the most parameter configurations shown in Table 4.3, the partial sum has a lower error than Cesàro's mean. One explanation is, that Cesàro's mean is more sensitive to the smoothing parameter m than the partial sum. However, if its smoothing parameter is well adapted Cesàro's mean lead to the best results for the function f 2 (t). Clearly, for the partial sum, as well as for Cesàro's mean, the mean error decreases with an increasing scaling factor r or, equivalently, with an increasing sampling rate. On the other hand, increasing c results in lower errors for the partial sum, but causes higher errors for the Cesàro mean for r > 4. The last result is due to Cesàro's averaging which leads to larger estimation bias in cases of functions τ (t) with stronger slopes.

Fourier Analysis of the Discretization Noise

Now the effect of sampling a space curve on lattice Z 3 is considered in more detail. Let F (t) be a random function obtained by an IUR rotation and a random translation of a parametric function f (t) on the unit cell [0, 1) 3 of Z 3 , where the rotation and the translation are independent random variables. Furthermore, let ⌊x⌋ denote the floor function of a vector x = (x 1 , x 2 , x 3 ) ∈ R 3 defined as ⌊x⌋ = (⌊x 1 ⌋, ⌊x 2 ⌋, ⌊x 3 ⌋) with ⌊x k ⌋ = max{n ∈ Z : n ≤ k} for k = 1, 2, 3. Then the mean error ǫ(t) of the sampling of F (t) on Z 3 can be defined as

ǫ(t) = F (t) -⌊F (t)⌋, t ∈ [0, 2π]. (4.19)
One could expect that the spectrum of random function ǫ(t) only contains high frequencies, but this is not the case. In fact, the spectrum of ǫ(t) also contains low frequencies with considerable amplitudes. This is illustrated in Fig. 4.6. The reason for that is that the sampling of a realization of F (t) on Z 3 depends on the direction of pieces of this realization. For example, when the curve piece is parallel to the coordinate axes, it is sampled differently than if the same curve piece is oblique to the coordinate axes.

The random function ǫ(t) has a broad spectrum and the shape of this spectrum depends on F (t) itself. These facts have consequences for the possibility of separating noise from outer Jordan discretization of the random function F (t). Even if the realizations of F (t) are smooth enough and, thus, its spectrum would consist of only low frequencies, the noise induced by the outer Jordan discretization cannot be perfectly separated from the low frequencies belonging to the realizations by any kind of low pass filtering. Using the partial sum of the Fourier approximation, the higher frequencies which significantly affect the results are truncated such that the accuracy of the torsion estimation is preserved. The binomial filter does not completely suppress the higher frequencies which yield to a significantly worse estimation. Cesàro's mean should be the best approach because it attenuates the noise part corresponding to the low frequencies. However, for a small value of m the attenuation of low frequencies is too strong, such that this low pass filter yields in wrong derivations causing overestimated torsion, in particular for the helix as noted in Section 4.2.1. A larger smoothing parameter m can compensate this effect but introduces also a higher variance of the estimates. One alternative could be to design a low pass filter adapted to the discretization noise. The accuracy of Cesàro's mean depends strongly on the smoothing parameter and there are functions where a systematically overestimating cannot be avoided. Therefore, the partial sum will be used in all following sections for the Fourier approximation.

This section has presented the study required to both assess the effects of the curve discretization on the torsion estimation and to understand the effect of different boundary handling methods. The results were obtained for the Fourier approximation. A similar investigation has to be done for the method based on discretization of the differentialgeometric formulas. The investigation is presented in the next section.

Evaluation of the Discretization of the Differential-Geometric Formulas

In the following, the impact on the torsion accuracy of the discretization of the differential-geometric formulas given by Eqs. (3.14) and (3.15) is investigated regarding the sequences of particle positions. A discretization using the second discretization scheme with the step ∆t > 0 induces systematic errors in estimating τ , and one probably would expect that the systematic error increases with increasing ∆t. Furthermore, it must be taken into consideration that in the presented application the points on f (t) are particle positions estimated from 3D image sequences. These positions cannot be determined exactly, i.e. there are random measurement errors (displacements) of the true particle positions which can considerably increase the estimation errors of κ and τ . Clearly, random displacements lead to statistical errors, but the consequences on the torsion estimation accuracy of discretization effects (expressed in terms of ∆t) and of the errors in determining the particle positions (i.e. the displacement variation) are intertwined with one another in a complex manner.

In order to study the accuracy of torsion estimation, it should be assumed that τ is not estimated from the exact positions f (k∆t), but rather from the blurred positions x k = f (k∆t) + ε k , where ε k are independent and identically distributed displacements (i.e. localization errors). The independence of these displacements is justified by the fact that the distance of subsequent particle positions in two 3D images is much larger than the pixel size. Now, given that the components of the ε k are Gaussian distributed with expectation 0 and variance σ 2 , the error of an estimate τi of the true torsion τ depends on the true curvature κ, the true torsion τ , the step ∆t and the variance σ 2 of the displacements. The accuracy of estimating the local torsion τ is mainly discussed in terms of the relative error

δ = 1 nτ 2 n i=1 (τ i -τ ) 2 , (4.20)
where n is the number of torsion estimates at a given point on the curve. Notice that for curves with constant curvature and torsion, the relative error δ is independent of t and, thus, it is sufficient to estimate the torsion at one curve point n times. If κ and τ depend on t, then δ also depends on t, and n is the number of (noisy) functions to be simulated. Furthermore, since particles with speed close to the critical one (determined by the scanning rate) appear more or less blurred, the accuracy of the measured particle positions usually depends on the speed. However, the particles considered in this thesis are much slower, and thus, one can assume a constant σ.

Clearly, σ decreases with increasing lateral resolution, but this relationship is more complex and described in the second discretization scheme in Section 3.1.

First, the helix f 1 is considered. The local curvature and torsion of the helix are independent of t, which makes f 1 (t) interesting for studying errors of estimation depending on the true local curvature and torsion. complexity of the relationships between relative estimation error δ 1 of f 1 (t), the true curvature κ and torsion τ , the step ∆ς and the standard deviation σ. Notice, that the relative error δ 1 for -τ is the same as for τ , i.e. the diagrams shown in Fig. 4.8 can be reflected at the κ-axis. The admissible region in the κτ -space, i.e. the region for which the relative error is less than 20 %, can become very small. It decreases with decreasing ∆ς and variance σ 2 . Of course, these effects are less distinctive in the curvature estimation than in torsion estimation. A large ∆ς can lead to a high error of the τi , but in cases of sufficiently smooth curves (of small curvature and torsion), the choice of a moderate ∆ς can considerably reduce the error induced by the random displacements ε k , but in applications considered in this thesis, ∆ς cannot be chosen independently of the experimental framework. It is determined by the frame rate of the camera and the speed of the particles to be observed.

By contrast to the helix, the local curvature and torsion of other parametric functions depend on t. This variability is useful to study the more realistic situation of estimating the torsion at a curve point from neighboring points with different torsion. Indeed, the change of torsion along the curve can have a significant impact on the estimation accuracy. An example of a curve of the parametric function f 2 (t) is chosen. Thereby, the scaling factor r > 0 and the constant c ∈ R can be varied in such a way that curve pieces of f 2 (t) fit a wide range of possible particle paths. The curve length is given by ς

(t) = r 2 E(2t, -c 2 )
, where E is the (incomplete) Legendre elliptic integral of the 2nd kind given by Eq. (4.9).

For a careful investigation of the estimation errors of τ at t, it is necessary to fix the step ∆ς and, consequently, to switch to the length parameterized version of f 2 (t). This can be done by computing ∆t numerically from ∆ς based on the inverse t(ς) of the function ς(t). For this purpose the approximation

ς(t) ≈ rE(-c 2 ) 2t π + E π 4 , -c 2 2E(-c 2 ) - 1 4 sin 4t (4.21)
is considered which can simply be inverted by Newton's method, where E(x) = E π 2 , x is the complete elliptic integral of the 2nd kind defined in Eq. (4.10). The right-hand side of (4.21) forms a partial sum of the Fourier series of ς(t) completed by a linear term. This approximation is sufficiently accurate for ς(t) as 0 ≤ c ≤ 1. Fig. 4.9b shows that the estimation errors depend on t.

For a more detailed investigation of estimation errors of the local torsion, it should be referred to [An et al., 2011], where further classes of parametric functions are considered. In principle, however, it might be reasonable to restrict those investigations on the helix, since every sufficiently smooth function f (t) can locally be fitted with a helix f 1 (t), where the parameters r and c of the helix are chosen in such a way that its curvature and torsion are the same as that of f (t). As a consequence, the relative error δ of torsion estimation of f (t) should be close to that of the locally fitted helix. To make this clear, the projection of the parametric function f 2 (t) in the κτ -diagram of helices is considered, i.e. the curve κ(t), τ (t) , see Fig. 4.9a for an example. The relative error δ 2 of τ for f 2 (t) can roughly be estimated from the κτ -diagram by following the projection of f 2 (t). In Fig. 4.9b the relative error δ 1 along the projection of f 2 (t) is compared with δ 2 directly obtained from f 2 (t) by Eq. (4.20) using the estimator τ given by Eq. (3.15). It is noticeable that in this example both error estimates are close to each other. The difference is a consequence of the change of the true torsion along the curve of f 2 (t).

In summary, from the data presented in Fig. 4.8 it is possible to estimate the error of curvature and torsion estimation by the use of the discretization of the differential-geometric formulas nearly independently of the investigated function. In contrary to the Fourier approximation, there is no need for a boundary treatment. It is appropriate to exclude the boundary torsion values because they have no influence on the further estimation.

A Comparison Study

After investigating the curvature and torsion estimation methods regarding their boundary handling, sensitivity towards discretization noise and expectable errors, this section now focuses on the comparison of the performances of the proposed methods with those of the reference spline approximation. For this purpose, it is necessary to carefully select the smoothing parameter for the proposed methods, as well as for the spline method. This choice was already done in Section 3.5 for the discretization of the differential-geometric formulas through the adjustment of the consistency and the step width of the gradient filter.

On the Choice of Smoothing Parameters

An important aspect of this work is the choice of the smoothing parameters m and p of the Fourier approximation and the spline approximation, respectively. For their choice, no objective criteria exist and, therefore, they must be intuitively set. Nevertheless, there are significant differences in the use of the smoothing parameters. In 3D images, high frequencies, usually associated with the noise spectrum, can be separated from low frequencies of the spectra of sufficiently smooth particle paths with the careful choice of the Fourier smoothing parameter. In other words, the choice of an optimal m would be possible, if the path and noise have separate spectra.

As illustrated in Fig. 4.10a, for the discretized curve of the parametric function f 2 , the error of estimating the torsion is widely independent of m even if parameters r and c are varying in a wide range.

Parameter p of the spline approximation cannot afford the same as m. For all p > 0, high frequencies induced by discretization effects cannot be completely suppressed by a smoothing with splines, and the retained amplitudes of high frequencies have a great influence on the torsion estimation accuracy. Clearly, with decreasing p the amplitudes of high frequencies caused by the discretization noise decrease too but, simultaneously, the amplitudes of lower frequencies are reduced. As a consequence, the optimal p which minimizes the estimation error highly depends on the curve shape. This is visible in Fig. 4.10b showing that the range of p values minimizing the errors for the settings referenced by C1 and C2 is quite different from the p values leading to an accurate torsion estimation in parameter configuration C3. On the contrary, in Fig. 4.10a it can be noticed that all m values ranging approximately in [4, 9] lead to small errors for all three parameter settings C1, C2 and C3.

The Influence of Sampling

In this section, an instance of helix f 1 (t) is used to simulate sequences of particle positions classically occurring and extracted from sequences of 3D images. Such a discrete path is characterized by a small number of 20 3D points, a large sampling step width ∆ς ≥ 1 and a small displacement error ε < 1 • 10 -3 . The parameters of helix f 1 (t) are adjusted according to typical torsion and curvature values extracted (measured) from sequences of 3D images : τ = 0.1 and κ = 0.11. In the first simulation, the step width variance and displacement error ε are set to zero, which correspond to a particle motion with constant speed. For such curves, natural cubic splines can be used to pass through all discrete points, such that the adjustment of the parameters of smoothing splines can be avoided. In the second simulation, the step width is varied with a standard deviation of 0.5 and a displacement error of ε = 0.0005 is introduced. The results of both simulations are shown in 4.11a and 4.11b, respectively.

It has to be noted that boundary effects arise if a cyclic curve is approximated by a non-periodic cubic spline, where the second derivatives of the spline at t = 0 and t = t n-1 are set to zero (natural condition). These constraints can have a considerable impact on torsion estimation as shown in Figure 4.11: Torsion estimation of a helix with the discretization of the differential-geometric formulas (DD) and the spline approximation (SA). The real torsion values are depicted by (T). The helix is defined by its curvature κ = 0.11 and torsion τ = 0.1. In (a) the constant position distance of ∆ς = 8 is without displacement error and variation. In (b), the position distance is ∆ς = 1 with a standard deviation of 0.5. The displacement error is ε = 0.0005.

A similar effect can be observed for the Fourier approximation for which the small point number and the curve discontinuity result in an inaccurate approximation. The graph of the Fourier approximation is not shown in Fig. 4.11a because its error is multiple times larger than the error of the two other estimation methods. As illustrated by Fig. 4.11b, for all methods, the torsion estimation gets much worse with a varying step width. The Fourier approximation is not suited for varying step width and therefore, it is again not considered in Fig. 4.11b. It is very difficult to give an overview about the complex interdependencies between step width, step width variation, displacement error and curve shape represented by true curvature and torsion values. Therefore, only one curve example is shown in Fig. 4.11b.

A large smoothing value was used in Fig. 4.11b for the spline approximation (p = 1 • 10 -4 ) to attenuate the errors due to the varying step width. However, the discretization of the differential-geometric formulas offers the best results for both shown cases, without the need of adapting a smoothing parameter.

Sequences of particle positions are not useful to compare the discretization of the differential-geometric formulas and spline approximation with the Fourier approximation. For a comparison of all three estimation methods, a sampling scheme based on the outer Jordan discretization has to be used. Clearly, the discretization of the differential-geomet-ric formulas is motivated by the estimation of curvature and torsion from sequences of 3D images. However, single pixels have to be selected from the outer Jordan discretization to allow for a comparison with the Fourier approximation. As explained in Section 3.2.2, it is important to reduce the noise by regularization. Therefore, some pixels carrying the relevant information on the original curve shape (so-called informative points) are selected from the outer Jordan discretization. The selection rule is as follows: Initially, starting point x 1 and end point x n are connected by the line segment

[x 1 , x n ] = {x = qx 1 + (1 -q)x n , 0 ≥ q ≥ 1}. (4.22)
Then, one searches for a skeleton point x k for k = 2, . . . , n -1 with the largest distance between this segment and

x k d = dist([x 1 , x n ], x k ) ≥ dist([x 1 , x n ], x i ), ∀i = 2, ..., n -1.
The skeleton is divided into two skeletons from x 1 to x k and from x k to x n , respectively. Then, this subdivision is applied to both skeletons, etc. The recursive algorithm stops when the distance d is less than or equal to the pixel size. The splitting points are the searched points x k with concentrated curvature [START_REF] Usamentiaga | Fast and robust laser stripe extraction for 3D reconstruction in industrial environments[END_REF]. The torsion is only defined if κ(t) > 0. Only these points need to be considered for curvature and torsion estimation by the discretization of the differential-geometric formulas.

In the following, parametric functions f 1 and f 2 , both with the parameter r = 32 and c = 1, are considered. Helix f 1 has a mean torsion of τ = 15.62 for 2048 instances of the helix with about 830 points each. The spline approximation systematically underestimates the true torsion. This suggests that the smoothing is too strong, but lower smoothing results in a higher standard deviation and in an increasing mean error. For instance, a smoothing parameter of p = 0.2 • 10 -3 leads to a mean torsion value of 13.9 • 10 -3 instead of 11.52 • 10 -3 . However, the standard derivation of the mean torsion increases from approximately 0.02 to over 3. The discretization of the differentialgeometric formulas has the largest standard deviation and therefore, the largest error, but the mean torsion value of the estimation fits with the real torsion value.

An example of torsion estimation for the function f 2 is shown in Fig. 4.12. The spline smoothing parameter of p = 4 In Fig. 4.12a the confidence band of the discretization of the differential-geometric formulas gets slightly smaller than that of the spline approximation (see 4.12b) and that with nearly similar mean values. However, the Fourier approximation has the lowest mean error with by far the smallest confidence band. These results can also be verified in Fig. 4.10, where the error of the Fourier approximation with smoothing parameter values between 4 and 9 are far below the smallest errors of the corresponding spline approximation.

Conclusion

In this chapter, the two torsion estimating methods presented in Chapter 3 were investigated regarding their accuracy and boundary handling. They were compared to the commonly used spline approximation regarding the outer Jordan discretization and sequences of particle positions. The results have shown that the Fourier approximation is best suited for connected particle paths obtained from the outer Jordan discretization. Notice, that the discrete Fourier transform used in the Fourier approximation needs constant distances and for this reason, the Fourier approximation is not suitable for the estimation of the local curvature and torsion from sequences of non-equidistant particle positions.

Clearly, the partial sum is a typical low-pass filter applied in the frequency domain. Thus, one could expect that the use of any other low pass filter, such as binomial filters of higher order or Cesàro's mean, would lead to similar errors of torsion estimation, when strong smoothing with large filter masks can be performed effectively via frequency domain, only. However, the results mainly depend on the boundary handling and ability to cut off high frequencies. The binomial filter cannot completely suppress high frequencies and, therefore, it causes the highest mean estimation error. On the one hand, Cesàro's mean has a clear advantage when the parametric function is reflected on t = 0 but, on the other hand, the highest accuracies over a wide range of test functions is offered by the partial sum with the restriction on f * 1 (t). Moreover, for the most functions investigated with the restriction on f * 1 (t) the smoothing parameter can be fixed to m = 4 which leads to far better estimation accuracies than that of the discretization of the differential-geometric formulas and the spline approximation. In specific cases, Cesàro's mean with a well-chosen smoothing parameter can still increase this accuracy. However, in other cases in particular for the helix, Cesàro's mean systematically overestimates the curvature and torsion. A more theoretical advantage of Cesàro's mean is, that in addition to Fejér's theorem, the curvature and torsion estimation is multigrid convergent.

Other smoothing techniques like the use of smoothing splines do not completely remove the high frequencies induced by sampling the curve on a lattice. This observation is of importance since even small amplitudes of high frequencies, which are not removed by smoothing, can have a considerable impact on the third derivative used in the torsion estimation.

Generally, the estimation of the torsion of particle paths from sequences of particle positions has been proven to be a difficult problem in image analysis. It needs a careful choice of estimation methods depending on the specific kind of data sampling, as well as a forecast of the estimation errors to be expected. The core of the estimation method well suited for this kind of data is a new discretization scheme of the differential geometric formulas. This discretization can be directly induced by observing moving particles by time-resolved ➭CT, where the distances between particle positions (i.e. the step width) in two subsequent 3D images are large compared to the accuracy of measuring these positions which mainly depends on lateral resolution. Furthermore, the values of curvature and torsion have themselves a considerable impact on the estimation errors. Simulation studies are helpful to get an overview of the very complex interdependencies between these different factors influencing the estimation accuracy. For example, one can see from the diagrams in Figure 4.8 that in cases of low curvature and torsion a moderately increased step width ∆ς can compensate less measurement accuracy. Furthermore, no dedicated boundary treatment is needed for this estimation method. In opposite, the spline approximation wrongly estimates the real value because the boundary condition of natural splines affects the estimation accuracy of the whole curve. Moreover, the smoothing parameter of the spline approximation must be carefully adapted to the curve shape, whereas the used discretization of the differential-geometric formulas needs no smoothing parameter.

One aspect which was not considered so far is the algorithm speed or complexity of the estimation methods. In particular, for the 3D image processing, it is very important to reduce the amount of data to increase the processing speed. This reduction is mainly done during the particle path extraction and therefore, the time for torsion estimation makes only a small contribution to the overall processing time. Nevertheless, the computation time of the discretization of the differential-geometric formulas increases linearly with the number of points. The cubic spline approximation has linearly increasing time, too but as pointed out by [START_REF] Toraichi | Computational complexity of spline interpolation[END_REF] it needs 9 processing steps per sampling point and is, therefore, 9 times slower than the discretization of the differential-geometric formulas. The Fourier approximation has the highest complexity with O(n log n), where n is the number of sampling points.

In the next chapter, the two evaluated and well-investigated estimation methods are applied to real 3D data, where the true torsion is, of course, unknown. The information of the torsion distribution is used to evaluate the materials regarding their suitability for the chromatographic process.

Chapter 5

Experimental Results

In this chapter, two different kinds of open porous media are evaluated regarding their suitability for chromatography. Therefore, two experiments are designed to acquire the necessary 3D images and 3D image sequences. With this data, the well-established characteristics like porosity, specific surface area etc. as well as the curvature and torsion of the particle path can be computed and they are used to compare the different porous media. The content of this chapter is partly presented in [START_REF] Blankenburg | Torsion estimation of particle paths through porous media observed by in-situ timeresolved microtomography[END_REF].

Experimental Setup

In this thesis, two kinds of discretization schemes are used. Both schemes have their origins in different acquisition setups, which are described in the following sections. In the static experiment, 3D images are acquired without a suspension and, therefore, without any moving particle. The results are used to obtain fast particle paths regarding the outer Jordan discretization scheme. The setup for dynamic experiments is more complex and involves a suspension flow during the acquisition of 3D image sequences. In particular, for the second experiments, special measures in the CT setup are necessary to detect the moving particles in the suspension without disturbing the suspension flow. The resulting sequences of 3D images are used regarding the second discretization scheme.

Setup of Static Experiments

The collimated beam ➭CT setup of the ESRF beamline ID19 was used to acquire images of open porous media. Fig. 5.1 shows the experimental setup for the static measurement. The open porous media samples are woven and non-woven nylon fiber fleeces with a fiber radius of 20 ➭m, as well as polyurethane foams with a PPI value (see Chapter 2.2.2) of about 60 inch -1 (2.36 mm -1 ). The samples are irradiated by X-rays with a fixed energy of 19 keV, shown as a green beam in Fig. energy is λ = hc 0 /E = 65 pm. The Fresnel diffraction intensity pattern is projected onto a phosphate scintillator, which converts the X-ray into visible radiation. This light is collected, transmitted and focused on a 7 megapixel (MP) "Pco Edge Cameralink monochrome CCD" camera with 32 bit image depth. The field of view is 1.67 × 1.67 × 1.4 mm 3 giving an effective pixel size of 1.1 ➭m. Notice that the pixels in the acquired 3D images are cubic, so it is common only to give the length of one side as effective pixel size.

In further experiments, two cameras of the Pco company with 4 MP and 1.4 ➭m and 3.5 ➭m effective pixel size were used. The distance z between sample and scintillator is 30 mm. Consequently, the Fresnel number is d 2 /(λz) = 2.6 • 10 6 ≫ 1, where d = 1.1 ➭m•2048 is the effective detector size. This geometry is appropriate for the acquisition of phase contrast images. Table 5.1 gives an overview on the acquired 3D images, which are used in the following sections. A complete list of all 3D images can be found in appendix A.

Setup of Dynamic Experiments

A representation of the experimental setup for the dynamic experiments is given in Fig. 5.2. It consists of a porous media filling a cylindrical cartridge (column) which is connected to a syringe pump (Pharmacia Peristalic P-1). During in-situ experiments, a particle suspension was pumped into the column. The flow rate of the suspension was about 0.15 mL s -1 and the cross-section area of the specimen was 78.5 mm 2 , which yields to a rough estimate of the speed of the particles moving through the pore space. The mean vertical particle speed was about 1.9 mm s -1 . The suspension consists of spherical silica gel particles (SG20, Merck KGaA) with a diameter of 20 ➭m suspended in an ethanol solution, where the silicate particle serve as substitutes of living mammalian cells. A Heidolph Unimax 1010DT shaker was used in order to keep the silica particles in solution with a solid to liquid ratio of about 5 %. The displacement rate of the particles is relatively high, so an observation of fast particles in the pore space is not possible. Table 5.2 shows a summary of the acquired 3D image sequences. These image sequences were used for determining the sequences of particle positions. Notice that besides the described setup, further setups were investigated and the results obtained with them will be briefly outlined. For all experiments with suspension flow, it should be taken into consideration that the high energy of the radiation can lead to gas bubbles in the suspension, in particular, if a water-based solution like phosphatebuffered saline (PBS) is used. These gas bubbles considerably disturb the suspension flow. Moreover, in tests with living cells in a PBS solution, the cells were invisible in the X-ray images because the diffraction indexes between PBS and cells are too similar. Even adding dimethyl sulfoxide (CH 3 ) 2 SO or potassium ferricyanide K 3 [Fe(CN) 6 ] into beads with a diameter of 10 ➭m instead of cells do also not lead to better results in terms of contrast. For these reasons spherical silica gel particles were chosen, both for their cell-like shape and for their contrast with the ethanol solution.

Moreover, it is quite difficult to control and maintain the predefined flow conditions during image acquisition while ensuring high contrast between cells and solution. Therefore, samples of fiber fleeces were prepared before image acquisition. They were flown through by a suspension in controlled conditions and after a specific time a fixative was subsequently admixed until the cells were fixed and shape preserved. Afterwards, the whole sample was dehumidified to guaranty a high contrast between fixed cells and surrounding air. Unfortunately, the cells either did not attach themselves to the fibers or they were not affected by the fixative. Consequently, these cells were destroyed during dehydration or the adhesion between cell and fiber was nullified during fixation or dehydration. However, in the prepared samples only big clusters of fixative including multiple fibers are observable. It was not determinable if the clusters contain any cells. Nevertheless, all these challenges have been solved with the use of silica gel, an ethanol solution and a carefully designed ➭CT setup, which is presented in the next section.

Time-resolved Microcomputed Tomography

Time-resolved ➭CT (see Chapter 1.5) involving liquids is generally a challenge due to the high X-ray dose delivered to the sample. The protocol used in the frame of this work combines a well-balanced choice between data acquisition rate, X-ray photon energy used and spatial resolution reached. In order to reduce the dose to the sample while maintaining a homogeneous wavefront suited for the X-ray inline phase contrast [START_REF] Cloetens | Phase objects in synchrotron radiation hard X-ray imaging[END_REF]], beamline's u13 single-harmonic undulator (peak photon energy at 26.3 keV) was used at minimal gap (11.5 mm) and combined with a 1 mm-thick diamond and 0.7 mm-thick aluminum filter to suppress the softer parts of the spectrum. The absence of further optical elements besides the beamline exit window (made of 0.5 mm-thick beryllium) ensures high imaging sensitivity. In order to increase the photon flux density at the position of the detector (approx. 150 m downstream of the source) six so-called compound-refractive lenses made of beryllium were used to collimate the beam onto the field of view of the detector. The latter consisted of an indirect system combining a 25 ➭m-thick LuAG:Ce (Ce-doped Lu 3 Al 5 O 12 , Crytur, Czech Republic) single-crystal scintillator whose luminescence image is projected via a 10×/0.3 NA objective onto the chip of a CCD camera of type pco.dimax (PCO AG, Germany) having 2 016 × 2 016 pixels, 11 ➭m pixel size, 50 % peak quantum efficiency at 500 nm, 36 GB on-board memory for fast intermediate storage, and 1 279 full images per second (fps) maximum frame rate [START_REF] Douissard | A versatile indirect detector design for hard X-ray microimaging[END_REF]. The effective pixel size of the system is 1.1 ➭m, which yields to a field of view of 2.2×2.2×1.65 mm 3 . The configuration allows one to acquire 1 000 projection angles in 1 s with a contrast sufficient for tomographic reconstruction. Hence, 500 projection angles covering an 180 degree rotation were used to perform ➭CT with a temporal sampling of 0.5 s. A direct recording of the images on the hard disc is not possible at the used image acquisition rate. Therefore, the camera's on-board memory was utilized, leading to acceptable limits of tomographic scans which can be acquired in a row (500 scans/row). Tomographic scans were acquired every 20 s in order to cover a representative time interval appropriate to the dynamic changes of the cell flow. The contrast in the images was increased by using single-distance phase-retrieval with an adapted version of Paganin's approach (assumed ratio of 300 for the real to the imaginary part of the refractive index, 100 mm sample-detector distance, i.e. the Fresnel number is d 2 /λz = 1.2 • 10 6 ≫ 1, where λ is the wavelength, z is the propagation distance and d denotes the effective detector size) [START_REF] Weitkamp | ANKAphase: software for single-distance phase retrieval from inline X-ray phasecontrast radiographs[END_REF], Mirone et al., 2014].

To adjust the photon flux density of the time-resolved ➭CT setup one have to consider the attenuation coefficients of silica gel and ethanol at the used photon energy of 26.3 keV which are ➭ Sil (26.3 keV) = 4.00 cm -1 and ➭ Eth (26.3 keV) = 0.16 cm -1 [Sanchez del Río and Dejus, 2004], respectively. The corresponding difference in the coefficients makes it possible to derive the total number of 1 • 10 10 photons [Brooks andDi Chiro, 1976, Greaff and[START_REF] Greaff | [END_REF], which are required to resolve the two materials with the absorption contrast. Combining this value with the detector size corresponding to one tomographic slice with 2 016 pixels of size 1.1 ➭m, division of the photons over 500 projections, a detector efficiency given by the stopping power of the scintillator (25 %), a transmission dominated by ethanol, and a total acquisition time of 0.5 s per tomographic scan one ends up with a needed photon flux density of 10 12 photons/(s mm 2 ). This order of magnitude matches with the photon flux density of the above mentioned configuration of beamline ID19. Nevertheless, it is a number which only states a minimum requirement for density resolution but does not reflect the much higher demands on the signal-to-noise ratio in order to perform a segmentation required for subsequent image analysis methods. Due to the non-linear dependence of the absorption contrast to the number of photons required, a substantially higher contrast based on the attenuation behavior would already exceed the available photon flux at beamline ID19. Here, the above mentioned phase-contrast imaging approach in combination with phase-retrieval techniques allows one to boost the contrast by roughly an order of magnitude [START_REF] Sanchez | Three-Dimensional Synchrotron Virtual Paleohistology: A New Insight into the World of Fossil Bone Microstructures[END_REF].

About 70 3D image sequences containing each twenty 3D images with a size of 16 GB were acquired with the setup described in the previous paragraph. Theses sequences were used to investigate the motion of slow particles which can be tracked in successive 3D images. Together with the 50 3D images obtained from the static experiments, the partially open foams and the fiber fleeces can now be evaluated regarding their suitability for the chromatographic filtering.

Characterization of Partially Open Foams

In this section, partially open foams provided by Merck KGaA are investigated and characterized regarding their suitability for the chromatographic process. Therefore, the characteristics presented in Chapter 2 were estimated. The errors of the estimates are computed by varying the binarization threshold and the lateral resolution.

The estimated porosity of the foam is about 97 %. This value is widely independent of the lateral resolution. The mean volume V , the mean surface area S and the mean size b of the cells are V =0.142 mm 3 ± 0.091 mm 3 , S =1.702 mm 2 ± 0.813 mm 2 and b =700 ➭m ± 150 ➭m, respectively, where the mean size, is the average of the mean width b over all cells of the image.

The mean diameter of the spheres with equivalent volume is 600 ➭m ± 190 ➭m. The reciprocal value in inch -1 is called PPI value see Eq. (2.6), which yields PPI= 63 inch -1 (2.48 mm -1 ). With the alternative relationship N L = S V /4 one gets PPI= 62 inch -1 (2.44 mm -1 ) and using the mean cell volume and the surface area PPI= N L = S/4 V = 75 inch -1 (2.95 mm -1 ). These PPI values are within the tolerance of the manufacturer's instructions. This observation shows that the binarization and the image enhancement (with a 3 × 3 × 3 median filter) leads to appropriate 3D images for fast path extraction.

Estimates of the surface density strongly depend on the lateral resolution. These The mean size expressed in terms of the average of the mean width of the removed facets is 388 ➭m.

The histogram of the spherical contact distribution is shown in Fig. 5.3a. The peak of the distribution is approximately at the radius of 150 ➭m. This value can be interpreted as the mean size of the strut system including nodes. The histogram of the granulometry distribution is shown in Fig. 5.3b. The peak of the distribution is approximately at the radius of 300 ➭m, which corresponds to the mean diameter of the spheres with equivalent volume. This means that the pores are more or less spherical. From the percolation probability shown in Fig. 5.4, it can be seen that particles smaller than 1985 mm completely percolate the pore space. However, the partially open foam becomes impermeable for a particle size larger than 3.308 mm. These thresholds are widely independent of the percolation length.

The shape of cross-sections of the struts is sketched in Fig. 5.5, which shows that the cross-section is concave almost everywhere. The typical minimum principal curvature is two or three times larger than the PPI value. It ranges from -12.6 mm -1 to -7.4 mm -1 . The maximal principal curvature of the struts is close to zero. The mean principal curvature is -6.2 mm -1 ± 2 mm -1 . Close to the nodes of the foam, the principal curvatures are approximately the same, so a mean principal curvature of -6.4 mm -1 ± 0.8 mm -1 can be archived. The fast particle paths are obtained by an adapted version of the FMM, which is applied to a single 3D image as described in Section 3.2.1. The principal motion direction is chosen to be the z-direction. The extracted paths correspond to an outer Jordan discretization of a continuous curve. Therefore, the Fourier approximation with a partial sum as low pass filter and a smoothing parameter of m = 4 are selected to estimate the curvature and torsion of the particle paths.

Paths of Fast Particles in 3D Images

The local curvature and torsion distribution of the particle paths are shown as histograms in Fig. 5.7. The total number of values was about 12 500 for both the curvature and torsion estimates. Points with zero curvature were not considered in the torsion histogram. The mean curvature is usually greater than zero since the curvature is an unsigned characteristic. The mean curvature is about 1.2 mm -1 with a standard deviation of 0.7 mm -1 . The torsion distribution is approximately symmetric; the mean torsion is 0.148 • 10 -3 ➭m -1 and the standard deviation is 8.12 • 10 -3 ➭m -1 . Notice that the symmetry of the torsion distribution is a necessary condition for the distribution of macroscopically homogeneous structure to be invariant with respect to reflection ac- cording to plane z = 0. In other words, a flow along a principal direction parallel to the z-axis (from "top to bottom" or vice versa) is neither left nor right rotating.

A low torsion of the particle path results in a low microcentrifugal effect of the liquidsolid flow, which has an adverse effect on the chromatographic separation ability of the foam. On the other hand, the foam is featured with a concave surface and has, therefore, a large contact area for the convex cells. The influence of this surface on the paths of slow particles near the surface is investigated in the next section.

Path of Slow Particles

The positions of particles were iteratively tracked from one 3D image to the next in order to extract the particle paths of slow particles near the surface from a 3D image sequence. Fig. 5.8 gives an example of a sequence of particle positions. Only slow particles were observable due to the moderate temporal resolution of 0.5 s. The majority of the particles were fast and, consequently, they are invisible in these sequences of 3D images. For moving particles a maximum of twenty different positions per particle was detectable. However, for most particles this number is much smaller due to the fact that particles can partially move outside the field of view (boundary effects). Furthermore, some of the low energy particles deposited on the inner surface during the experiment or just fixed particles can start to move. , where the mean distance between the particle positions was 8.34 ➭m and a mean particle speed was 16.7 ➭m s -1 . The latter is less than 1 % of the maximum particle speed that can be reached far from the inner surface. As pointed out in Section 5.2, the partially open foam investigated in this study is non-oriented, i.e. neither left nor right rotating and, therefore, the torsion distribution is symmetric (in particular with zero mean torsion). As a consequence, it is reasonable to depict the histogram of the absolute values of torsion estimates. The torsion distribution represented in the histogram given in Fig. 5.9a has a mean value and an empirical standard deviation of 0.0022 ➭m -1 and 0.1792 ➭m -1 , respectively. An impression of the relative errors of the torsion estimates τi is given in Fig. 5.9b, where the points (κ i , τi ) are plotted in a κτdiagram introduced in Chapter 4.3 with appropriately chosen ∆ς and σ.

Until now, particles with a speed higher than about 50 ➭m s -1 are not observable with the time-resolved ➭CT technique. Therefore, fast particle paths were investigated in the previous section. It turns out that the torsion of the fast particles is half a magnitude smaller than that of slow particles. This means that the shape of the inner surface induces high path torsion. This shows that the microcentrifugation resulting from a high torsion has a considerable impact on the chromatographic efficiency of the porous medium. Indeed, the empirical mean curvature was 0.16 ➭m -1 for the slower particles (mean speed of 16.7 ➭m s -1 ) and 0.013 ➭m -1 for the faster ones (with a mean speed of at least 1.9 mm s -1 ). Hence, the magnitude F = mκ ḟ of the centripetal force F that makes a fast particle of mass m follow a curved path f is much higher than that for slower ones. Nonetheless, the moderate increase of the curvature (and torsion) for decreasing speed is of high importance for the microcentrifugation effect. It partially compensates the effect of slowing down the centripetal force such that significant microcentrifugation probably appears as long as the particle is not deposited on the inner surface.

Characterization of Fiber Fleeces

An example of the fiber fleece with fast particle paths is shown in Fig. 5.10. Fiber fleeces are porous media with a geometry differing completely from that of the partially open foams. However, most characteristics of these two kinds of porous media can be calculated in the same way. The porosity of the investigated fleece is about 59.7 % ± 0.2 % and the specific surface area is 84.3 • 10 -3 ➭m -1 ± 2 • 10 -3 ➭m -1 . The fibers themselves have a mean width of about 41 ➭m and their cores have a mean curvature of 0.05 ➭m -1 .

The first principal curvature of the filter surface is positive almost everywhere and the (Schladitz density) with the parameter β ≈ 2. Here ϑ and ϕ are azimuth and zenith, respectively, such as (1, ϑ, ϕ) is the random direction vector given in spherical polar coordinates. This type of a density function is used to characterize the fiber direction distribution of fiber fleeces produced by a melt-blown process. The parameter β was roughly estimated as described in [Ohser and Schladitz, 2009, p. 253ff]. More details about the direction distribution are shown in Fig. 5.11a. Here, a discrete version of the inverse cosine transform is used for estimating the histogram of the fiber direction distribution. The corresponding directions are visualized in Fig. 5.11b. Clearly, due to the production process, there are no fibers along the z direction. About 75 % of all fibers are lying in the xy-plane, which explains the invariance with respect to rotations around all straight lines parallel to the z-axis.

The histogram of the spherical contact distribution of the solid matter is shown in From the percolation probability shown in Fig. 5.13, it can be seen that the particles, which completely percolate the pore space of the fleece, are much smaller than in the partially open foam. For particles larger than 0.491 mm a strong dependency on the percolation length can be seen. Particles larger than 0.819 mm cannot percolate the fleece. However, the investigated particles have a size of 20 ➭m and can, therefore, completely percolate the pore space.

The curvature and torsion distribution of the particle paths through the pore space are shown in Fig. 5.14. The total number of curvature estimates was about 49 000. The torsion was estimated from the same path elements, but those paths elements with nonexistent torsion (zero curvature) have not been taken into consideration for estimating the histogram of torsion. The mean curvature is about 7.75 mm -1 ± 4.26 mm -1 , which is seven times higher than that of the partially open foam (1.2 mm -1 ). The mean torsion is -3.26 • 10 -3 ➭m -1 and its standard deviation is 59.18 • 10 -3 ➭m -1 .

Generally, the geometry of the random porous media influences the curvature and torsion distribution of particle paths through the pore space. Thus, it is not very surprising that the curvature and torsion distribution of the fiber fleece differ from that of ). First, it should be noticed that the variance of the curvature of the particle paths through the pore space of the fiber fleece is more than 35 times higher than that of the partially open foam. The variance of the torsion of particle paths through the pore space of the fleece (59.18 mm -1 ) is about 7 times higher than that of the partially open foam (8.12 mm -1 ). These results can be explained physically with the much smaller porosity and a much higher specific surface area of the fiber fleece compared to the partially open foam. On the other hand, one should take into consideration that, for a (theoretical) system of parallel fibers with directions parallel to the flow direction, the curvature and torsion would be zero everywhere even in cases of low porosity. Thus, the variance of curvature and torsion depends on the fiber direction distribution as well. Furthermore, the torsion distribution is not symmetric (i.e. not invariant with respect to a reflection at the origin). In Fig. 5.14b one can observe a significant shift of the fiber histogram to negative torsion values, while that of the partially open foam (Fig. 5.7) is nearly symmetric. In Fig. 5.10, the fluid flow through the pore space from top to bottom is left rotating. This means that the fiber fleece is not invariant with respect to reflection at the xy-plane. This effect is caused by the melt-blown production process of the fiber fleece, where the lattice of fiber spinnerets is rotated. Even if the fibers are strongly swirled during their down-laying, the fleece inherits the rotation direction of the fiber spinnerets' lattice. The consequences are twofold:

❼ The torsion distribution is a quantity which allows for the detection of asymmetries of macroscopically homogeneous structures. ❼ The torsion distribution of paths through the pore space of porous media can be influenced significantly by their production processes.

A high specific surface area and a high torsion of particle paths can significantly increase chromatographic cell separation. On the other hand, the convexity of inner surfaces can decrease chromatographic cell separation, since only surface elements with low principal curvature (seen from the solid matter) are preferred sites of cells.

To evaluate which porous media is better suited for the chromatographic process, the most relevant characteristics should be recalled: ❼ A higher surface area provides a higher number of potential particle sites to be deposited at the inner surface.

❼ Surface elements with a negative mean curvature are preferred sites of deposited particles. The density of the particles deposited on the inner surface increases with a decreasing mean curvature. The reason for that is that the electrostatic interaction and the short-range interaction between the spherical particle of a given radius r and a chemically activated surface element of principal curvature -1/r are very high. From the sequences of 3D images, it was observed that close to the nodes of the foam (where the two principal curvatures are negative) the density of deposited particles is much higher than on the struts (where only one principal curvature is negative and the other one is about zero).

❼ High local curvature and torsion of particle paths induce centripetal forces that push particles to the inner surface where they can be deposited.

The multiple interactions between the area of the inner surface, its principal curvatures and the local curvature and torsion of particle path are very complex. It is not yet known, which combination of these characteristics leads to high chromatographic effects. It is expected that a fiber fleece with a fiber cross-section as shown in Fig. 5.15c (but with the same mean width of the cross-section, fiber density and fiber direction distribution of the current fleece) could have a better filter effect than the partially open foam investigated in Section 5.2.1.

A Random System of Overlapping Fibers

The aim of this section is to show that the asymmetry of the torsion distribution of particle paths of the fiber fleece considered in Section 5.3 can be induced by the microstructure of the porous medium. For this reason, a simple geometric model is constructed, from which one could expect such an asymmetry. A random porous medium is considered, where the solid matter forms a macroscopically homogeneous system of overlapping fibers.

Let Φ = {C k } k∈N be a random system of helices with the random parametric function

g k (t) = f 1 (t -t k ) + x k , k = 1, 2, ...
with f 1 as given in Eq. (4.1) and parameter r > 0 and |c| > 0. The t k are random variables uniformly distributed over interval [0, 2π] and the x k are the points of a macroscopically homogeneous (stationary) Poisson point process with the constant parameter λ in the plane z = 0, where λ is the intensity, i.e. λ is the mean number of points per unit area.

The system Φ is macroscopically homogeneous and invariant with respect to rotations around all straight lines parallel to the z-axis, but it is not isotropic and also not invariant with respect to the reflection of the plane z = 0. In the following, macroscopically homogeneous but anisotropic random structures fulfilling these properties are referred to as "orientated random structures" (left-rotating). From Eq. (4.1) it follows immediately that the specific curve length L V of Φ is

L V = λ √ 1 + c 2 |c|
(mean total length of the curves per unit volume). Furthermore, for each t, the fiber directions

ḟk (t) = r    -sin(t -t 0 ) cos(t -t 0 ) c   
are uniformly distributed on a circle of the upper half sphere. More precisely, the altitude ϑ of ġk (t) is constant,

ϑ = arccos |c| √ 1 + c 2
while the longitude ϕ of ġk (t) is uniformly distributed on [0, 2π].

The dilation Ξ = Φ ⊕ B ̺ of Φ with a ball of radius ̺ > 0 forms a fiber of radius ̺ > 0. In order to exclude self-intersection of the parallel sets of the random helices it is assumed that ̺ is less than the curvature radius r/(1 + c 2 ) of the helices and the (deterministic) distance

d(t) = g k (t) -g k (0) = f 1 (t) -f 1 (0)
is at least 2̺ for all t > 2π. Thus, one gets

d(t) = r 2 -2 cos t + c 2 t 2
(5.1) and

ḋ(t) = r(c 2 t + sin t) √ 2 -2 cos t + c 2 t 2 .
This means that the function d(t) takes its minimum if

c 2 t + sin t = 0 or c 2 = sinc t, 2π < t ≤ 5 2 π.
(5.2) Now, inserting (5.2) into (5.1) yields the minimum distance

d(t) = r √ 2 -2 cos t + t sin t
at t = t 0 , where t 0 is the solution of Eq. (5.2) for a given c with 0 ≤ |c| < 2/(5π). For |c| ≥ 2/(5π) it is sufficient to suppose that ̺ is less than the curvature radius of the helices.

Finally, the volume fraction V V of the random fiber system Ξ is given by

V V = 1 -e -π̺ 2 λ √ 1+c 2 |c| .
Fig. 5.16 shows a realization of a random fiber system Ξ with λ = 2.5 • 10 -3 ➭m -2 , r = 192 ➭m, c = -0.4 and ̺ = 5 ➭m. The specific curve length of the core system Φ is L V = 7.25 • 10 -3 ➭m -2 and the volume fraction of the fiber system Ξ is 43.4 %. The measured specific surface density is 133.11 • 10 -3 ➭m -1 .

The fact that the random structure is left-rotating is inherited to the flow through the pore space. This can be shown by simulating fast particle paths using the FMM for flow from top to bottom. For example, the torsion distribution of the image shown in Fig. 5.16 is presented in Fig. 5.17, which demonstrates that in mean the particle paths are left-rotating, too. Figure 5.17: Histogram of the torsion distribution of fast particle paths through the pore space of the random fiber system shown in Fig. 5.16, where the particle paths are simulated by FMM.

Conclusion

In-situ time-resolved ➭CT has been proven as a promising technique for observing moving particles in a liquid-solid flow of a suspension through porous media. As far as known, the scan rate combined with high spatial resolution and sensitivity for time-resolved ➭CT reached at the beamline ID19 of the ESRF is currently one of the best. Using this technique it is possible, for example, to track the motion of particles with diameters larger than 20 ➭m and a speed less than about 50 ➭m s -1 . Future developments in this field will lead to much higher scan rates allowing also for tracking of faster particles. Clearly, simulation of liquid-solid flow is a well established supplementary method of investigating chromatographic filtering. Such flow models can, for instance, include the formulation of blood rheology depending on the local concentration of the cells, as well as its change by selective deposition [START_REF] Ju | A review of numerical methods for red blood cell flow simulation[END_REF], which can notably be of interest in the leukemia context. Nonetheless, experiments with real suspensions pumped through porous media and direct observation of particle motion and deposition will also be indispensable methods to get a deeper insight into filtration processes in the future.

Generally, the estimation of the torsion of particle paths from time-resolved ➭CT images has been proven to be a difficult problem of image analysis. It needs a careful design of the experimental setup, a choice of estimation methods depending on the specific kind of data sampling, as well as a forecast of the estimation errors to be expected.

The torsion distribution of particle paths is closely related to the orientation of the pore space, where "orientation" has the following meaning in the frame of this thesis: Assume that the distribution of the medium is invariant with respect to rotation around the z-axis. In this case, each section of the medium with a plane of the normal direction parallel to the z-axis forms random planar structure, which is macroscopically homogeneous, as well as isotropic. If the principal flow direction is parallel to the z-axis (flow "from top to bottom"), then orientation means that, on average, the flow through the medium is left or right rotating. Currently, the torsion distribution of particle paths through a macroscopically homogeneous porous medium is the only geometric quantity for detecting such an orientation.

The foam structure investigated in this chapter is completely isotropic. More precisely, the distribution of the foam is invariant with respect to spatial rotations and, thus, it is also free of orientation. However, orientation can be detected for the investigated fiber fleece produced by a melt-blown process with a rotating lattice of fiber spinnerets. This is an important aspect for the development of porous media applied in cell chromatography since one can considerably increase the microcentrifugation effect when inducing orientation of the pore space by varying the production process.

The centripetal acceleration caused by the specimen rotation during image acquisition can have a huge effect on curvature and torsion estimation. In the presented experimental setup for dynamic experiments, the specimen's radius and the angular speed were r = 4 mm and ω = 2π s -1 , respectively, which yields to a magnitude of the centripetal acceleration of at most 158 mm s -2 . On the other hand, the magnitude of the centripetal acceleration of a particle moving along the path f is equal to m f . For the slow particles, the mean of this magnitude was about 200 ➭m s -2 only. The impact on torsion estimation is not clear, but it can be argued as follows: In the case of a significant influence of sample rotation on the torsion of particle paths through the foam, the torsion distribution would be shifted (to left or right, depending on the orientation of the specimen rotation). The observed torsion distribution, however, is symmetric (and in particular, the estimated mean torsion is close to 0). Nevertheless, in future experiments with higher temporal resolution, the effect of the specimen rotation on the observations must be taken into consideration. Finally, it should be referred to dynamic experiments observed by laboratory time-resolved ➭CT, where the radiation source and the scanner are rotated around the specimen [Cnudde andBoone, 2013, Dierick et al., 2014], but until now, lateral and temporal resolution, as well as contrast reached by those setups, are much lower than in the experiments described in this chapter.

Publication relating to the work described in Chapter 5 [START_REF] Blankenburg | Torsion estimation of particle paths through porous media observed by in-situ timeresolved microtomography[END_REF]. Torsion estimation of particle paths through porous media observed by in-situ time-resolved microtomography. Journal of Microscopy, 266(2):141-152

Conclusion and Perspectives

The main aim of this thesis was to develop methods to characterize porous filters regarding their suitability for cell chromatography according to the geometry of the pore space. Therefore, experiments were designed at the European Synchrotron Radiation Facility in Grenoble. Particular attention was paid to the time-resolved micro-tomography. The proposed setup is probably one with the highest temporal resolution so far. An essential point in the experiments was the use of an alcohol suspension with silica particles. They offer an appropriate contrast with the used radiation of 26.3 keV and the application of inline phase contrast. The Paganin's phase retrieval method was applied to the 3D data reconstructed with the inverse Radon transform.

The curvature and torsion of particle paths through the media are two important characteristics of porous media regarding their suitability for the chromatographic process. The extraction of the fastest particle paths from a single 3D image can be performed by an adapted version of the geodesic dilation. In particular, the geodesic dilation was applied after the Euclidean distance transform of a 3D image. Therefore, the Euclidean distance of the inner surface is interpreted as a scalar speed function. It was shown that the geodesic dilation can be seen as a solution of the eikonal boundary value problem. Therefore, the Fast Marching Method can be used to increase the computation speed of the simulated particles paths through the medium.

The accurate estimation of local curvature and torsion of the fastest paths without prior knowledge about the curve shape is a challenging task. With the proposed Fourier approximation it was possible to decrease the torsion estimation error of at least 65 % compared to the commonly used spline approximation. Therefore, the discrete curve is decomposed into a sum, where the start and end points of the first summand are equal and the second summand represents a straight line. The Fourier approximation is restricted to the first summand, whereas the derivatives of the second summand can be analytically determined. The derivatives of the discrete curve are simply the sum of both summand derivatives. It was shown that this boundary treatment, together with the partial sum as low pass filter, produces the smallest mean estimation error. Depending on the curve shape and the selected smoothing parameter, the errors of this method are multiple times lower than that of the reference spline approximation. Moreover, the smoothing parameter m = 4 of the Fourier approximation can remain unchanged for both a wide range of lateral resolutions and curvatures and torsion values, whereas the smoothing parameter of the spline approximation must be adjusted to the resolution and curve shape.

While the Fourier approximation suits for particles paths extracted from single 3D images, it cannot lead to accurate results when points are non-uniformly sampled during a tracking of particle position in sequences of 3D images. Therefore, a second curvature and torsion estimation method based on the discretization of the differential-geometric formulas was developed. The discretization of the differential-geometric formulas can be directly applied to sequences of particle positions. In opposite to the spline and the Fourier approximation, the discretization of the differential-geometric formulas needs no particular boundary treatment so that the boundaries do not longer affect the estimated curvature and torsion values. The spline approximation and the discretization of the differential-geometric formulas led to similar mean errors. But the major difference between the two methods is that the filter mask remains unchanged for the discretization of the differential-geometric formulas, whereas the smoothing parameter of the spline approximation must be adapted to the curve shape and resolution of the function. It has been shown that the accuracy can be improved in particular for high step widths with the choice of a higher consistency. Unfortunately, this implies the use of 3 additional positions, which is inappropriate in the presented application. The curvature and torsion values of the curve have themselves a considerable impact on the estimation errors. Studies based on simulated data are helpful to get an overview of very complex interdependencies. For example, one can see from the diagrams in Fig. 4.8 that, in cases of low curvature and torsion, a moderately increased step width can compensate less measurement accuracy. It has been shown in Section 4.3 that these diagrams are valid for arbitrary functions as long as the function can be locally approximated with a helix.

The investigation of the fiber fleece revealed that the pore space of the fiber fleece is left orientated. "Orientation" means that, on average, the flow through the medium is left or right rotating. Currently, the torsion distribution of fast particle paths through a macroscopically homogeneous porous medium is the only geometric quantity for detecting such an orientation. Furthermore, compared to the partially open foam, the fiber system is characterized by a higher variation of the curvature and torsion values, as well as by a larger inner surface. Therefore, it can be expected that the investigated fiber fleeces are better suited for the chromatographic processes. A further improvement would be the use of a system of fibers with almost everywhere concave cross-sections.

Perspectives

Improvements of rapid prototyping via 3D printing open up a new range of possibilities to design chromatographic filters [START_REF] Fee | 3D printed porous media columns with fine control of column packing morphology[END_REF]. A first step to a filter with high filter efficiency could be the design of a Weaire-Phelan foam with a high surface area, high curvature and torsion values of simulated particle paths.

In this thesis, the simulation of the particle paths is based on geodesic dilation. However, despite the immense computational effort and some open problems implied by this technique, it would be very helpful to create a two-phase flow simulation with high temporal resolution through the acquired 3D images. The resulting discrete particle paths could provide a complete overview of the flow conditions. On the one hand, the proposed simulation of fast particle paths could be evaluated regarding accuracy and estimation error and, on the other hand, the influence of curvature and torsion of the particle paths on the filter efficiency could be quantified. Moreover, with such simulations it would be possible to determine the complex interaction of specific surface area, principle curvature of surface elements and the curvature and torsion of particle paths. However, such a simulation should be accompanied by an optimization of the algorithm regarding the computational time.

Clearly, discretizations of curves and, in particular, the outer Jordan discretization induce high frequencies in the spectrum of the discrete curves. However, as the investigations of the presented thesis show, discretization can also induce low frequencies of considerable amplitudes. So it could be expected that a low pass filter, which attenuates also the low frequencies lead to lower estimation errors. Such a low pass filter is realized by Cesàro's mean and indeed Cesàro's mean has a very tight confidence band. However, the drawback of Cesaro's mean low pass filter is the overestimated mean value, which can probably be corrected with the choice of more appropriate weighting factors.

In this thesis, the investigation of the accuracy of estimating curvature and torsion is mainly focused to helices. It was pointed out that these results can be extended to functions with changing curvature and torsion by approximating the functions locally with helices. However, this only applies to functions with small changes of curvature and torsion. The investigation of stronger changes, which was not a topic of the present thesis, should be a subject of future work.

Further Fields of Application

Besides the context of cell chromatography, where the curvature and torsion have a considerable impact on the filter efficiency, there are further fields of applications, where the accuracy of the presented curvature and torsion estimation methods can also have beneficial effects.

In materials sciences, it is well known for a long time that the torsion of dislocation lines in a crystal is closely related to the density of inflections that are optically discernible in transmission electron microscopy [Rhines, 1977]. This fundamental perception traces back to [Fullman, 1953]. Since this time there were published a continuously growing number of articles on the investigation of the impact of torsion of dislocations on macroscopic tensile strength [START_REF] Zhao | Asymptotic behaviors of the stress fields in the vicinity of dislocations and dislocation segments[END_REF] and, vice versa, of the torsion tensor of the crystal on dislocation motion. Furthermore, the curvature and torsion of paths through the pore space are of particular interest in nuclear research. As an example, on can consider the porous beryllium used in the helium cooled pebble bed (HCPB) blankets investigated in the framework of the European Fusion Technology Programme. The long-term tritium and helium accumulation in HCPB blankets is crucial for the reliable and safe operation of fusion reactors. It depends on the percolation, the curvature as well as on the torsion of the pore space [START_REF] Möslang | Gas bubble network formation in irradiated beryllium pebbles monitored by X-ray microtomography[END_REF].

As pointed out in [START_REF] Bae | Design of Battery Electrodes with Dual-Scale Porosity to Minimize Tortuosity and Maximize Performance[END_REF], low torsion can have a considerable impact on carrier mobility and, thus, on cell charging and discharging of batteries with porous electrodes. Finally, it should be remarked that mechanical properties like strength and stiffness of fiber reinforced materials depend on the torsion of the fibers' center lines.

Many applications in medicine are focused on the geometrical characterization of arteries. The relationship between the torsion of the internal carotid artery in human infants and the brain temperature is studied in [START_REF] Meng | Three-dimensional description and mathematical characterization of the parasellar internal carotid artery in human infants[END_REF]. A large absolute torsion correlates with the occurrence of vascular pathologies. Similar investigations of the torsion based on tomographic data are made for the coronary artery [START_REF] Pao | Bending and twisting of an in vivo coronary artery at a bifurcation[END_REF], Puentes et al., 1998, Medina et al., 2004, Zhenga and Qib, 2011, Strandmark et al., 2013] and for the superficial femoral arteries [START_REF] Wood | Curvature and tortuosity of the superficial femoral artery: a possible risk factor for peripheral arterial disease[END_REF].

The presented curvature and torsion estimation methods, can make a valuable contribution to all these different applications and in particular to the design of highly efficient chromatographic filters for leukemia treatment. A.5 Fibers of Filter Cake

In this section, the data obtained from a woven "filter cake" is summarized. The aim was to qualify if the filter characteristics are changing on different positions. Therefore, three pieces at three different positions were cut from the filter cake. The pieces were marked with a color code to identify the origin position in the "filter cake". Each piece was investigated on at least three different positions. Mots-clés: estimation de la torsion, chromatographie cellulaire, tomodensitométrie à résolution temporelle, géométrie différentielle, approximation de Fourier, simulation de trajectoire 3D.

Abstract

The extraction of cancerous cells from body fluids is an important procedure in clinical diagnostics and therapy. Notably, when the separation technique is based on cell chromatography, it is important to have precise knowledge about binding capacities of target cells in porous media. Therefore, experiments using time-resolved micro-computed tomography were designed and carried out at the European Synchrotron Radiation Facility. The curvature and torsion distributions of cell paths in a two-phase flow through a porous medium are valuable information to characterize the efficiency of chromatographic processes. However, the computation of torsion is very challenging, since it is based on higher order derivatives which are very sensitive towards discretization noise. In this thesis, two new curvature and torsion estimation methods of particle paths are presented. The first method is based on a Fourier approximation. Case studies showed a decrease of the torsion estimation error of at least 65 % compared to the commonly used spline approximation. Moreover, the smoothing parameter of the Fourier approximation can remain unchanged for both a wide range of lateral resolutions and curvatures and torsion values. Since this Fourier approximation approach cannot be applied at non-equidistant points, a second method based on the discretization of the differential-geometric formulas (DDGF) was developed. The spline approximation and the DDGF led to similar mean torsion errors. However, the filter mask remains unchanged for the DDGF, whereas the smoothing parameter of the spline approximation must be adapted to the curve shape and discretization.

Keywords: torsion estimation, cell chromatography, time-resolved micro-computed tomography, differential geometry, Fourier approximation, 3D cell path simulation.

Figure 3 :

 3 Figure 3: Dispositif pour détecter les trajectoires des particules. (1) Erlenmeyer avec la suspension, (2) Incubateur Heidolph Unimax 1010DT, (3) Pompe de type "Pharmacia Peristalic P-1" et (4) support de l'échantillon avec le milieu poreux.

Figure 4 :

 4 Figure 4: Exemple d'extraction des données 4D (séquences d'images 3D) des chemins des particules lentes. La particule suivie est en rose et est montrée à différents pas temporels.

  existent, òu c ℓ = (c ℓ,1 , c ℓ,2 , c ℓ,3 ) ∈ C 3 et C est l'espace des nombre complexes. Par ailleurs,

  des méthodes et analyse des résultats pour m ≤ n/2 (approximation de Fourier).

  Figure 1.2: Edge system of a Weaire-Phelan foam. a) the bounded periodic unit; b) the composed foam.

Figure 1

 1 Figure 1.3: ESRF synchrotron in Grenoble with: 1 the linear accelerator, 2 the booster synchrotron and 3 the storage ring [Facility (ESRF), 2009].

  Figure 1.4: Radiation propagation and corresponding diffraction pattern on different distances.

  Figure 1.5: Scheme of the near field interference of a sample, showing where interference occurs.

Figure 1 . 6 :

 16 Figure 1.6: Application of the Paganin filter on a slide of a reconstructed image from the dataset 110. The resolution is 1.1 ➭m. (a) reconstructed Laplacian image; (b) modified using Paganin filter.

Figure 2

 2 Figure 2.1: Volume rendering of the tomographic data set "A003" showing an open foam 512 3 out of 2048 3 pixels with a size of 1.1 ➭m. This volume was acquired by means of phase-contrast ➭CT using synchrotron radiation.

Figure 2

 2 Figure 2.2: Cross sections of a foam cell; the black point in the upper right corner shows the relative size of the particle which passes through the medium. Cross section of a Weaire-Phelan foam with marked mean cell radius and pore width (a). Real cross section at a resolution of 3.5 ➭m (b).

•

  Figure 2.3: Volume rendering of a tomographic data set "0019" showing a fiber fleece.

  Figure 3.2: A sequence of outer Jordan discretization with growing pixel number and increasing values of the integral of torsion. The tangent direction of the underlying line segment is (1, -1, 1).

  Figure 3.3: Transfer functions of filters: partial sum (PS), Cesàro's mean (CM) and binomial filter (B) of a discrete curve with n points and a smoothing parameter of m.

  Figure 3.5: Estimates of the third derivative of the cosine function for a blurring with σ = 0.001 and a sampling width ∆t = 0.19: true (T), cubic spline approximation (SA) with p = 5, Fourier approximation (FA) with m = 9 and application of the discretization of the differential-geometric formulas (DD). Notice that curve (T) is overlapped by curve (FA).

Figure 3 Figure 3 . 7 :Figure 3 . 8 :

 33738 Figure 3.6: Estimation of the mean error δ for the third derivative estimated by the Fourier approximation based on the partial sum for different smoothing parameters m. The mean error depends on the step width ∆t and the variance σ 2 , where the σ axis is in a logarithmic scale. Each image point stands for a specific combination of the step width ∆t and the variance σ 2 of Gaussian distributed position error ǫ. The mean error δ is encoded by colors according to the color scale (right). In the black area the mean error δ exceeds 1.

  Figure 3.9: Estimation of the mean error δ for the third derivative estimated with the discretization of the differential-geometric formulas for different filter masks. The mean error depends on the step width ∆t and the variance σ 2 , where the σ axis is in a logarithmic scale. Each image point stands for a specific combination of the step width ∆t and the variance σ 2 of Gaussian distributed position error ǫ. The mean error δ is encoded by colors according to the color scale (right). In the black area the mean error δ exceeds 1.

Figure 4

 4 Figure 4.1: The curves of the parametric functions f 1 (t) with r = 0.8 and c = 0.1 (a) and f 2 (t) with r = 1 and c = 1 (b) defined in Eq. (4.1) and Eq. (4.2), respectively.

Figure 4

 4 Figure 4.2: Estimates of mean local curvature κ (t) (a) and mean local torsion τ (t) (b) of curve f 1 (t). Calculated from the partial sum s 4 (t) (PS) and Cesàro's mean σ 4 (t) (CM) of 1 024 instances of the helix with the parameters r = 256 and c = 0.5. The blue lines (T) depict the theoretical values.

δ4 = 7 Figure 4

 74 Figure 4.3: Estimates of mean local curvature κ (t) (a) and mean local torsion τ (t) (b) of the curve f 1 (t) under reflection condition. The estimates are based on the partial sum s 8 (t) (PS) and Cesàro's mean σ 8 (t) (CM) of 1 024 instances of the helix with the parameters r = 256 and c = 0.5. The blue lines (T) depict the theoretical values.

Figure 4

 4 Figure 4.4: Estimates of mean local curvature κ (t) (a) and mean local torsion τ (t) (b) of the curve f 1 (t) with the restriction of the Fourier approximation to f * (t). The estimates are based on the partial sum s 4 (t) (PS) and the Cesàro's mean σ 8 (t) (CM) of 1 024 instances of the helix with the parameters r = 256 and c = 0.5. The blue lines (T) depict the theoretical values.

  An overview about the numerical values of the mean errors δm of curvature and torsion estimation for different low pass filter and different boundary treatments are given in Tables 4.1 and 4.2, respectively. Table 4.1: Numerical values of the mean error δm in units of 1 • 10 -3 of curvature estimation for different low pass filters and boundary treatments. Tested with 1 024 instances of IUR rotations of the parametric function f 1 with r = 256 and c = 0.5. Boundary treatment Partial sum Cesàro's mean Binomial filter None (m = 4)

Figure 4

 4 Figure 4.5: Mean of estimates of local torsion τ (t) of curve f 2 (t). The estimates are based on the partial sum s m (t) (PS) and the Cesàro's mean σ m (t) (CM) of 1 024 instances of f 2 (t) with the parameters r = 32, c = 0.5 and the smoothing parameter m = 8 (a) and m = 4 (b). The blue lines (T) depict the theoretical values.

  Table 4.3: Numerical values of the mean error δm of torsion estimation from 1 024 instances of IUR rotations of the curve f 2 (t) in units of 10 -3 . partial sum with m = 4 Cesàro's mean with m = 8 r c = 0.5 c = 1.0 c = 1.5 c = 0.5 c = 1.0 c = 1.5

Figure 4

 4 Figure 4.6: Power spectrum of the discretization error for (a) f 1 with r = 8 and c = 0.5 and (b) f 2 with r = 32 and c = 2.5. The spectra are estimated from 1 024 realizations of the random function ǫ(t).

Figs

Figure 4 . 7 :

 47 Figure 4.7: The relative error δ 1 of the curvature estimation from helices f 1 (t) depending on the true κ, the true τ , the step ∆ς between the sample points and the variance σ 2 of Gaussian distributed displacements. The relative error is encoded by colors according to the color scale on the right figure side. In the black area the relative error δ 1 exceeds 20 %.

Figure 4 . 8 :

 48 Figure 4.8: The relative error δ 1 of the torsion estimation from helices f 1 (t) depending on the true κ, the true τ , the step ∆ς between the sample points and the variance σ 2 of Gaussian distributed displacements. The relative error is encoded by colors according to the color scale on the right figure side. In the black area the relative error δ 1 exceeds 20 %.

Figure 4

 4 Figure 4.9: Torsion estimation error for parametric function f 2 (t) with r = 2.5, c = 1.0, ∆ς = 1 and σ = 0.001: (a) κτ -diagram of the relative estimation error δ 1 for the local torsion τ of helices f 1 (t) superimposed by the projection (grey curve) of the parametric function f 2 (t); (b) the relative estimation error δ 1 taken from the κτ -diagram along the projection of f 2 (t)(AP) and δ 2 directly estimated from points x k on f 2 (t) (DE) .

Figure 4 . 10 :

 410 Figure 4.10: The influence of the choice of the smoothing parameters m and p of the Fourier respectively spline approximation on the mean error δ of estimating torsion; (a) relationship between m and δ; (b) relationship between p and δ. The results are obtained from 1 024 discretizations of IUR rotations of f 2 with r = 32, c = 1.0 (C1), r = 32, c = 0.5 (C2) and r = 128, c = 1 (C3) for t ∈ [0, 2π].

Fig

  • 10 -3 and the mean estimates of the different methods are τ = 15.85 • 10 -3 ± 4.21 • 10 -3 for the Fourier approximation with m = 4, τ = 15.13 • 10 -3 ± 37.26 • 10 -3 for the discretization of the differential-geometric formulas, τ = 11.52 • 10 -3 ± 19.52 • 10 -3 for the spline approximation with p = 4 • 10 -5 ,

Figure 4

 4 Figure 4.12: Comparison of torsion estimation of f 2 by the spline approximation (SA), the Fourier approximation (FA) and the discretization of the differential-geometric formulas (DD). The example curves f 2 have the radius r = 32 and the slope parameter c = 1. The mean values are computed for 20 000 independent instances of f 2 with about 750 sample point each. The confidence bands include 90% of the values. The spline smoothing factor is p = 4 • 10 -5 and the Fourier approximation uses m = 4. The true torsion value is depicted by T

  Figure 5.1: Sketch of the setup for time-resolved ➭CT with the interior setup of the image acquisition device. The distance between the sample and X-ray source was about 150 m.

  Figure 5.2: Setup to detect particle trajectories. (1) Erlenmeyer flask with suspension, (2) Heidolph Unimax 1010DT shaker, (3) Pharmacia Peristalic P-1 syringe pump, (4) porous media inside a cylindrical cartridge.

Figure 5

 5 Figure 5.3: Histograms of the spherical contact distribution (a) and of the spherical granulometry distribution (b) of the solid matter and of pore space, respectively, of the partially open foam "A003". Both estimations were done from a subvolume consisting of 1024 × 1024 × 1080 pixels.

Figure 5

 5 Figure 5.4: Percolation probability of the pore space of the partially open foam "A003" of a subvolume consisting of 1024 × 1024 × 1080 pixels of the 3D image.

Figure 5

 5 Figure 5.6: 3D image examples of the partially open foam with some fastest particle paths computed with the adapted FMM (see Chapter 3.2.1). This figure represents a subimage of 1448 × 1448 × 1024 pixels out of the full image "A005" consisting of 2560 × 2560 × 2160 pixels. The pixel size is 1.4 ➭m.

Fig. 5 .

 5 Fig.5.6 shows a visualization of the investigated foam. Examples of fast particle paths, which are used to calculate the curvature and the torsion, are shown in red. The fast particle paths are obtained by an adapted version of the FMM, which is applied to a single 3D image as described in Section 3.2.1. The principal motion direction is chosen to be the z-direction. The extracted paths correspond to an outer Jordan discretization of a continuous curve. Therefore, the Fourier approximation with a partial sum as low pass filter and a smoothing parameter of m = 4 are selected to estimate the curvature and torsion of the particle paths.

Figure 5 . 7 :

 57 Figure 5.7: Curvature (a) and torsion (b) distribution of fast particle paths through the foam. These results were obtained with m = 4 and applying the partial sum in the Fourier approximation.

Figure 5 . 8 :

 58 Figure 5.8: Sequence of subimages with a slow moving particle marked with a red arrow, where the subimages are presented as a volume rendering. In this example the marked particle is slowed down with increasing time. The pixel size of the 3D image is 1.1 ➭m.

Figure 5

 5 Figure 5.9: (a) Histogram of distribution of the absolute values |τ | of the local torsions τ of particle paths extracted from sequences of 3D images obtained by in-situ time-resolved ➭CT (low energy particles). (b) The κτ -diagram of the relative estimation error δ h for the local torsion τ of helices for ∆ς = 8.3 and σ = 0.000 5. This diagram is superimposed by estimates (κ i , τi ) of curvature and torsion values extracted from the sequences of 3D images (small gray dots).

Figure 5 . 10 :

 510 Figure 5.10: Subset image (700 × 700 × 750 pixels) of fiber fleece "053" with a pixel size of 0.65 ➭m. The shortest paths through the fiber system are in red.

Fig. 5 .Figure 5 . 11 :

 5511 Fig. 5.12a. The peak of the distribution is at the radius of 23.4 ➭m, which corresponds

Figure 5 Figure 5 . 13 :Figure 5

 55135 Figure 5.12: Histograms of the spherical contact distribution (a) and of the spherical granulometry distribution (b) of the fibers and of pore space, respectively, of the fiber fleece. Both estimations were done from a subvolume consisting of 1024 × 1024 × 1080 pixels of the 3D image.

Figure 5

 5 Figure 5.15: Cross sections through fibers: (a) the real cross section, (b) the theoretical cross section which is convex almost everywhere, (c) proposed ideal cross section, which is concave almost everywhere.

Figure 5 .

 5 Figure 5.16: Realization of a random system of overlapping fibers with a core system of random helices characterized by λ = 2.5 • 10 -3 ➭m -2 , r = 192 ➭m, c = -0.4 and ̺ = 5 ➭m. The volume consists of 512 3 pixels of size 1 ➭m and the fastest paths are given in red.

Table 4 .

 4 2: Numerical values of the mean error δm in units of 1 • 10 -3 of torsion estimation for different low pass filters and boundary treatments. Tested with 1 024 instances of IUR rotations of the parametric function f 1 with r = 256 and c = 0.5.

	Boundary treatment Partial sum Cesàro's mean Binomial filter
	None (m = 4)	13.46	6.88	11.63
	Reflection (m = 8)	2.71	0.87	3.79
	Restriction on f * (t)	1.12	0.66	4.32
	(m = 4)
	Restriction on f * (t)
	(m = 8)

Table 5 .

 5 1: List of investigated open porous media, which will be characterized in this chapter.

	ID	Material	Kind of	Resolution Total volume [pix]
			media
	053 019 A003 polyurethane Foam nylon Fibers nylon Fibers A005 polyurethane Foam	0.65 ➭m 0.65 ➭m 3.5 ➭m 1.4 ➭m	2560 × 2560 × 2160 2560 × 2560 × 2160 2560 × 2560 × 2160 2560 × 2560 × 2160

Table 5 .

 5 2: List of investigated sequences of 3D images with suspension flow. Each sequence consists of 20 3D images with a temporal resolution of 0.5 s.

	ID Material	Kind of media Resolution Total volume [pix]
	110 polyurethane Foam 112 polyurethane Foam 113 polyurethane Foam 126 polyurethane Foam	1.1 ➭m 1.1 ➭m 1.1 ➭m 1.1 ➭m	2016 × 2016 × 2016 2016 × 2016 × 2016 2016 × 2016 × 2016 2016 × 2016 × 2016

  Dry Partially Open Foams in AirThese 3D images were taken 2012 for a preliminary investigation. The same foam is taken for different resolutions, but the center position remains the same. Partially Open Foams with Si Particle Suspension These sequences of 3D images are obtained with time-resolves ➭CT. The pump was sometimes turned off in order to slow down the particles. In the last sequence of 3D images, the time resolution was decreased to 17 ➭s. The 3D images with ID 55 until 59 were used for adjustment.The 3D images acquired with laminography represent a thin mat of glass fibers.

							A.1. Dry Fleeces with Cells A.4. Laminography Data A.5. Fibers of Filter Cake
	2.185 13.617 -0.515 13.617 -3.215 13.617 2.185 13.117 -0.515 13.117 2.185 13.8 ID Z-pos Y-pos Discarded Description 5 x 5 x 5 x 5 x 5 x 5 x Value 0 0 dry foam 0 0 foam with water Si suspension 0 0 x foam with water Si suspension 7.53 -4.49 4.53 -4.49 3.03 -4.49 3.03 -4.49 7.53 -5.99 pump turned off mediumA_a 62.21➦ same medium as mediumA_b mediumA_b 62.21➦ slice 1567 correspond to slice A.2 Parameter [mm] [mm] 7.53 -4.49 506 in mediumA_a
		0 7.53	-3.215 13.117 Image dimension 0 x -5.99 x	5	x 2048x2048x2048 foam with water Si suspension
		0 4.53	-3.215 12.117 Total angle 0 -5.99		5	x foam in alcohol without Si particle x 180➦
		0 4.53	-0.515 12.117 Angles number 0 -0.49	x	5	x new sample 1000
		0 4.53	2.185 12.117 shutter time per angle 5 0 x 1.51	x foam with alcohol Si suspension 1 ms
		0 6.53	0.185 18.867 Distance sample scintillator 300 mm 15 x -3 1.51	x
		5 8.53	0.185 14.5 Energy -3 2.01	x	15	x 19 keV
		5 10.53	0.185 13.367 0 2.01	x	15	x new sample prepared
		1 12.53	0.185 12.117 0 2.01	x	15	x	x
		1 14.53	0.185 12.817 0 2.01		15	x pump turned off	x
		2 14.5	0.185 0 4.01	13.95 ID		15 Resolution [➭m] x pump turned off
		4.2 10.53	3.435 13.364 A0003 -1.2 -5.99	x	15	x new sample 3.5	x
		4.2 5.53	3.435 1.2 -5.99	13.93 A0004		15	x new sample prepared 0.7 new sample	x
		4.2 7.53	2.81 13.364 A0005 1.2 -5.99	x x	15	1.4	x
		4.2 9.53	2.81 1.2 -5.99	13.93 A0006		15	0.35	x
		3.1 5.53	0.185 -5.25 0	14.5		15	x pump turned off
		3.15 5.3	1.5 2.75 0	14.5		15	x
	1.5 1.5 13.364 19.93 -0.303 28.06 -0.303 15.064 -0.303 21.564 -5.303 21.564 -3.303 21.564 -2.903 15.064 -2.903 28.06 -1.603 28.06 -1.603 21.564 4.65 2.26 6.25 0 1 0.5 0.25 0.25 0.25 0 0 0 4 4 0 0 A.3 Parameter 0.65 5.28 1.78 1.78 1.78 5.28 9.03 3.53 9.53 5.53 5.53 5.53 3.53 3.53 0 0 Resolution 6.53 -2.25 4.53 2.25 A.4 Laminography Data x x	15 15 15 15 15 15 15 15 15 15 15	x x pump turned off pump turned off new sample time resolution 17 s Value 1.1 ➭m pump turned off	x x x
		-1.603 0.897 Image dimension 15.06 24.76 4.53 2.25 Total angle 3.53 -5.75	x	15 10	2016 × 2016 × 2016 180➦
		0.897 Angles number 18.76 3.03 4.75	x	10	500 pump turned off
		0.897 Shutter time per angle 12.76 Distance sample scintillator 100 mm 10 1 ms Energy 26.3 keV Video frames 20 3.03 4.75 3.03 -5.25 Parameter Value x 3.03 Resolution 0.33 ➭m -5.25 3.03 -5.25 Image dimensions 2560 × 2560 × 2300
		Time resolution 3.03 -4.49	x		0.5 s new sample
	ID	3.03	-4.49 Laminography angle Description
	4.53 mFilterVlies_C 63.32➦ -4.49

  RésuméL'extraction des cellules cancéreuses d'un fluide corporel est une procédure importante lors d'un diagnostic clinique et d'une thérapie. En particulier, lorsque la technique de séparation est basée sur la chromatographie cellulaire, il est important de disposer de connaissances précises sur les capacités de liaison des cellules cibles avec le milieu poreux. Pour cette raison, des expériences utilisant la tomodensitométrie à résolution temporelle ont été conçues et réalisées à l'Installation Européenne de Rayonnement Synchrotron. Les distributions des courbures et des torsions des trajectoires de cellules situées dans suspension s'écoulant à travers un milieu poreux sont des informations précieuses pour caractériser l'efficacité des procédés chromatographiques. Cependant, le calcul de la torsion est un défi car étant basé sur des dérivées d'ordre supérieur qui sont très sensibles au bruit de discrétisation. Cette thèse présente deux nouvelles méthodes d'estimation des courbures et des torsions de trajectoires de particules données respectivement sous la forme de points discrets connectés ou non connectés. La première méthode est basée sur une approche dite d'approximation de Fourier. Des études de cas ont mis en lumière une diminution de l'erreur d'estimation des torsions d'au moins 65 % par rapport à la méthode de référence d'approximation par les splines. Par ailleurs, le paramètre de lissage de l'approximation de Fourier peut rester constant pour une large plage de résolutions latérales et pour différentes valeurs de courbures et de torsion. La méthode dite d'approximation de Fourier n'étant pas applicable à des courbes échantillonnées avec un pas variable, une deuxième méthode basée sur la discrétisation des formules géométriques différentielles (DDGF) a été développée. L'approximation par les splines et la DDGF conduisent à des erreurs moyennes similaires. Cependant, le masque filtrant reste inchangé pour le DDGF, alors que le paramètre de lissage de l'approximation par les splines doit être adapté à la forme ainsi qu'au pas d'échantillonnage de la courbe.

	Parameter	Value
	Resolution	0.65 ➭m
	Image dimension Total angle	3586 × 3586 × 2160 360➦
	Angles number	6000
	Shutter time per angle	100 ms
	Distance sample scintillator 30 mm
	Energy		19 keV
	Stitching	activated
	ID	Y-pos [mm] Z-pos [mm] Color ID
	0015	-5.815	14.117	pink gold
	0016	-5.815	23.117	pink gold
	0017	-5.815	3.117	pink gold
	0017b	0.685	3.117	light blue
	0019	0.685	13.117	light blue
	0020	0.685	23.117	light blue
	0021	1.685	23.117	dark blue
	0022	1.685	13	dark blue
	0023	1.685	3.117	dark blue
	0052	-0.815	7	pink gold
	0053	-0.815	11	pink gold
	0053b	-0.815	17	pink gold

Acknowledgement

First and foremost, I offer my deepest gratitude to my supervisors Prof. Christian Daul and Prof. Joachim Ohser for their constant support and extraordinary advice. Without their encouragement throughout my thesis, with their patience and knowledge, this thesis would not have been possible. The door to Prof. Ohser's office was always open whenever I ran into a trouble spot or had a question about my research and writing. The distance between Germany and France was at no point an obstacle, Prof. Daul helped to overcome borders at any time. I would like to thank Dr. Bernd Stanislawski (Merck KGaA, Germany) for many valuable hints on the chromatographic background of this work and for providing the partially open foam, the fiber fleeces, the cells and the silica gel suspension. Moreover, I'm grateful for the opportunity to prepare my samples and practice the laboratory handling of cells and cell suspensions at Merck. Special thanks go to Dr. Alexander Rack for his valuable support during the image acquisition of the porous media at beamline ID19 of the European Synchrotron Radiation Facility (ESRF), Grenoble. His advice on time-resolved microtomography was a great help for me. Furthermore, I'm grateful to Katharina Losch, Jenny Hon, André Liebscher, Jan Niedermeyer and Christopher Helbig for assisting me with the experiments made at the ESRF. I would like to thank my colleagues at the "Mathematics and Science Faculty" at the Hochschule Darmstadt for being wonderful, pleasant and helpful. Many thanks to Christine Pierson from CRAN for her help regarding all administrative questions I had. Without her, I would have been lost. I am very thankful to all the jury members and the invitees for agreeing to participate in my thesis defense. I would like to thank the German Federal Ministry of Education and Research (BMBF) for their funding of this project under the grant MNT/05M13RCA/AniS. I also wish to thank all members of the Anis project (BASF SE, Fraunhofer ITWM, HARTING KGaA, Hochschule Darmstadt, MANN+HUMMEL GmbH, Merck KGaA, Ruprecht-Karls-Universität Heidelberg, TU Kaiserslautern, Universität Ulm) for their support and for inspiring me with new ideas in our discussions. In particular, I wish to thank Dr. Katja Schladitz (Fraunhofer ITWM) for providing me Mavi for my Linux PC and the support for this program. Finally, I must express my very profound gratitude to my parents and to my girlfriend for providing me with unfailing support and continuous encouragement throughout my years of study and through the process of researching and writing this thesis.

Bibliography

Publications relating to the work described in Chapter 3 [START_REF] Blankenburg | Torsion of particle trajectories through pore space and its estimation using information on local pixel configurations[END_REF]. Torsion of particle trajectories through pore space and its estimation using information on local pixel configurations. In Pirard, E., editor, Proceedings of 14th International Congress for Stereology and Image Analysis (ICSIA), page WA04, Liège. ISS Blankenburg, C., Daul, C., and Ohser, J. (2016a). Estimating torsion of digital curves using 3D image analysis. Image Analysis andStereology, 35(2):81-91 Publication relating to the work described in this chapter [START_REF] Blankenburg | Parameter free torsion estimation of curves in 3D images[END_REF]. Parameter free torsion estimation of curves in 3D images. In Proceedings 23rd IEEE International Conference on Image Processing (ICIP), pages 1081-1085, Phoenix, Arizona USA Blankenburg, C., [START_REF] Blankenburg | Torsion estimation of particle paths through porous media observed by in-situ timeresolved microtomography[END_REF]. Torsion estimation of particle paths through porous media observed by in-situ time-resolved microtomography. Journal of Microscopy, 266(2):141-152

Appendix A

3D Data Summary

The following sections describe the 3D data, which was taken at the beamline ID19 of the ESRF. The first table in each section summarizes important technical data used for the image acquisition. The second table gives an overview of the obtained data.

A.1 Dry Fleeces with Cells

The fleeces in this section were flown through by a suspension in controlled conditions and after a specific time a fixative was subsequently admixed until the cells were fixed and shape preserved. Afterward, the whole sample was dehumidified to guaranty a high contrast between fixed cells and surrounding air.