Nicola Viganò

Wolfgang Ludwig

Kees Joost Batenburg

Alexandre Tanguy

Simon Hallais

Alexandre Dimanov

Michel Bornert

Joost Kees

Batenburg

Discrete Representation of Local Orientation in Grains using Diffraction Contrast Tomography

Keywords: from plane normals to grains, tomography, optimization, 1, 3, 2 Sparsity and suitable sparse representations, 1, 3, 3 Tomography from diffraction of distorted volumes, 1, 3, 4 Advanced projection matrix modelling, 1, 3, 5 Reconstruction algorithms instances, 1, 4 Other existing approaches to the grain reconstruction, 1, 4, 1 Adaptive near-field HEDM forward modelling, 1, 4, 2 GrainSweeper and GrainSpotter, 1, 4, 3 Discrete tomography and homogeneity constraints, 2 Reconstruction discarding ω information 3

I would like to thank:

First of all my supervisors, Wolfgang and Joost, who definitely were fantastic role models for me and my intellectual growth, and who had so much to teach me during these years we worked together. More specifically, thanks to Joost because of his great ideas and pragmatism, that heavily contributed to shape my thesis. And many thanks to Wolfgang, because he was constantly and consistently guiding me throughout my entire PhD and at the same time he was able to put up with me and give me his friendship, regardless of all my bad jokes and oddities, to the point that Stefan used to tell a joke about "him almost adopting me as his seventh son".

All the graintracking crew members, like Stefan, Yoann, Andy, Lorenzo, Péter and Laura, for the very useful interaction at work, but mostly for the great time spent together both during and outside the working hours. Moreover, I should also mention Marco, my office mate for almost three years, for the very nice environment in the office and the many scientific and non-scientific discussions.

The Antwerp and Amsterdam crews

"Talk is cheap. Show me the code."

Linus Torvalds "Do or do not, there is no try." Yoda "La sconfitta è un'eleganza per chi si arrende in partenza." (Losing is an elegant way out for those who surrender without even fighting) Subsonica

Linus Torvalds is known to not be an easy person to deal with, and surely he had arguments with many people over the last two decades. Nonetheless, he is someone to be get inspiration from, and he surely has some wisdom to teach with his famous sentences. Among the many he said, Linus' sentence quoted on this page is my all times favourite, because it embodies my spirit in life. Only who gets his hands dirty and does the work that he envisions, proves that he deserves to talk about it. Moreover, in life there are only two results that matter: success and failure. We would like to never fail, and it doesn't matter how hard we tried, if that happens, there is no excuse that can change or alleviate that. However, this very fact is sometimes hiding also those cases where the will to succeed was simply not strong enough. Only who puts his maximum effort in what he does has a chance to achieve something, and maybe one day to succeed. I still have in Italy, who seem to still remember about me, and the new friends I made living in Grenoble and traveling around Europe. They are probably too many to be called out one by one, but they know who they are, and without them, it would have been impossible to make it to the end of the PhD.

Lucia, who has been there for me even when she wasn't, who loves me and unconditionally supports me, and is able to stand me, no matter how annoying I can be sometimes. 6 included in chapters 4 and 5 are about the application of the model to experimental data, with a validation against electron back-scatter diffraction (EBSD) measurements. Moreover, the problems analyzed in chapters 4 and 5 represent both milestones in the evolution of DCT towards the ability to analyze more complicated and interesting scientific cases. Chapter 4, more specifically, analyses the case of a NaCl sample, where the grains have the tendency to form sub-grains, making the indexing and reconstruction of its poly-crystalline structure more challenging than the samples affected by smooth orientation gradients. The results from this article show an excellent agreement of the local orientations with the EBSD map, and some of the reconstructions were even performed from raw-data, without any type of segmentation of the diffraction spots. While the single grain reconstructions using the six-dimensional framework (6D-DCT) show its ability to improve the shape of the reconstructed grains, compared to the previous 3D-DCT, it was also demonstrated in chapter 4 that the algorithm can reconstruct clusters of grains, by defining a bounding box in the six-dimensional space and using the raw-data associated to it. This fact suggests that not only we can surpass the single grain reconstructions, by being able to reconstruct entire regions inside a volume, with a clear benefit for highly textured samples, but more importantly that indexing is not a limiting factor any more for the reconstruction of deformed samples. The indexing of all the grains in a textured region or the indexing of all the sub-grains of a bigger grain is not needed anymore, because it is possible to reconstruct and successively index those regions by first running the reconstruction algorithm over their containing bounding box in the six-dimensional space. Chapter 5, instead, analyses the case of a 1% plastically deformed microstructure in Ti 4% Al alloy, giving rise to the formation of deformation twins. As it will be mentioned in chapter 5, twins are crystallographic domains that are very difficult to deal with, because they tend to be thin plate-like regions, which usually don't give rise to high signal intensities, and when they do, their shape is blurred by other non idealities of the acquisition system, like the point spread function (PSF) of the detector itself. The six-dimensional framework was able to improve the reconstruction of twinned grains, compared to the previous three-dimensional implementation, by introducing the possibility to reconstruct the twins together with their parent grain, and other technical improvements, like the possibility to estimate the structure factor of the analyzed materials, and to introduce a model for the PSF of the detector in the reconstruction problem.

Chapter 6 provides an outlook on the possible future extensions and applications of the developed reconstruction framework and the remaining limitations of the technique. When first published, the results on synthetic data from the third chapter highlighted some key features of the proposed framework, and showed that it was in principle possible to extend DCT to the reconstruction of moderately deformed materials, but it was unclear whether it could work in practice. The following chapters instead confirm that the proposed framework is viable for reconstructing moderately deformed materials, and that in conjunction with other techniques, it could also overcome the limitations imposed by the grain indexing, and be applied to more challenging textured materials.

Abstract (Français)

Contexte Cette thèse de doctorat introduit un modèle et un algorithme sixdimensions pour la reconstruction des orientations cristallines locales dans les matériaux polycristallins. Le modèle s'applique actuellement aux données obtenues avec un rayonnement synchrotron (faisceau parallèle et monochromatique), mais il est également possible d'envisager des extensions aux instruments et sources de laboratoire (polychromatique et divergent). Le travail présenté est principalement une extension de la technique connue sous le nom de "Diffraction Contrast Tomography" (DCT) qui permet la reconstruction de la forme et de l'orientation cristalline des grains dans des matériaux polycristallins (avec certaines restrictions concernant la taille et le nombre total de grains ainsi que la mosaïcité intragranulaire). La DCT utilise un rayonnement X dur (haute énergie) et monochromatique, associé à un dispositif expérimental simple pour l'acquisition des données : comme en tomographie X classique (par absorption), l'échantillon est placé sur une platine tournante et un détecteur bidimensionnel est positionné juste après l'échantillon dans l'axe du faisceau. Quand les grains entrent en condition de diffraction, pour les différents angles de rotation ω, ils produisent des tâches de diffraction sur le détecteur. La forme et l'orientation des grains sont ensuite déduites par une chaîne d'analyse qui passe par la segmentation des tâches, suivie par l'identification des paires de Friedel. Chaque paire de Friedel fournit un vecteur diffraction et une recherche systématique parmi cette liste de vecteurs permet alors de déduire les positions et orientations des grains présents dans l'échantillon (indexation). D'autre part, en DCT on suppose toujours qu'il n'y a que peu ou pas de désorientation intragranulaire et qu'il est donc possible de considérer les tâches de diffraction comme des projections géométriques du volume des grains. La forme tridimensionnelle de chaque grain peut donc être reconstruite par des méthodes algébriques. En effet, en l'absence de désorientation intragranulaire, les tâches de diffraction ne s'étalent que peu dans l'espace réciproque, et leur projection (2D) dans la direction ω est directement accessible. Or cette approche simple, dans laquelle chaque grain est décrit par son orientation moyenne, perd sa validité quand le matériau s'éloigne des caractéristiques idéales, notamment si la désorientation intragranulaire dépasse quelques dixièmes de degrés. À chaque différente (sous-)orientation correspond alors une géométrie 10 de diffraction légèrement différente. Par conséquent, pour un grain déformé avec des orientations locales très proches, on obtient des tâches déformées qui s'étalent et se superposent sur plusieurs images consécutives de l'acquisition, produisant ainsi des volumes tri-dimensionnels nommés diffraction blobs. L'information de la désorientation locale à l'intérieur d'un grain est alors encodée dans la distribution 3D d'intensité de l'ensemble des tâches associées à un grain. Dans ce mémoire nous proposons un modèle et des algorithmes qui permettent de remonter à cette information par la résolution du problème inverse associé. La méthodologie développée permet de s'intéresser à des micro-structures contenant des grains présentant une désorientation intragranulaire qui peut alors atteindre plusieurs degrés.

Formulation six-dimensionnelle Pour décrire l'orientation locale des grains nous avons choisi un modèle qui incorpore à la fois l'information spatiale et l'information d'orientation. Comme suggéré dans [24], la description la plus complète de ce type de problème nécessite une représentation dans un espace à six-dimensions qui donne la possibilité d'associer une distribution d'orientation (ODF) à chaque position dans l'échantillon. L'espace six-dimensionnel proposé est simplement le produit dyadique d'un espace réel à trois dimensions (de type euclidien), et d'un espace des orientations à trois dimensions représenté par un sous-volume de l'espace de Rodrigues (et englobant l'ODF du grain en question). Le paramétrage de Rodrigues a beaucoup d'avantages, et il permet notamment une approximation euclidienne pour des petites régions de l'espace. Pour travailler avec cet espace à six-dimensions nous avons choisi de faire un échantillonnage régulier sous forme de voxels cubiques dans les deux sous-espaces associés à la position et à l'orientation. Plus précisément : à chaque orientation échantillonnée dans l'espace des orientations, nous associons un volume dans l'espace réel ou, réciproquement : à chaque voxel dans l'espace réel, nous associons un volume (ODF) dans l'espace des orientations. La quantité scalaire associée à un voxel dans l'espace six-dimensionnelle correspond au "pouvoir diffractant" de cet élément. En utilisant cette représentation discrète de l'espace à six-dimensions, il est possible de représenter mathématiquement la reconstruction comme un problème inverse : on résout un système linéaire. Par contre, pour atteindre une résolution suffisamment fine, l'échantillonnage dans l'espace des orientations peut contenir de 10 3 à 10 4 degrés de liberté. On se retrouve alors avec une énorme quantité d'inconnues sans pour autant avoir augmenté le nombre de mesures. Dans la pratique, il est plus convenable de transformer la résolution d'un tel système linéaire fortement sous-déterminé par un problème de minimisation d'une fonctionnelle incluant d'une part la norme l 2 du résidu du système linéaire (data fidelity term) et d'une autre part un nouveau terme qui permet de régulariser le problème en sélectionnant parmi l'infinité des solutions celle qui minimise la fonctionnelle combinée. Deux exemples pour ce deuxième terme, utilisés dans ce travail, sont la norme l 1 du vecteur des inconnues (le volume de reconstruction en six-dimensions), et la "variation totale" (TV) du volume tri-dimensionnelle obtenue par la somme de toutes les orientations pour chaque voxel dans l'espace réel. L'ajout de ces contraintes permet de sélectionner des solutions physiquement significatives. Notons toutefois que le choix de la fonctionnelle peut dépendre de caractéristiques micro-structurales des échantillons. Par exemple dans le cas de la norme l 1 on va favoriser des reconstructions avec un caractère plus parcimonial (peu d'orientations actives par position dans l'espace réel), alors que dans le cas de la variation totale du volume réel on va favoriser des reconstructions avec une distribution d'intensité homogène. Pour la solution numérique de ce type de problèmes, il est important de choisir un algorithme adapté à ce type de fonctionnelles. Pour la solution des problèmes contenant la norme l 1 , on a plusieurs possibilités. Un des algorithmes le plus connus est FISTA [2], dont une implémentation est présentée dans le deuxième chapitre. Une autre possibilité est d'utiliser des algorithmes de type Chambolle-Pock [7], qui offrent la flexibilité d'adapter, et donc d'optimiser, l'algorithme à la fonctionnelle choisie.

Structure du manuscrit Le premier chapitre de ce mémoire sert comme introduction et rappel d'un certain nombre de concepts nécessaires à la compréhension des chapitres suivants, présentés sous forme de collection d'articles scientifiques. On peut y retrouver une introduction à la description mathématique des rotations et orientations cristallographiques dans les matériaux, les détails essentiels de la DCT, le formalisme utilisé en optimisation mathématique et appliqué aux problèmes de reconstruction tomographique, ainsi qu'une revue succincte d'autres approches expérimentales et algorithmiques pour la caractérisation de l'orientation cristalline dans les matériaux polycristallins. Dans les chapitres suivants le modèle six-dimensionnel est introduit et d'abord évalué sur des données synthétiques (chapitres 2 et 3) avant d'être appliqué à des données réelles provenant de micro-structures de plus en plus exigeantes (chapitres 4 et 5). Chacun de ces chapitres comprend une introduction succincte qui décrit le contenu et les résultats principaux de l'article associé. Dans ces parties on peut également trouver des résultats et des détails supplémentaires qui n'ont pas pu être inclus dans les versions soumises aux journaux scientifiques pour des raisons de restrictions de longueur et/ou nombre de figures.

Chapitre 2: Le modèle initial Le deuxième chapitre présente la première formulation du modèle qui n'était pas encore complète, et qui essaie de résoudre le problème de l'orientation locale à partir de projections bi-dimensionnelles (résultant de la sommation des tâches de diffraction). Comme l'étalement des tâches de diffraction en ω encode l'information de la désorientation dans la direction de l'axe de rotation, il s'avère difficile de résoudre par grains cette composante spécifique de la déformation. Même si cette première tentative n'a donc pas pu aboutir aux résultats recherchés, elle a démontré qu'il était possible d'améliorer les reconstructions obtenus avec le modèle tri-dimensionnel, utilisé jusqu'à présent. Nous avons donc pu valider l'idée de base que des algorithmes visant à résoudre (régulariser) des problèmes fortement sous-déterminés en exploitant le fait que la solution recherchée doit avoir un caractère parcimonial, étaient applicables au problème inverse de la reconstruction d'orientations locales à partir des tâches de diffraction.

Chapitre 3: Le modèle complet Dans cet article, publié dans Journal of Applied Crystallography, nous avons amélioré le modèle introduit dans le chapitre précédent et nous démontrons sur des données synthétiques qu'il est effectivement possible de reconstruire l'orientation cristalline locale dans un grain. Le modèle utilisé est similaire à celui introduit dans le chapitre précédent mais en y ajoutant l'information de rotation ω de l'échantillon, dans la structure des tâches de diffraction des grains. En effet, utiliser l'information codée dans la rotation ω permet d'augmenter la quantité d'information (le nombre de mesures) d'un ordre de grandeur. Grâce à cette information supplémentaire, il devient possible de remonter à la composante de désorientation dans la direction de l'axe de rotation. Dans cet article nous avons légèrement changé le type de régularisation utilisé, en passant d'une minimisation de la norme l 1 de la transformé en ondelette de Haar du volume six-dimensionnel du chapitre précédent, à la minimisation directe de la norme l 1 du volume six-dimensionnel. Nous avons choisi ce nouveau type de fonctionnelle pour favoriser des reconstructions à caractère parcimonial, en accord avec la nature de la déformation de l'échantillon modèle (phantom) utilisé dans ce chapitre. Nous avons testé l'algorithme et le modèle dans plusieurs conditions, en essayant d'identifier et sonder ses limites en termes de ressources informatiques et de la résistance au bruit. Les résultats ont été très convaincants et encourageants.

Chapitre 4: Une première application à des données expérimentales et extension aux données brutes Le chapitre quatre propose une première application du modèle six-dimensionnel à des données réelles : il peut être considéré comme le résultat le plus important de ce travail de thèse. L'article inclus dans ce chapitre a été soumis à Physical Review Letters, et traite de la reconstruction d'un échantillon de halite (NaCl). La micro-structure de cet échantillon, caractérisé par la présence de sous-joints de grains à faible désorientation (de 1 à 5 degrés) rend l'indexation (i.e. le processus qui identifie le centre de masse et l'orientation des grains à partir de tâches de diffraction) beaucoup plus difficile que dans d'autres types de microstructure. Les différents sous-domaines peuvent selon les cas, générer des clichés de diffraction qui sont parfois séparés (plusieurs tâches segmentées), parfois superposées (une seule tâche connectée) sur le détecteur. De ce fait, certains sous-grains ne peuvent être trouvés par la procédure d'indexation DCT classique, et la reconstruction (assemblage) de l'ensemble de l'échantillon présente des régions vides. Comme nous le montrons dans ce chapitre, il est possible d'identifier et reconstruire ces zones manquantes en sélectionnant une région d'intérêt dans l'espace six-dimensionnel et en y associant la région correspondant dans l'espace détecteur. Notons que l'identification de la région d'intérêt est facilitée par la sous-structure mentionnée de l'échantillon qui présente des parties dans l'espace réel caractérisées par des orientations très proches dans l'espace des orientations. En utilisant les positions et orientations des sous-grains indexées il était possible de définir une région contenant toutes les orientations concernées, incluant aussi les orientations des sous-grains non indexées. La comparaison avec une cartographie EBSD acquise sur une surface latérale confirme bien que cette nouvelle approche permet de reconstruire et identifier des sous-grains manquants. Soulignons que dans ce cas on utilise les images brutes, sans segmentation ni indexation préalable pour procéder à la reconstruction de la région d'intérêt. Une autre différence par rapport à la formulation du chapitre précédent est l'utilisation d'une fonctionnelle modifiée visant à favoriser l'uniformité de l'intensité reconstruite. L'amas de sous-grains le plus grand reconstruit dans l'article joint à ce chapitre, si considéré comme un seul grain, avait des dimensions considérées prohibitives lors de la rédaction de l'article du chapitre précédent, et grâce à l'avancement des technologies des ordinateurs modernes, sa reconstruction est devenue possible un an après. En considérant que les technologies des ordinateurs sont en constante évolution et que le modèle proposé ici est capable de reconstruire la micro-structure dans des sous-régions confinées sans avoir recours à la phase de segmentation et d'indexation utilisée en DCT conventionnelle, on peut penser que dans un futur proche il sera possible de reconstruire des régions encore plus grandes et plus déformées, contournant ainsi certaines des limitations intrinsèques à la DCT avant l'introduction du modèle six-dimensionnel.

Chapitre 5: Reconstruction en présence de macles de déformation Le cinquième chapitre présente un article traitant de la reconstruction d'un échantillon d'alliage de titane contenant 4% d'aluminium, de structure hexagonale, caractérisé par une micro-structure et un type de déformation assez différent par rapport au cas traité dans le chapitre précédent. Les grains présentent des gradients d'orientation et dans certains grains favorablement orientés par rapport à l'axe de déformation macroscopique, on observe la formation de mâcles (i.e. domaines avec une relation d'orientation avec le grain parent) sous forme de fines plaques. En optimisant les paramètres d'acquisition, on montre qu'il est possible d'identifier ces mâcles à partir d'épaisseurs de l'ordre du micron. La reconstruction des mâcles pose des problèmes supplémentaires par rapport au cas analysé dans le chapitre précédent. Bien que le grain parent et ses mâcles soient caractérisés par des orientations moyennes qui sont liées par des opérations de symétrie, ces orientations ne se rapprochent jamais dans l'espace des orientations. La reconstruction simultanée nécessite alors la gestion de régions d'intérêt séparées, ce qui ne pose effectivement pas de problème majeur dans le cadre du modèle six-dimensionnel. Néanmoins, la reconstruction simultanée du parent et de la mâcle nous a obligé à baser la normalisation des intensités des différentes réflexions sur des bases physiques (tels que le facteur de structure et la géométrie de diffraction). De plus, les faibles épaisseurs et intensités des tâches liées aux mâcles nécessitent un traitement qui tient compte de la réponse impulsionnelle du détecteur. Moyennant ces améliorations, nous démontrons qu'il est possible d'aboutir à une reconstruction en bon accord avec la cartographie EBSD, acquise sur le même échantillon déformé à un pour-cent, présentant des désorientations intragranulaires de l'ordre du degré et des mâcles de déformation d'une épaisseur (sub)-micrométrique.

Conclusions et perspectives Le modèle six-dimensionnel et les algorithmes d'optimisation présentés dans ce mémoire transforment la technique de la tomographie par contraste de diffraction, auparavant limitée à des microstructures caractérisées par une faible désorientation intra-granulaire, en la rendant applicable à des matériaux présentant des microstructures plus complexes et/ou difficiles à reconstruire. La description de la microstructure par "grain" (décrit par son orientation moyenne) est remplacée par une description six dimensionnelle, i.e. une fonction de distribution d'orientation tri-dimensionnelle pour chaque voxel réel du volume reconstruit. Les résultats des chapitres 4 et 5 montrent qu'il est effectivement possible de reconstruire la déformation locale dans des échantillons modérément déformés dans le cadre de la tomographie par contraste de diffraction. Cette variante de la méthodologie de diffraction X tri-dimensionnelle (3DXRD) est basée sur une illumination en champ plein (2D), ce qui rend l'acquisition au moins un ordre de grandeur plus rapide comparée à d'autres variantes de cette méthodologie. Nous avons également démontré que le modèle six-dimensionnel peut reconstruire et indexer des grains, en utilisant les images brutes associés à une région d'intérêt dans l'espace à six-dimensions. Ce dernier résultat est intéressant dans le sens où il offre une possibilité d'éviter les étapes de segmentation et d'indexation, concepts connus pour être vulnérables en présence de désorientations intragranulaires importantes. A titre d'exemple nous mentionnons le cas de matériaux fortement texturés, où l'on peut trouver de nombreux grains avec des orientations très proches, difficiles a traiter dans le cadre de la DCT classique. Nous terminons par indiquer de possibles pistes de développement, qui consistent pour la plupart à l'adaptation du modèle six-dimensionnel à d'autres types et géométries d'acquisition :

Multirésolution Il semblerait intéressant d'associer des mesures "far-field" à chaque acquisition "near-field". Les mesures "far-field", où le détecteur se trouve à une distance plus grande mais possède des pixels plus larges, sont caractérisées par une résolution angulaire supérieure et une résolution spatial inférieure par rapport aux mesures de "near-field". Cela permettrait d'évaluer plus facilement l'ODF de chaque grain, et d'utiliser cette information comme contrainte supplémentaire dans la fonctionnelle utilisée lors de la reconstruction six-dimensionnelle. Vu le nombre de degrés de liberté réduit dans l'espace réel, il serait également envisageable, de considérer l'extension du modèle sur tout l'espace d'orientation, au moins pour le cas d'une illumination réduite (tranche par tranche, donc un problème 5D). Cette stratégie rendrait les étapes de segmentation et indexation obsolètes et rejoint les idées et efforts actuellement poursuivit par nos collègues Danois.

Sources de laboratoire Il a déjà été possible d'appliquer la DCT traditionnelle dans ce contexte, mais uniquement pour des matériaux modèles (absence de désorientation intragranulaire). Le récent développement d'un projecteur polychromatique compatible avec les algorithmes ici présentés [START_REF] Wim Van Aarle | An accurate projection model for diffraction image formation and inversion using a polychromatic cone beam[END_REF], devrait rendre possible l'application de la DCT six-dimensions à des échantillons modérément déformés.

Topo-tomographie Dans cette configuration un vecteur de diffraction est aligné avec l'axe de rotation, qui lui-même est incliné en position de Bragg, ce qui permet d'acquérir un grand nombre de projections pour le grain sélectionné -mais avec une résolution spatiale bien supérieure à la DCT [22,[START_REF] Simons | Dark-field X-ray microscopy for multiscale structural characterization[END_REF]. En utilisant le modèle proposé ici il sera désormais possible de s'intéresser à des grains déformés.

twikkeld en gevalideerd op synthetische diffractie-data van fantomen. Hoofdstuk 4 en 5 bevatten twee artikelen die zich nog in het review-proces bevinden. In deze hoofdstukken wordt het werk van hoofdstuk 2 en 3 uitgebreid naar het werken met experimentele diffractie-data, waarbij de techniek wordt gevalideerd met metingen verkregen uit electron back-scatter diffraction (EBSD). De problemen die worden geanalyseerd in hoofdstuk 4 en 5 kunnen worden gezien als mijlpalen bij het ontwikkelen van DCT naar een breder toepasbare techniek, die geschikt is voor het beantwoorden van nieuwe interessante en gecompliceerde wetenschappelijke vragen. Meer specifiek gaat het in hoofdstuk 4 om de analyse van een NaCl sample, waarbij de korrels de neiging vertonen om deel-korrels te vormen, waardoor de indexatie en reconstructie van de poly-kristallijne structuur uitdagender is dan bij samples die slechts gladde oriëntatie-veranderingen bevatten. De resultaten uit dit artikel laten een uitstekende overeenkomst zien met de resultaten verkregen uit de EBSD-techniek. Bovendien werden enkele reconstructies uitgevoerd op de ruwe meetdata, zonder eerst segmentatie van de diffractie-spots toe te passen. De reconstructies van individuele korrels met de zes-dimensionale techniek (6D-DCT) laten zien dat de techniek de vorm van de korrels beter reconstrueert dan met de eerdere 3D-DCT techniek. Daarnaast wordt in hoofdstuk 4 aangetoond dat de 6D-DCT techniek kan worden toegepast op clusters van korrels, waarbij een deelgebied in de zes-dimensionale ruimte wordt gedefinieerd en de ruwe meetdata hiermee wordt geassocieerd. Deze resultaten suggereren dat we niet alleen de techniek kunnen opschalen van individuele korrels naar de reconstructie van een volledig gebied (met duidelijke voordelen voor samples die veel textuur bevatten), maar, nog belangrijker, dat de indexatie-stap niet langer een beperkende factor vormt voor de reconstructie van gedeformeerde samples. De indexatie van alle korrels in een gebied dat textuur bevat, of de indexatie van alle sub-korrels in een grotere korrel, is niet langer nodig omdat het mogelijk wordt om eerst de reconstructie uit te voeren en pas daarna die gebieden te indexeren, door de reconstructies uit te voeren in de betreffende deelgebieden van de zes-dimensionale ruimte.

In hoofdstuk 5 wordt een casus uitgewerkt van een 1% plastisch gedeformeerde microstructuur in een Ti 4% Al legering, waardoor gedeformeerde tweelingstructuren ontstaan. Deze tweeling-structuren zijn kristallografische domeinen die uitdagend zijn om mee te werken, omdat ze bestaan uit dunnen plaat-achtige structuren die weinig signaal opleveren. Zelfs als er wel voldoende meetsignaal is, worden de metingen vaak uitgesmeerd door andere verstorende factoren, zoals de puntspreidingsfunctie (PSF) van de detector. Met de nieuwe zes-dimensionale techniek is het mogelijk om de reconstructie van de tweeling-structuren te verbeteren ten opzichte van de eerdere 3D-DCT techniek, omdat het mogelijk wordt om de twee deelkorrels samen met hun "ouder-korrel" te reconstrueren. Daarnaast zijn er andere technische verbeteringen zoals de schatting van de structuurfactor van de geanalyseerde materialen, en de introductie van een model voor de PSF in het reconstructieprobleem. Hoofdstuk 6 blikt vooruit op de mogelijke toekomstige uitbreidingen en toepassingen van het voorgestelde raamwerk voor zes-dimensionale reconstructie en be-handelt de beperkingen van de techniek. Toen de resultaten van hoofdstuk 3, met synthetische data, voor het eerst werden gepubliceerd, werd aangetoond dat het in principe mogelijk was om DCT uit te breiden naar een reconstructietechniek voor matig-gedeformeerde materialen. Het was toen echter onduidelijk of deze techniek ook in de praktijk goed zou werken. In de daaropvolgende hoofdstukken wordt bevestigd dat de techniek inderdaad bruikbaar is in de praktijk voor het reconstrueren van matig-gedeformeerde materialen. Bovendien wordt aangetoond dat als de 6D-DCT techniek wordt gecombineerd met andere technieken, de beperkingen van de indexatie-stap kunnen worden omzeild en toepassing op uitdagende materialen, met textuur, mogelijk worden.

22

CONTENTS

A Prediction of scattering ωs 175

List of Figures Orientations are a key concept in the modelling and description of physical properties and quantities like distortion and deformation. The representation of orientations and rotations is the same problem, because orientations are rotations with respect to a reference initial state [12]. The first ever parametrization of orientations was presented by Euler in his paper "Formulae Generales" [10] (1776), and it goes by the name of Eulerian angles. This representation is still used nowadays, especially in crystallography, where the three angles that define an orientation are also used as settings of a threecircle goniometer, and the reference orientation is simply the one for which these settings are (0, 0, 0). However, this first eulerian representation suffers from at least two non negligible disadvantages: the asymmetry of the parameters, since they are in no sense equivalent to each other, and the pathological behavior of the parametrization close to the origin, where all the problems arising from the non-commutativity of the rotations should disappear.

In the very next published paper, called "Nova Methodus", Euler described the first of these two disadvantages and presented a new parametrization based on the couple axis of rotation and associated angle of rotation [11]. This new representation is actually the base for a whole new class of parametrizations. If we now call the unitary axis of rotation n and the corresponding rotation angle φ, we would then have four parameters to describe rotations, with an extra equation that connects three of those parameters: the one enforcing the axis modulus equal to one. It is then an obvious choice to combine them together in one entity, by multiplying n by some function of φ like:

1. nφ 2. n tan 1 2 φ = R 3. n sin 1 2 φ = Q 25 CHAPTER 1. INTRODUCTION
Expression (1) is the simplest but also the least suitable for any kind of calculation, while the (2) is the representation given by Rodrigues in 1840, and R is nowadays called Rodrigues vector. Finally the (3) is the definition of the Quaternion parameter vector, named Q, which represents the vector part of a quaternion, while the scalar part would be Q 0 = cos 1 2 φ. One interesting difference between the mappings (2) and ( 3) is that the latter is finite, enclosed within the sphere of radius one in a four-dimensional space, while the former extends to infinity. However, for crystal lattices with more than triclinic symmetry, this is not a problem, because there is a certain degree of redundancy in the representation, that makes every remote point equivalent to much closer points to the origin. Finally, one common disadvantage of all the representations (1), ( 2) and ( 3), is that they are not homochoric, which means that for a completely random distribution of orientations we don't have equal probabilities in equal volume elements in different points of the orientation space defined by these orientations. In the end, this is not a big problem because they are sufficiently nearly homochoric close to the origin. Having introduced a few of the possible orientation parametrizations, we will now take a step back and focus a bit more in depth on another great tool know as Orientation Matrices which are, as we said before, nothing else than rotation matrices. While orientation matrices are not particularly useful for the purpose of orientation mapping, they are instead very precious when it comes to computations in the orientation space, and especially when dealing with high symmetry crystal lattices.

Rotation matrices

We know from elementary physics that while rotations about the same axis do commute, rotations about different axes don't. For instance, given a coordinates reference system XYZ, a rotation about the z -axis by 30 • , represented by R z (π/6), followed by a rotation along the same axis by 60

• is equal to R z (π/3)R z (π/6) = R z (π/2) = R z (π/6)R z (π/3
). But if instead we consider two 90 • rotations about the x -axis and the z -axis, R x (π/2)R z (π/2) will not be equal to R z (π/2)R x (π/2). To give a more quantitative demonstration of the previous statement, we will now work out quantitatively the manner in which rotations fail to commute when they are not about the same axis. Let's first recall that in three-dimensions we represent rotations using 3 × 3 real, orthogonal matrices. Taking a vector v, when we rotate it, the new vector will be called v , and the old and new components will be related by the orthogonal matrix R, in the form of: v = Rv, having RR T = R T R = 1, where the superscript T stands for a transpose of a matrix. Being R an orthogonal matrix, the l 2 -norm of v , defined as ||v || 2 will be equal to ||v|| 2 . If we now take a rotation about the axis z by an angle φ, we will now assume that the rotation applied to an object will affect only the object and leave the reference system unchanged, while the angle φ is positive when the rotation is 1.1. ORIENTATIONS 27 counterclockwise along the rotation axis direction. We then verify that:

R z (φ) =   cos φ -sin φ 0 sin φ cos φ 0 0 0 1   (1.1)
The rotations which follow the approach used to build equation 1.1 go sometimes under the name of "active rotations", while the ones where it is the reference system that rotates, are called "passive rotations". Let's now consider an infinitesimal rotation around the z -axis:

R z ( ) =   1 - 2 2 - 0 1 - 2 2 0 0 0 1   (1.2)
where the terms of order 3 and higher are considered negligible. We will also have:

R x ( ) =   1 0 0 0 1 - 2 2 - 0 1 - 2 2   (1.3)
when rotating along the x -axis and:

R y ( ) =   1 - 2 2 0 0 1 0 - 0 1 - 2 2   (1.4)
when rotating along the y-axis. If we now compute R x ( )R y ( ) and R y ( )R x ( ), we will obtain:

R x ( )R y ( ) =    1 - 2 2 0 2 1 - 2 2 - - 1 - 2 2    (1.5a) R y ( )R x ( ) =    1 - 2 2 2 0 1 - 2 2 - - 1 - 2 2    (1.5b)
From equations 1.5a and 1.5b, we obtain:

R x ( )R y ( ) -R y ( )R x ( ) =   0 -2 0 2 0 0 0 0 0   = R z ( 2 ) -1 (1.6)
from which we finally prove that rotation matrices about different axes do not commute, unless the rotation is so small that terms of order 2 are negligible.

CHAPTER 1. INTRODUCTION

Let's now consider the fact that while rotation matrices have 9 entries, the orthogonality constraint R T R = RR T = 1 corresponds to 6 linearly independent equations, being the matrix RR T symmetric, and so we have that we only need three independent parameters to define a rotation matrix. This means that for each point in one of the previously introduced representations of orientations we have one and only associated rotation matrix, and as a consequence one and only orientation matrix.

If we now consider that the effect of successive rotations can be obtained by the multiplication of the corresponding rotation matrices, we can observe that the set of all multiplication operations with these orthogonal matrices forms a group. This means that the following requirements for orthogonal matrices like orientation matrices are satisfied:

1. The product of two orthogonal matrices is still an orthogonal matrix:

(R 1 R 2 ) (R 1 R 2 ) T = R 1 R 2 R T 2 R T 1 = 1 (1.7)
2. The associative law holds:

R 1 (R 2 R 3 ) = (R 1 R 2 ) R 3 (1.8)
3. The identity matrix 1 -physically corresponding to no rotation -is defined by:

R1 = 1R = R (1.9)
and it is a member of the class of all orthogonal matrices.

4. The inverse matrix R -1 -physically corresponding to rotation in the opposite sense -is defined by:

RR -1 = R -1 R = 1 (1.10)
and it is also a member.

This group goes under the name of O(3), where O stands for orthogonal and 3 for three-dimensions [14].

From these properties, we deduce that we can easily compute the result of successive rotations by simple matrix multiplication of their rotation matrices. So, if a rotation R 2 is applied after an initial rotation R 1 , the result will be

R 12 = R 2 R 1 .
On the other hand, if we now consider two rotation matrices R 1 and R 2 as the orientation matrices that bring an object from the same initial orientation to two different orientations, the misorientation between the two final orientations is defined as

M 12 = R T 1 R 2 ,
and it is the rotation matrix that brings from the orientation associated with R 1 to the orientation associated with R 2 . We conclude by introducing a parallelism with the euclidean spaces, to better understand these key operations with the orientation matrices. We could indeed 1.1. ORIENTATIONS 29 see a similarity between the successive multiplications of rotation matrices and the sum of two elements in an euclidean space, with the important difference that the operation among rotation matrices is not commutative, and so, given

R 1 , R 2 ∈ O(3) and r 1 , r 2 ∈ R n we have that: r 1 + r 2 =⇒ R 2 R 1 and r 2 + r 1 =⇒ R 1 R 2 , but while r 1 + r 2 = r 2 + r 1 , R 2 R 1 = R 1 R 2 .
For what concerns misorientations instead, the parallelism looks stronger, because the subtraction operation is not commutative in the euclidean spaces, which brings us to:

r 1 - r 2 =⇒ R T 2 R 1 and r 2 -r 1 =⇒ R T 1 R 2 .

Rodrigues space vs. Euler space

We will now first look at the type of orientation space generated by the parametrization using Euler angles, and expose its shortcomings. By then going in more depth with the description of the Rodrigues type of parametrization, it will be clear why we will use it throughout the rest of the book. The Euler angles define three rotations, which perform the desired rotation when applied in the correct order. The most common convention is the so called ZXZ, which performs the three rotations, with associated angles φ 1 , ψ, and φ 2 , first around the z -axis by the angle φ 1 , then around the rotated x -axis by the angle ψ, and finally around the rotated z -axis by the angle φ 2 .

The first observation about the Euler-space is that it is periodic, with period 2π:

R (φ 1 + 2π, ψ + 2π, φ 2 + 2π) = R (φ 1 , ψ, φ 2 ) (1.11)
but also:

R (φ 1 + π, 2π -ψ, φ 2 + π) = R (φ 1 , ψ, φ 2 ) (1.12)
which defines a glide plane in ψ = π (a reflection plane, with a simultaneous displacement of π in both φ 1 and φ 2 ). The immediate result from this is that the determination of the Euler angles from the corresponding rotation matrix is not straightforward, because it depends on the considered range of angles that define. In fact, the asymmetric unit of the Euler space, which is the smallest volume of the space that contains all of the structural and symmetry information of the space itself, and which is in bijective correspondence with the rotation matrices, depends on the choice of the allowed ranges for the Euler angles. Another disadvantage of the parametrization given by the Euler angles is that different rotations about the same fixed axis lay on curved lines in the orientation space. These disadvantages sum with the other problems identified in section 1.1. The orientation space generated by the Rodrigues parametrization, following its definition, is not periodic and it has only one singularity for the rotation angle close to ±π. Moreover, all the rotations with the same axis of rotation lay on straight lines, which are called geodesics [12]. Now, the rotation matrix R associated with a given Rodrigues vector r = 30

CHAPTER 1. INTRODUCTION
n tan 1 2 φ can be obtained by the formula:

R = 1 1 + r • r (I (1 -r • r) + 2 (r ⊗ r + I × r)) (1.13)
where r ⊗ r is the outer product of the vector r with itself, and I × r is the so called tensor product of a tensor with a vector, which means that (I × r) ij = k l k ljI ik r l , where i kj is the permutation symbol or Levi-Civita symbol [19]. The Rodrigues vector associated to a rotation matrix can instead be computed using the following formula per each component i:

r i = - 3 j=1 3 k=1 ijk R jk 1 + 3 j=1 R jj (1.14)
These two transformations between the orientation matrices and the Rodrigues vectors will be denoted by the operators R = L (r) and r = L -1 (R). As for the rotation matrices (section 1.1.1), successive rotations are also represented in Rodrigues space. Given the vectors r 1 , r 2 ∈ R 3 , where R 3 is the Rodrigues space for rotations in three-dimensions, the result of successive rotations r f = (r 1 ; r 2 ) is given by the following formula:

r f = (r 1 ; r 2 ) = 1 1 -r 1 • r 2 (r 1 + r 2 + r 1 × r 2 ) (1.15)
which is the result of the application of a rotation represented by r 2 after the rotation represented by r 1 . An obvious observation is that rotations still cannot be simply added in the Rodrigues space, but also that for small rotations, where at least one of the two vectors r 1 or r 2 tends to zero, the dot product and the vector product in equation 1.15 tend to zero as well, making the result of the operation tend to the sum of the two vectors. This last observation suggests that the operations observed for the rotation matrices in section 1.1.1, find their equivalent in Rodrigues space, thanks to equation 1.15, so that the misorientation between two orientations represented again by the vectors r 1 and r 2 can be computed as m 12 = (r 2 ; -r 1 ). Another interesting observation is obtained if the vector r 2 is decomposed in the two parts r 2 = tan 1 2 φ 2 n2 = λn 2 , which results in the possibility for equation 1.13 to be rewritten like the following:

r f = (r 1 ; r 2 ) = r 1 + λ 1 -λr 1 • n2 • (r 1 • n2 • r 1 + n2 -r 1 • n2 ) (1.16)
which is obviously divided in two parts, one that is a scalar function of λ, and the other one that is a vector completely independent of λ. Moreover, this also confirms that the result of successive rotations about the same axis stay on that very same axis, which in the Euler space was only true for orientation trajectories passing through the origin.

If we now suppose to apply a translational coordinate change in the Rodrigues space, which corresponds to a rotation of the reference system in real space, 1.1. ORIENTATIONS and let the Rodrigues vector r z be the one corresponding to the shift, we would bring any vector r ∈ R 3 to a new position r ∈ R 3 with the simple formula r = (-r z ; r). This lets us conclude that straight lines transform into straight lines and that planes transform into planes. Finally, the angular distance of two points in Rodrigues space will be represented by the following scalar function:

φ = 2 tan -1 | (r 2 ; -r 1 ) | (1.17)

Crystals and crystal symmetries

We now apply the concepts introduced above to the study of crystallographic orientations. A crystal is mathematically defined as a Bravais lattice that has a certain basis associated to each point of the lattice, where the Bravais lattice is defined as an infinite set of points that follow the same geometric pattern through the whole space. By the definition of Bravais lattice, the points of a crystal lattice can be described in three dimensions by the formula l = n a a + n b b + n c c, where l gives the position of each point of the Bravais lattice that supports the crystal lattice, the numbers n a , n b and n c are integer numbers and the vectors a, b and c are the basis vectors of the lattice. Even if the Bravais lattices are supposed to be infinite, if we compare the size of the atoms that build a crystal lattice, which for our materials are the basis of the lattice, and the minimum size of volume elements allowed by the resolution used in the experiments, the approximation of infinite lattice can still be considered valid.

Then, the points of the basis, being located on a regular grid, for certain directions will appear to be organized on two-dimensional planes, called lattice planes. The vectors perpendicular to the lattice planes will then be called plane normals. To define a lattice plane, it is enough to identify three lattice points that belong to the plane, and a common way to represent a family of lattice planes all parallel to each others, is to use the so called Miller indexes. These indexes are simply the inverse of the three intercepts with the coordinates axes of the closest plane to the origin of the coordinates for a given family of lattice planes. Following the Miller indexes convention, sometimes lattice planes are identified by the three indexes (h, k, l), and the negative version of the indexes is usually represented with an "over-line": ( h, k, l). All the different permutations and negations of the three indexes will always generate lattice planes of the same kind, that will go under the name of an hkl-family.

Broadly speaking, crystalline objects can be seen divided in different categories, depending on the structure of the crystal: they are called mono-crystals when the crystal is the same everywhere and with the same orientation through-out the entire object, while if the crystal is divided in domains with different crystal orientations, they would be called poly-crystals, and the domains would be then called grains. At the same time, if the volume inside the grains has everywhere the same crystal orientation, the grains will be defined as undeformed, while if 32 CHAPTER 1. INTRODUCTION they are affected by some sort of deformation, it is possible to quantify this deformation by the means of orientation spread. Moreover, if the samples present multiple types of crystal structures, and so different types of lattices l, they will be identified as multi-phase materials. Given a reference coordinate system, and a crystal lattice that is described by three basis vectors, the orientation of the lattice (read its basis vectors) with respect to the reference coordinate system can be represented using the tools introduced in section 1.1.2. However, crystal lattices can also be characterized by the presence of specific properties called symmetries, which in turn will be influencing the properties of the used orientation space. Since both the determination of the crystal structure and the crystal orientation can be determined by means of diffraction measurements (i.e. X-rays, electrons, and neutrons), we will now briefly introduce the basic properties and implications of the different symmetry operations on the types of crystals and the diffraction patterns they will generate. Two objects are defined as congruent if there is a correspondence between each point of one object and one point in the other object, and if the distances among the points are the same for both objects. As a consequence, the corresponding angles will have the same absolute values. Congruences are then divided in two different categories: direct and opposite, where for the first, the corresponding angles have same signs, while for the second they have opposite signs [13].

The operations that belong to the first category are the following:

1. translation: an equal displacement applied to all the points in the same direction.

2. rotation: a rotation of the same angle applied to all the points about the same axis.

3. rototranslation (or screw movement): the combination of the two previous movements.

The operations that belong to the second category are the following:

1. inversion: all points are reflected with respect to the same point.

2. reflections: all points are reflected with respect to the same plane.

3. rotoinversion: the combination of a rotation and an inversion with respect to one point of the axis.

4. the combination of a reflection and a translation along the reflection plane; the plane is known as glide plane.

5. rotoreflection: the combination of a rotation and a reflection with respect to the plane perpendicular to the rotation axis.

If now we apply one of those operations to a crystal lattice, and the lattice stays the same, that operation will be called a symmetry operation for the given 1.1. ORIENTATIONS 33 lattice, while the points, axes and planes with respect to which the symmetry operations are performed, will be called symmetry elements. Not all these operations will be influencing our description of the crystallographic orientation, because the translation (1) in itself doesn't influence the crystallographic orientation. Moreover, rotoreflections (5) are always equivalent to other rotoinversion (3) operations, and reflections (2) are equivalent to rotoinversions of angle π. Even if rototranslations (3) or the presence of glide planes (4) are not distinguishable in orientation space from rotations (2) or reflections respectively (2), they will however produce systematic absences in their diffraction patterns. This means that only rotation and inversion symmetries could affect the representation of the orientation space. We will now define a point group as the set of symmetry operations that keep at least one point fixed, and for threedimensional crystals we can observe up to 32 point groups. However, diffraction is a phenomenon that is based on the distance among atoms, which is an inversion invariant property, and so in general it will not be possible to distinguish rotoinversions (3) from the corresponding simple rotations (2). This might become false for specific crystals, that break the inversion symmetry by having a chiral basis structure, but in this thesis we will assume to always work with materials that don't break it.

In the end, this means that out of the 32 point groups, many will appear to have the same properties, and they could be grouped under new bigger groups. These groups that can be observed by the means of diffraction experiments are called Laue classes, and their number is greatly reduced from 32 to 11. Moreover, if the basis of the crystal is taken into account, those 11 Laue classes can be further grouped into only 7 crystal systems which means that the orientation space will have 7 different structures, depending on the underlying crystal system symmetry. The 7 crystal systems are called: triclinic, monoclinic, orthorombic, tetragonal, trigonal, hexagonal, and cubic. All these different crystal systems, through their symmetry operations will make certain regions of the orientation space redundant, just like the way it was happening for half of the Euler space. This in turn means that, apart from triclinic, the asymmetric unit of the orientation space, will be smaller than the full space, and that the singularity for rotation angles of π will disappear, while it will still be an issue for some symmetry operations that involve rotations of π, which cannot be represented in Rodrigues space. On the other hand this also complicates the geometry of the orientation space, making the Rodrigues space not simply connected any more. The chosen asymmetric unit of the orientation space will be called fundamental zone.

By having lost the simple connectivity, computing distances and performing other operations in the space becomes more difficult. To understand how to perform those operations, we first have to understand how the mapping to the fundamental zone works. Let's call the set of symmetry operators S = {S i } where each S i is a different symmetry operator. Then, given an orientation R = L (r), we could get a symmetry equivalent orientation in the passive rotation convention by post-multiplying it by one of the symmetry operators, thus 34 CHAPTER 1. INTRODUCTION getting R i = S i R, while in the active rotation convention R i = RS i . We will now only use the passive rotation convention, because it simplifies the formulas and for the comparison of orientation from the very same crystal system, they are equivalent. At this point, if we define the fundamental zone as the asymmetric unit which is closer to the origin, the mapping of the given orientation R to the said fundamental zone, would simply be the R i * = S i * R for which:

i * = argmin i 2 tan -1 |L -1 (S i R) | (1.18)
We will then identify the mapping to the fundamental zone for both the orientation matrix and Rodrigues vector representation as R = F (R) and r = F (r) If no symmetry operation is taken into account (i.e. for the triclinic crystal system), the misorientation calculation on the orientation space would still be the same and so also equation 1.17 for the distance computation would stay the same, but for all the other crystal systems, the distance will instead have the following look:

φ * = 2 tan -1 |L -1 (S i * R 1 ) T (S j * R 2 ) | (1.19a) = 2 tan -1 |L -1 R T 1 S T i * S j * R 2 | (1.19b) (i * , j * ) = argmin i,j 2 tan -1 |L -1 (S i R 1 ) T (S j R 2 ) | (1.19c) = argmin i,j 2 tan -1 |L -1 R T 1 S T i S j R 2 | (1.19d)
which can be simply rewritten as:

φ * = 2 tan -1 |L -1 R T 1 S i * R 2 | (1.20a) i * = argmin i 2 tan -1 |L -1 R T 1 S i R 2 | (1.20b)
The misorientation will now be one of the intermediate results of equations 1.20 and be:

M 12 = R T 1 S i * R 2 (1.21a) i * = argmin i 2 tan -1 |L -1 R T 1 S i R 2 | (1.21b)

Coordinate systems

We already introduced the need for a reference coordinate system in section 1.1.3, for determining the orientation of a crystallographic domain. Normally only two types of coordinate systems could be used for fully determining the crystalline structure of a sample: a reference coordinate system attached to the sample and aligned with the main directions of the sample, and as many coordinate systems as the different crystallographic domains in the sample, each of them attached to one and only of the said grains, and aligned with its basis vectors. The first type of coordinate system will from now and on called the sample coordinate system, and the second will be called the crystal coordinate system instead. As introduced before, the relationship that allows us to change from the sample coordinate system to the crystal coordinate system and vice versa, is the orientation of each crystal domain.

Even if these two reference systems, coupled with the type and orientation information of each domain are enough for describing the crystallographic structure and orientation inside a sample, during an experiment the sample is likely to be rotating around one or more axes, thus a third reference system is need. This third type of coordinate system has to have the property of being fixed and so to be a perfect reference system for describing all the sample movements in space during the measurement, and to describe the equipment (i.e. detector) position. The axes of these coordinate systems will be called x l , y l and z l for the laboratory coordinates, x s , y s and z s for the sample coordinates, and x c , y c and z c for the crystal coordinates, while the lengths over these axes will be measured in millimeters [START_REF] Reischig | Advances in X-ray diffraction contrast tomography: flexibility in the setup geometry and application to multiphase materials[END_REF]. Finally there are two other coordinate systems that are worth mentioning for the understanding of the topic of this thesis: the detector coordinate system and the reconstruction coordinate system. The first coordinate system, as shown in figure 1.1, is attached to the detector, while the second coordinate system is attached and perfectly aligned to the sample coordinate system. The axes will be called x r , y r and z r for the reconstruction coordinates, and u and v for the detector coordinates. The length of these last five axes will be decided by the resolution of the experiment.

Orientation Distribution Functions and Inverse Pole Figures

We are now going to end this first section by briefly introducing some useful tools that are used in this thesis to represent the results and more broadly by 36 CHAPTER 1. INTRODUCTION the material scientists for representing some of the properties describing the crystalline nature of materials.

As mentioned in section 1.1.4, the orientation of a crystallographic domain is the relationship between the sample reference system and the crystal coordinate system of the given grain. If we now take all the grains in the volume, each of them will have a different orientation. Knowing the orientation of the grains in a sample, or even better, knowing the orientation of the crystal in every point in a sample is a piece of information that the material scientists can use to find the link between the microscopic properties and the observed macroscopic behavior of the studied materials. However, even if we achieve to gather such knowledge, we have to come up with a representation that is able to give some type of quantitative information.

If we now take dV as the fraction of volume in the sample that has a given orientation r within the volume element dr of the orientation space, and denote V by the total sample volume, then f (r) will be the orientation distribution function of the sample, and it will be defined as:

dV V = f (r) dr (1.22)
If we now imagine that all the grains inside of a sample would be perfect undeformed crystals, and so for each grain, all his volume elements, namely his voxels, would have exactly the same crystal orientation, the function f (r) would simply be a collection of Dirac deltas, centered on the specific orientation of the grains and their height would be proportional to the volume of the said grains. However, in the presence of sub-grain level misorientation, where the voxels of each grain have slightly different crystal orientations, those deltas become broader, while keeping the same total volume. Sometimes, materials can also show the presence of sub-grains, which are smaller portions of the grains that show small angle grain boundaries with the rest of the grain volume, and present themselves in the ODF as separate smaller volumes, which however are still close to the average orientation of the grain. For analysing these finer structures inside of grains it is then interesting to build the ODF of single grains instead of a more global sample ODF.

As the ODFs are an intrinsic three-dimensional object, it is not very suitable for visualization in two-dimensional figures. Since the direction of any threedimensional vector can be represented as a point on a unit sphere, in general, due to the symmetry of the crystals, if we take the points corresponding to the directions of a given set of plane normals for a given grain, they will be defining the orientation of the said grain on the unit sphere. A pole figure is the stereographic projection of these said points, by fixing the sphere coordinate system as the sample coordinate system. An inverse pole figure is instead the stereographic projection of the axes of the sample coordinate system, with respect to a unit sphere whose coordinate system is aligned with the crystal coordinate system. Following these principles, the inverse pole figure coloring (IPF color-coding) takes the stereographic projection of the orientation space fundamental zone,

DIFFRACTION CONTRAST TOMOGRAPHY

37

where a specific direction in the sample coordinates was chosen as the north pole, and it assigns a specific coloring to this two-dimensional region. Assuming then that all the voxels in a region of the sample can be approximately considered to have one and only orientation each, that region can be colored with the colors associated to their specific orientation.

Diffraction Contrast Tomography

Since the beginning of the last decade, big efforts have been put into the development of non-destructive three-dimensional imaging techniques for structural characterization of polycrystalline materials. These grain mapping techniques, which use X-ray diffraction as their mean for performing a description of a materials three-dimensional micro-structure, aim at resolving the local phase and crystal orientation, the so-called micro-texture. By being intrinsically nondestructive, they can also be used for studying the effects of temporal evolution of these micro-structures during processes like annealing and deformation. Different combinations of beams and detectors can be used. Beam shapes can be of three different types: zero-dimensional in the point beam configuration, one-dimensional in the line beam configuration, and two-dimensional in the configuration either known as full beam or extended beam. The beam can also be either mono-chromatic or poly-chromatic, while the detector can have high spatial resolution and be positioned a few millimeters away from the sample in the so called near-field configuration, or have low spatial resolution and be placed further away from the sample in the far-field configuration. Diffraction contrast tomography is a near-field monochromatic beam technique, usually using the extended beam approach [21], but interesting alternative configurations using a line beam and a different position of the detector have also been tested [START_REF] Reischig | Advances in X-ray diffraction contrast tomography: flexibility in the setup geometry and application to multiphase materials[END_REF]. It is however known to work under the restriction of dealing with undeformed or nearly undeformed materials, with intra-granular orientation spread of below one degree. Other techniques, discussed in section 1.4, have proven to work with more deformed materials, but they depend on the use of smaller beam profiles like the one-dimensional beam illumination, that requires longer scans for three-dimensional analyses, while DCT can be estimated to be at least an order of magnitude faster, because of its two-dimensional illumination.

This thesis is about improving the model behind DCT, to allow for the reconstruction of more deformed and challenging samples. But before diving in the technical aspects of the technique itself, we will first introduce the basics of the diffraction geometry.

Diffraction geometry

Having introduced the coordinate systems needed for describing all the quantities involved in the setup, we will now describe the diffraction geometry as a function of the three angles ω, θ and η. Later we will see how these parameters can contain important information about crystal orientations and orientation spreads inside the grains. First of all, given the interaction of the atomic planes in a crystal with an incoming beam, diffraction can be easily modelled as constructive interference, of the scattered beam from the atomic planes. For diffraction to occur at appreciable angles, the incoming beam (i.e. photons, electrons or neutrons), should have a characteristic wavelength that is comparable with the spacing among the points of the lattice. As it can be seen in figure 1.2, the effect of diffraction is similar to a reflection by the atomic planes of the incoming beam and the relationship that allows to predict when diffraction will happen is Bragg's law:

nλ = 2d sin θ (1.23)
where n is an integer number, λ is the incoming beam's characteristic wavelength, d is the distance between the atomic planes of the same type, and θ is the incidence angle between the direction of the incoming beam and the atomic plane. Taking now the experimental setup in figure 1.3, we suppose that the sample is sitting on a rotation stage, whose axis of rotation is perpendicular to the incoming beam direction, and the angle of rotation is called ω. We then usually align the laboratory coordinate systems to have the z -axis aligned with the rotation stage of the axis, the x -axis and y-axis of the laboratory system along the beam direction and perpendicular to the other axes respectively. The sample coordinate system will then be aligned to the laboratory coordinate system for ω = 0, and we will call Ω ω the rotation matrix that describes the orientation of the sample coordinate system with respect to the laboratory coordinate system. When meeting diffraction condition for one atomic plane, the angle between the plane normal of such atomic plane and the plane perpendicular to the incoming beam will be θ, and the angle between the incoming beam direction and the diffracted beam will be 2θ. Thanks to the perpendicularity between the incoming beam direction and the sample rotation axis, if a given atomic plane reaches diffraction condition at a specific rotation angle ω, it will also reach diffraction condition at ω + π, and the two reflections will form a Friedel pair.

Despite the fact of the lattice plane being exactly the same, the plane normal associated to the reflections in a Friedel pair will be opposite to each other: h and -h respectively, and the Miller indexes will also be (h, k, l) and ( h, k, l) respectively. Supposing that the grains have very little or even no intra-granular orientation spread, all the volume elements inside a grain will share the same orientation, and so the (θ, η, ω) angles of the diffracted beams for a given reflection will be the same for all of them. When such reflections hit the surface of the detector, they give rise to single images that are the geometric parallel beam projection of the grain volume, that are called called diffraction spots. On the other hand, in the presence of deformation inside the grains, each reflection hitting the detector will give rise to a stack of images very close in ω to each others, because the deformation will affect the (η, ω) angles, and the images of the said stack will no longer be geometrical projections, but will be deformed instead. These three-dimensional stacks of images generated by each reflection are called diffraction blobs. Since the goal of this thesis is the enabling of the grain structure reconstruction in presence of moderate levels of deformation, and the same reconstruction framework can always be used to reconstruct underformed grains, we will hereafter use the broader term "diffraction blob" also for the cases where "diffraction spot" would seem more appropriate. If now h is the normalized plane normal of an atomic plane in the crystal coordinate system, and g is the orientation matrix that defines the orientation of the crystal coordinate with respect to the sample coordinates, then we will call y = g -1 h = g T h the plane normal vector in the sample coordinate system. We can then compute Ĝl , the normalized plane normal in the laboratory coordinate system, using the following equation:

Ĝl = Ω ω g -1 h = cos θ   -tan θ -sin η cos η   (1.24)
where η is the angle between the rotation axis and the projection of the scattering vector on the plane perpendicular to the incoming beam direction, where the rotation axis resides.

If the orientation of the grain is known, it is possible to predict at which ωs 40 CHAPTER 1. INTRODUCTION it will get into Bragg condition for a given plane normal, using the following equation:

b • Ĝl = bT Ω ω g -1 h = ± sin θ (1.25)
where b is the normalized incoming beam direction vector, and the ± sign is due to the fact that each signed plane normal can give rise to two different reflections, as long as its negated version. In appendix A we find the four-fold solution to equation 1.25.

Orientations in DCT

In diffraction contrast tomography each grain represents a different entity with associated a grain center in real space and an average grain orientation in Rodrigues space. Following the discussion in section 1.2.1, we are now going to see how the coordinates (θ, η, ω) are connected to the grain orientation. Given a plane normal vector h in the grain coordinate system, we will have that y = g -1 h is the plane normal in the sample coordinate system. However, if the matrix g is not known, because the orientation of the said grain is not known, a minimum set of two couples (y, h) is in principle needed for retrieving the grain average R-vector, because, as mentioned in section 1.1.1, the matrix g is an orientation matrix, which in turn means that it has only 6 free parameters. Moreover, geometrically speaking, knowing the normal vector to an atomic plane, only fixes two degrees of freedom in determining that grain's orientation, because all the rotations around the said plane normal would be valid. We can retrieve y from the geometrical parameters of the recorded spots on the detector using the following equation:

y = Ω T ω cos θ   -tan θ -sin η cos η   (1.26)
and then each couple (y, h) will determine a line in Rodrigues space, that can be computed by the formulas:

r 0 = h × y 1 + h • y (1.27a) r = r 0 + t h + y 1 + h • y (1.27b)
where t is a free parameter. The crossing of two or more lines defined by equations 1.27, will be the average orientation of the related grain [START_REF] Henning Friis Poulsen | Three-Dimensional X-Ray Diffraction Microscopy[END_REF].

From the images to the Friedel pairs

The source of information gathered and used by DCT in reconstructing the grain structure of a sample, are the diffraction spots or blobs, that as we saw 1.2. DIFFRACTION CONTRAST TOMOGRAPHY before in section 1.2.1, are the result of each diffraction event happening in the grains. The goal of the initial steps of a DCT analysis is to isolate each diffraction blob from the background noise, and find the matching blobs that form a Friedel pair. These two steps are called in the DCT processing chain: blob segmentation and matching, respectively [START_REF] Reischig | Advances in X-ray diffraction contrast tomography: flexibility in the setup geometry and application to multiphase materials[END_REF].

The segmentation is based on a double threshold system, where a first intensity threshold finds the intensity peaks related to the spots, and then using a second threshold that is linearly proportional to the intensity of the peaks, the outline of the blob is determined. The segmentation process will then give rise to three-dimensional volumes in the (u, v, ω) coordinates, and the actual shape of the blob is identified by the use of morphological reconstruction, to avoid the partial inclusion of unconnected blobs, falling in the bounding box region of other blobs. What the segmentation cannot distinguish are the real blob overlaps, where different grains for two reflections would project to the same (u, v, ω) region. The size in ω of the diffraction blobs is called ω-spread.

The matching is based on the fact that blobs from a Friedel pair share many common features, and if the sample rotation axis is perpendicular to the incoming beam direction, the ωs at which they will be observed are connected by the law: ω b = ω a + π, where the two spots from a Friedel pair will be identified by a and b, being a the one with smaller ω. In particular, the blobs in a Friedel pair will have similar aspect ratio, identical ω-spread, and if the material doesn't break Friedel's law, also identical total intensity (except for small variations due to absorption and extinction). Finally, by having matched the spots belonging to Friedel pairs, it is possible to trace lines, in the sample coordinate system, from the centers of mass of one spot to the other one, and this line will pass through the center of mass of the grain that generated the spots. Using this piece of information to extract both η and θ associated to each pair and plotting a diagram of η vs. θ, allows us to optimize the geometry of the detector in a reliable and quantitative manner, making it possible to fit even small deviations from its ideal positioning and orientation for an actual dataset (i.e. no extra calibration measurement is needed).

Indexing: from plane normals to grains

The next step in the DCT processing route is called indexing, and it is the process that identifies grain orientations and positions, assigns labels (in our case unique IDs) to them. In DCT, the indexing is based on the Freidel pairs. The lines connecting the centers of mass of two spots from the same Friedel pair, can be used for determining the center of mass of the grains. In fact, the lines connecting the pairs belonging to one grain will all pass through the center of mass of the said grain. Moreover, by having identified these lines in real space, and being able to compute their (θ, η, ω) values, we can retrieve the associated y vectors to the pair, but also, apart from the pathological case with dissimilar h vector of the same norm (e.g. 333 and 511), the corresponding type of h vector from the angle θ. This means that the Friedel pairs will define both a line in real-space and a line in orientation space. More precisely, one has to take

CHAPTER 1. INTRODUCTION
into account all possible (symmetrically equivalent) variants of a given h vector (multiplicity) and draw the full set of lines for each pair. As a consequence, the indexing procedure is based on identifying the points of both the orientation and real spaces where these lines cross, and finding the ones that correspond [21] [26].

Once the "indexing" has produced these basic pieces of information about the grains, like grain center and grain average orientation, along statistical information about the size of the associated spots, the following step called forward simulation is responsible for searching through te list of all possible reflections hitting the detector and find the segmented blobs that possibly belong to these reflections, and in case also to reject the spots that have either too low intensities or that present overlaps with other grains. To determine whether to include and also enable these blobs for reconstruction, the forward simulation recently adopted a statistical analysis of the spots based on a two steps process, which first defines a likelihood for a given spot to belong to the grain, and then a measure of its goodness. These concepts are however restricted to "well behaved" micro-structures, and statistical measures are based on the validity of the initial set of pairs assigned during the previous indexing step.

Finally the grains, in case of undeformed materials, can be reconstructed using standard oblique-angles tomography algorithms (like SIRT), available from the ASTRA library [23]. The reconstruction, normally, has to take into account the fact that different hkl-families have different scattering intensities, due to the so called structure factors that differ from family to family, and the each spot's intensity may be influenced by other factors like incoming beam profile, absorption of the incoming beam and of the scattered beam, Lorentz factor and extinction. In practice (e.g. for metallic alloys), the structure factors may not be known a priori, and their exact determination is not easy for a near-field technique as DCT, however the total intensity of the blobs belonging to the same grain should never change because the projected volume is always the same. This means that it is possible to renormalize all the spots to any given number, and they would give rise to correct reconstructions. On the other side, being the intrinsic volume of each grain different, the expected gray levels for the regions belonging to the grains will not be related to each others.

From grains to grain-maps

While the raw detector images are usually corrected to account for possible distortion effects in the optics, the images are also usually affected by the so called detector point spread function, which is responsible for additional distortion and blurring of the images, and so of the diffraction spots as well. This means that the spots used in the reconstruction will be blurred to some degree, and so will be the reconstructions. As a consequence, the following step to the reconstruction will be the grain thresholding, and it is responsible for finding a decent threshold for segmenting the grain volumes and determine the individual grain shapes.

Once grains have been reconstructed and segmented individually, they have to

DIFFRACTION CONTRAST TOMOGRAPHY

43

be put together to form a grainmap. This process is called assembling, and it is very simple because the information on the center of mass of the grain is already known from the indexing step. Another complication arises from the fact that by blurring the reconstruction, the PSF also can affect the exact determination of the grain boundaries, and this results in conflicts during the assembling, where different grains would claim the same region at the boundaries. Up to recently the solution was to disable those regions for both grains, and then dilate the grain volumes until they touched each other, and letting them claim that region by assigning the contested voxels to the closest grain. Even if this last step looks quite rough and unphysical, it is actually following a simple principle: if some region of the grainmap is known to not be empty and if it is very close to one grain, it will be very likely that volume elements of the said region will belong to the closest grain, especially if the uncertainty over the grain boundaries is of the order of magnitude of the distance between the said voxels and the current boundary of the closest grain. However, recently the possibility for a rough estimation of the structure factors for the different families from the actual measured intensities was added (given in chapter 5), and in conjunction with a renormalization based on the Lorentz factor, that is responsible for broadening the blobs along ω by a factor close to 1/| sin η|, it was possible to rely on a grey level scale for all the grain reconstructions. This resulted in the possibility to resolve the boundary conflicts by choosing to assign a voxel in the region claimed by two or more grains, based on the gray level value that the voxel has in each of the reconstruction volumes.

Determining grain orientation spread from spots

It is possible to estimate a grain's orientation spread by looking at the spreads over the parameters η and ω of each detected diffraction blob. While the ωspread is obviously the size over ω of each blob, for what concerns the η-spread, a rough estimation could be given by the angle between the vectors located in the center of mass of the grain, and pointing to the extremes of the blob in the η coordinate. In practice this can be easily done, by dividing the estimation in two cases:

∆η = ∆v | sin η0| , if | sin η 0 | ≥ √ 2/2 ∆u | cos η0| , if | sin η 0 | < √ 2/2 (1.28)
where ∆u and ∆v are the sizes of the blobs in the u and v directions respectively, η 0 is the average η of the spot, and all the angles are supposed to be in radians.

Taking now equation 1.27, and computing the lines for the couplets: (-∆η/2, -∆ω/2), (∆η/2, -∆ω/2), (-∆η/2, ∆ω/2), and (∆η/2, ∆ω/2), these four lines in Rodrigues space, for small values of ∆η and ∆ω, will be almost parallel to each others, and we can approximate them to be the edges of a three-dimensional rod that defines a region of the orientation space where the ODF of the grain is confined in. By crossing all the rods from the available reflections would give an estimate of the three-dimensional orientation spread of the grain.

Mathematical optimization and tomography

In synchrotron imaging, computerized tomography is usually seen as a twodimensional technique, where filtered back projection is the most used tool for performing the inversion of the Radon transform. Indeed FBP works correctly for many projection geometries, but diffraction produces only a limited number of projections, that don't uniformly span the 2π range of rotation and the projections directions are usually inclined with the plane perpendicular to the rotation axis. If we also take into account the possibility for a projection geometry that depends on the orientation of the voxels, all these properties of the diffraction problem render the FBP algorithm unsuitable for its solution. Moreover, the FBP doesn't explicitly allow for introducing advanced priors in the reconstruction, that would enforce specific features of the solution. This is the reason why in this thesis we will cast tomography as an optimization problem, where the final reconstruction is the solution that minimizes a functional. The final goal will be to formulate more advanced models for the reconstruction of grains or even clusters of grains that are affected by non negligible intra-granular spread, and develop the algorithms able to solve the formulated minimization problems on a time scale that is compatible with the rest of a DCT analysis.

Tomography as optimization

In its discrete representation, standard three-dimensional tomography can be formulated as a linear algebra problem where the problem to be solved has the form:

Ax = b (1.29)
where x is the volume to be reconstructed, b represents the images recorded on the detector, and A is the so called projection matrix, which embodies the projection geometry of the tomographic problem. The meaning of equation 1.29 can be seen in figure 1.4, where each ray j corresponds to the line j in equation 1.29, and the coefficients A ji can have different definitions, and the two most popular are: the strip area overlap and the Joseph's interpolation. The first one assumes that the rays have a finite width, and the coefficients A ji are nothing else than the overlapping area between the ray j and the voxel i divided by the area of the voxels. Joseph's method instead samples the ray lines at a regular distance that is aligned with the voxel centers and the coefficients A ji are simply the distances between the said interpolation points and the center of voxel i [15]. For super-sampled formulations of Joseph's method, where intermediate points over the rays are also taken into account, we suggest the reading of previous work from the author of this thesis [START_REF] Viganò | Design and development of a tomographic library with physical corrections for quantitative analysis[END_REF].

Even if the formulation in equation 1.29 is very elegant, it is not very helpful when dealing with limited and noisy data. As the measurements are affected by noise, the vector b would contain slightly inconsistent entries, so a better solution to the problem would be to represent it in terms of a least square where from all the possible solutions we select the vector x * which minimizes the l 2 norm of the residual.

∇ ||Ax * -b|| 2 2 = 0 (1.31)
Unfortunately the matrix A doesn't usually have full rank, and so there is an infinite number of solutions for equations 1.29 and 1.30.

∇ ||Ax -b|| 2 2 = ∇ (Ax -b) T (Ax -b) (1.32a) = 2A T Ax -2A T b (1.32b) = 2A T (Ax -b) (1.32c)
In fact (1.30) implies that x * makes the equality in (1.31) true, and we see from (1.32) that (1.30) is equivalent to Ax = b. So every solution of Ax = b is also solution to (1.30), but the minimization formulation better deals with noise. Moreover, the functional in equation 1.30 is convex, and the space where it is defined, by being simply connected, is convex as well. The great advantage of convex problems defined over convex sets is that if they have a local solution, that solution will also be the global one [4]. One of the simplest but also most powerful methods for solving the functional in equation 1.30, is the SIRT algorithm (Simultaneous Iterative Reconstructive Technique) [16]. The SIRT algorithm is defined as an iterative first order gradient descend algorithm. An iterative algorithm usually takes n steps or CHAPTER 1. INTRODUCTION iterations for solving the given optimization problem, where for instance in the case of the problem in equation 1.30, it computes at each iteration i an estimate vector b that is the result of the forward projection of the current intermediate solution x i , so that b = Ax i , and computes the residual or distance from the real measurements b:

r (x i ) = b -b = b -Ax i (1.33)
Using the residual, the algorithm will compute an update for the solution x i , generating the current solution x i+1 for the next iteration. If the update is a simple back projection of the residual, with just an appropriate diagonal rescaling, the algorithm will belong to the class of gradient descend algorithms:

x i+1 = x i + P 2 A T P 1 (b -Ax i ) (1.34)
where P 1 and P 2 are diagonal rescaling matrices in the measurement space and the solution space respectively. Finally, the order of convergence will be a measure of how quickly the algorithm is able to approach the solution of a problem, and more precisely, an algorithm of order k is able to reduce the error of the current solution i to a precision of O 1/i k . So, as SIRT is a first order algorithm, the error of its current solution i will be of the order O (1/i).

Sparsity and suitable sparse representations

Having introduced in section 1.3.1 the formulation of tomographic reconstructions as minimization problems, we are now going to introduce the latest developments in signal processing theory, and give a brief introduction to methods that were developed over the years for solving the new functionals.

As mentioned before, equation 1.30 still suffers from undeterminacy due to the matrix A not being full rank. However, for quite some time Tikhonov regularizations have been used as the mathematical representation of a specific feature that the real solution was known to have and that could be enforced in the reconstruction. The introduction of this type of regularizations modifies the minimization functional like the following:

x * = argmin x ||Ax -b|| 2 2 + λ||Lx|| 2 2 (1.35)
where λ is a weight for the regularization term in the functional, and L is the preferred operator for enforcing a specific feature. Different operators have been used in the past, like I to help minimizing the energy of the reconstruction (read the square l 2 -norm of x: ||x|| 2 2 ), or some derivatives like D and D 2 to enforce some sort of smoothness or continuity in the reconstructed object. Another interesting property that could be enforced by different priors is sparsity. A sparse signal is mainly composed by zeros, and the preferred way to maximise the sparsity of a signal would be to minimize its l 0 -norm. By definition, the l 0 -norm is a semi-norm that if applied to a vector, it returns the number of non-zero components. However, using the l 0 -norm leads to a complex optimization problem which is not computationally tractable (in fact it is NP-Complete in some cases). As it can be seen from figure 1.5, for convex problems, an equivalent result can be obtained using the l 1 -norm, which is also convex, and the same conditions discussed in section 1.3.1 hold. It was demonstrated in [8] (1999-2001) that by choosing an overcomplete basis that would contain the features of the expected signal, and then performing a minimization over the l 1 -norm of the solution expressed in such representation, would lead to a robust recovery of the signals in the presence of noise. The key point of that publication was that the l 1 -norm would promote sparsity in the representation of the reconstructed signal, and let the reconstruction choose only the objects from the dictionary that would fit the signal better.

Later in [5] and [6] a whole new branch of the signal processing field, called Compressed Sensing, was started. The theory claims to lower the limits imposed by Shannon's theorem, for sparse signals, and for specific underdetermined problems, to be able to achieve perfect recovery. In CS the equivalent of the projection matrix A is called sensing matrix, and the most common are either random sensing matrices or randomly incomplete Fourier matrices. Nowadays, there is an ongoing work to incorporate the "tomographic sensing matrices" into the CS theory, so that new uncertainty principles similar to Shannon's theorem could be developed for tomographic problems. However, as discussed in [9], the properties of the tomographic projection matrix predict quite poor worst-case recovery results, while for the average-case more encouraging results have been found. Even if we cannot ensure any sort of perfect recovery at the level of underdeterminacy of normally treated tomographic problems, and so, at least for the moment, CS doesn't directly apply to tomography, it is possible to incorporate some of its ideas, for solving common problems within the two fields. A nice example is the possibility to improve the reconstruction quality of sparse objects, or find a suitable representation of the said objects that is sparse in the chosen domain, so that the l 1 -norm minimization techniques could still be 48

CHAPTER 1. INTRODUCTION used.
In fact, it is becoming increasingly popular to incorporate priors into the tomographic functional that would look like a Tikhonov regularization:

x * = argmin x ||Ax -b|| 2 2 + λ||Lx|| 1 (1.36)
where this time the l 1 -norm is used in place of the square l 2 -norm, and the L operator can be again the identity I for signals that are already sparse in their base representation, or the Haar wavelet transform H, for signals that are successfully sparsified by a wavelet representation, or finally by the generalization of the gradient operator D in more than one dimension:

||D (•) || 1 = ||∇ (•) || 1 = TV (•) (1.37)
where the l 1 -norm of the gradient is also known as the total variation operator. This operator is especially used for volumes that mostly exhibit flat regions with steep changes at their boundaries, because it enforces the sparsification of the derivative of the initial volume.

Another observation is that up to this point, we were only dealing with unconstrained optimizations, but intelligent constraints that preserve the convexity of the solution space, could both help to get better final reconstructions and also to speed up the process. For signals that only live in the non-negative orthant (generalization of the quadrant in n-dimensions), the non negativity constraint is a popular solution for constraining the solution. Indeed an example of reconstruction of a signal that is non-negative and that can be efficiently compressed (sparsified) in the wavelet domain, would look like the following Lasso (least absolute shrinkage and selection operator) formulation:

x * = argmin x ||Ax -b|| 2 2 + λ||Hx|| 1 (1.38a) subject to: x 0 (1.38b)
where the wavelet transform H is used in the regularization term, and is the element-wise "greater than" operator.

Tomography from diffraction of distorted volumes

Having introduced the standard three-dimensional tomography formalism and geometry in section 1.3.1, we are now going to extend the traditional tomography representation to the DCT geometry in presence of deformation in the grain volumes. The undeformed projection geometry depends indeed on the orientation of the grains, and it will be different for each grain because it will depend on six components: the three real space components x that determine the position of the grain, and the three orientation space components r that determine which reflections will fall on the detector and at which angle they will be coming out of the grain volume. For undeformed grains the position x i of the voxel i will change for every voxel, the orientation r will stay the same, and this translates in a rigid translation of the position on the detector from For what concerns deformed grains instead, the projection geometry will depend on each voxel's orientation. The easiest representation for this model, where we assume that each voxel has a local orientation, is a three-dimensional vector field that in each voxel stores the associated three-dimensional Rodrigues vector for the locally reconstructed orientation, and a fourth gray level value that indicates the scattering intensity coming from the said voxel. Each voxel i will then be a vector

x i ∈ V 4 = O 3 ⊗R
This model requires equation 1.29 to be modified into:

A(x)x = b (1.39)
where now the projection matrix A depends on the solution vector x, making equation 1.39 a non linear system of equations. As a consequence, formulations like the ones in equations 1.30 and 1.38, would be modified like the following:

x * = argmin x ||A(x)x -b|| 2 2
(1.40) and:

x * = argmin x ||A(x)x -b|| 2 2 + λ||Hx|| 1 (1.41a) subject to: x 0 (1.41b)
where, however, the non negativity and l 1 -minimization operations would only take place on the scattering intensity component. The three-dimensional vector field representation has one very big advantage due to its compactness, which results in it being extremely memory efficient, especially for big and heavily deformed volumes. However, it also has remarkable disadvantages that make it less than ideal for the mathematical properties, the algorithms and computer implementations, in practice. As stated before, the system of equations in 1.39 is not linear and so all the minimization formulations that derive from it are not convex any more. This means that any algorithm is not guaranteed to converge to the correct solution any more, and it might just end up finding local solutions. The forward projection matrix can be easily defined, but it will be much less computationally efficient than its simpler linear three-dimensional version, because a new projection geometry has to be computed at each projection for each voxel. For what concerns the back projection matrix instead, it will not be defined in any straight forward manner any more, because it should include some update mechanics for the solution vector x, that are not explicitly defined by the representation.

In this thesis we chose to adopt a more powerful representation, first presented in [24], that moves away from the single local average orientation representation for each voxel, and it introduces a local three-dimensional Rodrigues sub-space for the orientations of each voxel. We can then sample the orientation space with a fixed grid spacing, just like if it was a three-dimensional real-space. This translates in the description of the grain reconstruction solution space as a sixdimensional scalar field, that is the outer product of the cartesian real-space The solution vector in the six-dimensional formulation has two different interpretations that can make a significant difference when implemented. The most straight forward (figure 1.8a), which was just mentioned, would allow a three-dimensional orientation space for each real-space voxel, with the advantage of immediately having an ODF of each real-space voxel as the result of the reconstruction. However, computationally speaking, this implementation requires to compute the projection geometry for each voxel just in time each time that a voxel has to be forward projected or back projected. Even if the geometry could be stored for each voxel, occupying a big amount of memory, it wouldn't still allow to do vectorized and cache efficient operations, because the contiguity of the real-space part of the voxels would not be taken into account. The opposite implementation instead, which can be seen in figure 1.8b, associates a different real-space volume per each orientation that in turn has a fixed projection geometry. This results in the possibility to have fast projection and back projection routines that reuse the forward projection and back projection codes for conventional three-dimensional oblique-angle tomography, because each orientation-space voxel would just be a three-dimensional real-space volume, with a fixed geometry where only the offsets should be computed. Finally, even if the general look of equation 1.29 would stay the same, meaning that the minimization formulations would stay convex, the actual interpretation is the following:

   A 11 • • • A 1P . . . . . . . . . A M 1 • • • A M P       x 1 . . . x P    =    b 1 . . . b M    (1.42)
where each volume x i , which has a fixed orientation, would contribute to each diffraction blob b j on the detector, with i = 1, . . . , P , j = 1, . . . , M , P the total number of sampled orientations and M the total number of measured diffraction 52

CHAPTER 1. INTRODUCTION blobs.
More detailed discussions about the projection matrix and all its physical corrections and components, will be given in the next section.

Advanced projection matrix modelling

We will now present three different extensions and aspects of the modeling of the projection matrix that have to be taken into account when dealing with the analysis of real world use cases:

Blob scattering intensities Due to differences in structure factors, Lorentz factor, attenuation and incoming beam profile (as introduced in sections 1.2.4 and 1.2.5), each diffraction blob will have different scattering intensities. If these differences are known, we can introduce them in our projection model from equation 1.42 as the following:

Ax = C Ãx = C P i Ãi x i = b (1.43)
where the matrix C is a diagonal matrix, which predicts the scattering intensities, with Lorentz factor included for each of the reflections, and the matrix à assumes the same scattering intensity for all the blobs. As mentioned in section 1.2.4, we then want to renormalize the spots, and only keep the purely geometrical and non-invertible projection matrix à on the left side of the equation. We then obtain:

Ãx = N i Ãi x i = C -1 b = b (1.44)
where now b are the re-normalized diffraction blobs. Hereafter we will rename à into A and b into b, to simplify the notation.

Detector point spread function Another issue discussed in 1.2.5, was related to the detector PSF being responsible for some blurring of the reconstruction and loss in resolution. If the experimental PSF of the detector system used in a given experiment were to be determined, they could be modelled by the following equation:

P i F i A i x i = b (1.45)
where each matrix F i models the PSF for each orientation, and so it will differ for each incident beam direction over the detector. However it is a safe assumption to think that within a given diffraction blob, the difference in PSF due

MATHEMATICAL OPTIMIZATION AND TOMOGRAPHY

53

to the slightly different incident beam direction of each orientation will not be significant, and so we can rewrite equation 1.45, into the following:

F Ax = F P i A i x i = b (1.46)
where the matrix F is now the operator responsible for applying the PSF to the forward-projected images. The matrix F is usually not directly invertible, and the best approach is to simply treat it as an extension to the projection matrix A. This means that in the usual algorithmic iterations, where the backprojection matrix is usually modelled as the transpose of the forward projection matrix, extra care has to be taken to handle F . Generally speaking F is not symmetric, so the adjoint of F should be used, resulting in the algorithms trying to minimize:

A T F † F Ax = A T F † b (1.47)
where the superscript † stands for adjoint.

Orientation-space super-sampling Finally, as the projection matrix model is based on a discrete sampling of the six-dimensional space formed by the outer product of a three-dimensional real-space and a three-dimensional orientation space, the resolution in the sampling grid will play an important role for determining the accuracy of the tomographic model itself. While the real-space sampling resolution is given by the detector pixel size, which translates into the real-space voxel size, the orientation-space resolution is not fixed by the experiment itself. This means that if computer memory is a problem, or the sampled bounding box in orientation space is big, the orientation sampling could become not accurate anymore.

In real-space, a strip model of the projection rays could help with moderate dishomogeneity of the voxel size compared to the detector pixel sizes, especially when the pixels are larger than the voxels, because the strips would intercept all the voxels having an overlap with the strip. In this case, instead, Joseph's method wouldn't be able to associate some voxels to the related detector pixels, but various oversampling techniques exist to solve this problem, like the supersampling method used in [START_REF] Wim Van Aarle | Super-resolution for computed tomography based on discrete tomography[END_REF] for their super-resolution application.

In orientation-space each point determines a projection geometry in real-space, so it is not very obvious how to model a a projection geometry that is made by a solid three-dimensional interval of orientations.

It is in fact easier to start from the back-projection super-sampling technique presented in [START_REF] Wim Van Aarle | Super-resolution for computed tomography based on discrete tomography[END_REF], that is based on a sampling approach, and model a similar type of orientation super-sampling. The orientation space voxels, as shown in the inset of figure 1.9, instead of being represented by a single sampled point in the center of the orientation space volume forming the said voxel, they could be divided into S = s 3 sub-voxels, that would forward and back project the same associated real-space volume, but with slightly different orientations given by the centers of the new regions and projection coefficients of 1/S instead of 1.

The new projection matrix would then look like:

Ax = P i S j 1 S A ij x i = b (1.48)
where for a given volume i, the sub-matrices A i = {A i1 , A i2 , . . . , A iS } are the associated orientation-space sub-voxel sampled orientations.

For memory constrained configurations, the use of this type super-sampling can allow to greatly reduce artifacts due to otherwise poor sampling, because it virtually allows to increase by a factor s the resolution in orientation-space, without having an increase of memory occupancy by a factor S = s 3 . On the other hand, the increase in resolution is not as effective as if we were to use SP volumes, and the actual profile of the super-sampled projection matrix will be similar to a blurred non-super-sampled projection matrix.

Reconstruction algorithms instances

For solving the formulations in section 1.3.2, algorithms like SIRT are not enough any more, because even if they can be extended to incorporate non-negativity constraints or the six-dimensional geometry of section 1.3.3, they cannot naturally deal with extra terms in the functional.

The easiest choice to develop algorithms from these advanced functionals when they first appeared, was to use interior point algorithm instances that were tailor made for a given functional, just like in [5]. The advantage of those log-barrier Newton and primal-dual Newton methods is that a second order algorithm derives in a straightforward manner from the functional and constraints

OTHER EXISTING APPROACHES TO THE GRAIN RECONSTRUCTION55

formulations.

A log-barrier method can incorporate the non-negativity constraints into the functional with a weight that increases through the different steps of the newton method, and model them as logarithmic functions that tend to hard constraints when using increasingly larger scale factors. In the primal-dual methods, two problems are associated and solved together, where the functional of the primal problem is the one that we want to solve, while its constrains become the functional of the dual. Each iteration of the primal problem gives the update for the dual problem and vice versa.

Even if Newton algorithms have the great advantage of being second order algorithms, and to come straight from the functional formulation, constraints like non-negativity and sparsity, introduce increasingly larger and increasingly smaller eigenvalues respectively in the Hessian matrix used by the Newton algorithm. This causes the said Newton formulation to suffer from the poor conditioning of the problem during the later and fundamental stages of the iteration process.

For specific kinds of problems, another type of algorithm was later developed in 2009, and called FISTA [2]. This second order, gradient descent algorithm, has the advantage of having quadratic convergence, and the computational cost of almost a SIRT algorithm, while avoiding to build an increasingly badly conditioned Hessian matrix. However, it is only tailored to efficiently solve Lasso formulations 1.38, while the TV variations can be more tricky to be implemented [3].

In 2011, a new primal-dual second order framework for building algorithms directly from minimization functionals was developed by Chambolle and Pock in [7], and later polished for tomography problems in [START_REF] Sidky | Convex optimization problem prototyping for image reconstruction in computed tomography with the Chambolle-Pock algorithm[END_REF], which uses proximal mappings to solve the primal and dual problems, with a gradient descent based technique. In this thesis, the most prominent results were obtained using algorithm instances created from this powerful framework for solving each time the six-dimensional problem at hand.

Other existing approaches to the grain reconstruction

Polycrystalline materials are of great interest for material science and industries, so it comes to no surprise that various techniques exist already, and that various algorithms have been developed for those techniques. Most of the work done so far uses the three-dimensional vector field approach seen in section 1.3.3, while some full six-dimensional methods (i.e. covering the whole orientation-space) are under development.

In this section we are going to review the most prominent and inspiring methods in the field, where the first two sections describe some established methods, while the last section describes two similar methods for refining grain maps previously obtained with simpler methods.

56

CHAPTER 1. INTRODUCTION

Adaptive near-field HEDM forward modelling

The framework and algorithm developed in [START_REF] Suter | Forward modeling method for microstructure reconstruction using x-ray diffraction microscopy: Single-crystal verification[END_REF], and then improved in [20], is another near-field technique that has a completely different approach from DCT, and doesn't look for grains into the grain maps. It performs a completely independent voxel orientation optimization, of binarized detector images, and then assembles the output of this highly parallel optimization technique into a grain map. In [20], grain structures and some sort of multi-resolution technique can also be used to improve quality and computational times, but the core idea and final result are still based on the maximization of the following functional or confidence function for each voxel i at position x i :

C (q, x i ) = 1 N qual (x i ) h klN (q • G hkl , x i ) (1.49a) N (q • G hkl , x i ) = n χ [I s (j g , k g , z n , ω g ) , I e (j g , k g , z n , ω g )] (1.49b) χ (I s , I e ) = 1 if: I s ∧ I e 0 otherwise (1.49c)
where q ∈ SO( 3) is a given test orientation, C (q, x i ) is the actual confidence function. But before explaining the other objects in the set of equations 1.49, we make a step back and describe how an experiment of this kind is performed.

The setup is a regular near-field 3DXRD or near-field High Energy Diffraction Microscopy experiment, where, similarly to DCT, the sample is on a rotating stage, and a two-dimensional detector is positioned after the sample, with the direct beam hitting the lower central part of the detector, in the form of a slice (one-dimensional) beam. The laboratory reference system is identifying the direct beam direction as the z-axis, the sample rotation axis as the y-axis and the x-axis as the resulting one from a right-handed system. Given a rotation step size of δω, each rotation step l will be at ω l , and typical values would be δω = 1 and ω ∈ [0, π]. More than one scan are usually performed, with detector at positions L n , where n = 1, . . . , N , and N is the total number of scans. After a first step of background subtraction, which discriminates the noise from the diffraction spots, and it binarizes the images, the optimization is performed. Getting back to equation 1.49, given an orientation q, for a voxel i, it is possible to predict for every reflection of a plane hkl, its positions on the detector. The position of the reflection g on the detector at position L n will be identified by the tuple (j g , k g , z n , ω g ), where j g , k g are the column and the row on the detector matrix, z n is the coordinate position corresponding to the position L n of the detector, and ω g is the rotation angle for the said reflection. The values I s (j g , k g , z n , ω g ) and I e (j g , k g , z n , ω g ) will then be the experimentally observed signal and the predicted signal for a position predicted on the detector for the given reflection g. It is enough to say that, being the images binarized, the function χ only checks for the occurrence of intensity in the binarized images of the predicted peak positions.

If then we take N qual (x i ) as the number of possible peaks on the detector for

OTHER EXISTING APPROACHES TO THE GRAIN RECONSTRUCTION57

a given orientation, we see that the function N (q • G hkl , x i ) checks for the presence of all the reflections of a given hkl family and then computes the confidence function C (q, x i ) as a sort of completeness function for the given orientation q at the voxel i.

While this framework first published in [START_REF] Suter | Forward modeling method for microstructure reconstruction using x-ray diffraction microscopy: Single-crystal verification[END_REF] has a theoretical very sound base and has proven to reconstruct real data with deformations of up to 10 • across a single grain, it is still testing the complete fundamental zone of the orientation space for each voxel, while most of it is empty, and it doesn't take advantage of the fact that the interior of a grain is likely to have a similar orientation. In [20], it is indeed pointed out that the confidence function C (q, x i ) is sharply peaked in one small region of the orientation space, while it has just small fluctuations in the rest of the space, due to the noise. Even in the presence of pseudosymmetries, where some regions in the fundamental zone share a significant amount of reflections with the real orientation of the voxel, these regions show a sharply peaked profile, and most of the space remains practically empty. So the authors of [20] decided to act on three key points for reducing the impact of this morphology: (a) they would perform a multi-resolution analysis, were first the orientation space would be more coarsely tested, and then the resolution would become finer and finer in the later steps, when only certain sub-regions were identified as possible candidates, (b) in the earlier stages of the process fewer reflections would be tested, since most of the regions would show already lower confidence function values, and (c) they wouldn't test higher harmonics of a given family if the said family didn't match experimental data. Once few regions in orientation space are identified as final possible candidates, a local optimization based on a Monte Carlo approach would be performed to finally decide which region has the highest confidence function value.

Finally in [20], the possibility to seed the voxels with the orientation of their neighbours, and then try a local optimization from that point is also discussed.

As said earlier, this framework has been used at the Advanced Photon Source of the Argonne National Laboratory since its development, and has proven to work both on synthetic data and real data with significant deformation states. On the other hand, the great robustness of this method comes with the disadvantage of much slower three-dimensional acquisition compared to DCT, and a nonisotropic three-dimensional resolution.

GrainSweeper and GrainSpotter

Both GrainSweeper and GrainSpotter are algorithms developed by Søren Schmidt from DTU Physics, Denmark, and they work on near-field and far-field X-ray diffraction imaging data respectively [START_REF] Henning Friis Poulsen | Three-Dimensional X-Ray Diffraction Microscopy[END_REF], [START_REF] Schmidt | Direct observation of 3-D grain growth in Al-0.1{%} Mn[END_REF], [START_REF] Sø Ren Schmidt | GrainSpotter : a fast and robust polycrystalline indexing algorithm[END_REF]. GrainSweeper could also be applied to far-field data, but it would then work as an indexing algorithm, instead. On the other side, GrainSpotter is mainly thought as an indexing algorithm that can also fit grain center positions.

GrainSweeper was first introduced in [START_REF] Henning Friis Poulsen | Three-Dimensional X-Ray Diffraction Microscopy[END_REF], and reportedly used in [START_REF] Schmidt | Direct observation of 3-D grain growth in Al-0.1{%} Mn[END_REF] Given an experimental setup, very similar to the DCT setup, where a sample is placed on a rotating stage, and images on a inline two-dimensional detector are acquired at fixed angular steps, as the sample stage rotates. The images are then binarized, and the resulting spots are identified and their center-of-mass is recorded along with its spread. For each voxel in the volume, two hkl families are selected, one as primary, and one as secondary.

Since each of them will have spots appearing on circles (varying η, for a fixed θ) over the detector at certain ωs, the algorithm will try to find the orientation that satisfies the occurrence of spots on the primary and the secondary circles. This can easily done by taking all the primary spots, the spots falling under the primary circle at each ω, and for each of them searching at the ωs imposed by the crystallographic constraints, if secondary spots appear. For a couple made by a primary and a secondary spot, an orientation will be computed, and the full pattern will be compared against all the recorded spots. The orientation will then be considered a possible candidate only if a certain number of reflections matches with the recorded spots.

The limits and advantages of this approach are similar to the one seen in sections 1.4.1, because they both perform a voxel-wise optimization that doesn't take into account the intensities on the detector and possibly the crosstalk among different voxels.

GrainSpotter instead, as mentioned before, is an algorithm that operates on far-field data, where the pixel size of the detector is bigger than in the near-field case, but the detector is moved further away, so that the spatial information about the positions in the sample are almost lost, but the sensitivity on η is greatly increased, and the spots all lie on the rings at the specific θs of the hkl-families. By first assuming that all the scattering vectors come from the center of the sample, GrainSpotter performs an indexing step in the orientation space, by subdividing it in smaller volumes where the Rodrigues space is close to euclidean, and then identifies the crossing points of the lines (geodesics) generated by the (θ, η, ω) tuples of each identified spot. Each smaller volume is called "local Rodrigues space", and the identification of the vertexes where the geodesics cross, is performed by subdividing the local Rodrigues space in voxels, and counting how many geodesics pass by each voxel. Higher counts will mean higher probabilities for a correctly identified vertex. Only potential vertexes with a count of crossings higher than a given threshold will be kept. Once this step is done, and the grains are identified, GrainSpotter proceeds to fitting the grain orientation by minimizing the distance between the fitted average orientation and the geodesics. Optionally, it is also possible to carry out the same procedure for optimizing the grain center.

Discrete tomography and homogeneity constraints

In [1] and [18] two different types of algorithms are proposed for improving twodimensional grain maps, by incorporating morphology terms in the functional minimisation. Both methods use Gibbs priors, while in the case of [18] it also

OTHER EXISTING APPROACHES TO THE GRAIN RECONSTRUCTION59

an homogeneity prior that enforces specific types of connectivity. While the second approach is a greedy deterministic algorithm that takes advantage of concepts like tabu search to perform a global optimization, the first approach uses discrete tomography concepts to assign a grain label, and so an orientation to each voxel in a grain map. However, the main difference between the two is that the second approach makes explicit use in the homogeneity term of possible small misorientations among the voxels inside one grain, making it viable for deformed materials, while the first is only applicable to undeformed or nearly undeformed grain maps. The minimised functionals, look respectively like the following:

γ (f ) = exp [-β (α||P computed -P input || 1 + H 1 (f ))] (1.50a) γ (f , o) = exp [-β (α||P computed -P input || 1 + H 1 (f ) + H 2 (f , o))] (1.50b)
where f is a given grain map, o are the local orientations of the voxels, the term in the l 1 -norm is the data fidelity term on the detector, H 1 (f ) represents the Gibbs priors, and H 2 (f , o) is the homogeneity term. The Gibbs priors assume the existence of seven classes of 3 × 3 voxels features including the "grain interior" and the "grain boundaries", called G c , where c is the running index for counting the seven types. Each class has a predefined weight or preferences, called U c , so if N (G c , f ) is the number of features G c in the grain map f , the Gibbs functional will be:

H 1 (f ) = - 6 c=0 N (G c , f ) U c (1.51)
The term H 2 (f , o) is, as said before, the measure of homogeneity of the grain map, and it explicitly considers the local voxel average orientations o. If we define the following function:

Φ C (o) ij = exp - (d (o i , o j )) 2 2δ 2 (1.52)
as the likelihood of the orientation distance between two voxels i and j from the same grain, we can then define the homogeneity term as the following:

H 2 (f , o) = - n l=1   C∈C f ,l N N λ 1 Φ C (o) + C∈C f ,l N N N λ 2 Φ C (o)   (1.53)
where C f ,l N N is the set of vertically and horizontally neighbouring voxel pairs {i, j} in a grain l, and C f ,l N N N is the set of diagonally neighbouring voxel pairs {i, j} in a grain l. So far, both approaches have only been tested against synthetic data, and no real world study case is known.

Chapter 2

Reconstruction discarding ω information

In this chapter we are going to review the first publication in this thesis, that was about the first formulation of the model, where the core idea of a joint space between real-space and orientation-space was already in place, but some of its implications were not completely clear yet. This publication [START_REF] Viganò | Discrete Representation of Local Orientation in Grains using Diffraction Contrast Tomography[END_REF] was a conference proceeding presented in an international conference called "8th International Symposium on Image and Signal Processing and Analysis" in Trieste (2013). This article is relatively short because it was subject to a hard constraint of maximum six pages, but it is able to give a brief introduction to DCT, and then it effectively exposes the main ideas behind the six-dimensional representation and the advantages of using some additional constraints on top of the basic data fidelity term. After the description of the model it describes the used algorithm and its implementation, and it then introduces some interesting results obtained at the time. The algorithm used in this paper was still a modified FISTA instance that solved a rather interesting functional. Even if FISTA was not designed to solve functionals like the one in equation 1.38, it was very simple to adapt it. However, even if the reasons explained in the paper behind the use of an Haar transform into the l 1 -norm term were viable, the actual implementation was counterproductive for the smooth misorientation gradient case treated in this paper. The exact information on the distribution of a single orientation through the whole volume is not known a priori. What is know instead, is that the sum of all the orientations should be homogeneous across the real-space volume if the material is homogeneous itself. This means that the Haar transform should have been applied to the sum of all the orientations instead of the procedure chosen in this paper, where it was applied to each separate volume associated to all the sampled points in orientation-space. Nevertheless, this choice has little influence on the actual result, because the data fidelity term plays the biggest role with the chosen value of λ = 10 -3 , and some sort of smoothing is the only effect associated with the Haar transform on the reconstruction. For the undeformed twin reconstruction instead, the choice of the functional was completely legitimate, and the result is in line with the expectations. We only realized recently that the result of figure 7 is wrong, due to a missing square root in the calculation used to produce it. The correct result is shown in figure 2.1, where the distance between the real solution and the reconstructed solution is plotted as a whole in figure 2.1a, and then broken down into the different axis components in the other sub-figures. In this synthetic data reconstruction, the average orientation of the grain was centered in (0, 0, 0) in Rodriguez space, and in this particular case there is a correlation between the deformation along the z -axis and the ω-spread, while the Lorentz factor only acts as a multiplicative factor on top of that correlation. The biggest deviation that drives the global error of the reconstruction is the one along the z component of the Rodrigues space, and this means that having suppressed the ω information by condensing the diffraction blobs into summed diffraction spots, determines the impossibility to reconstruct the deformation along the rotation axis of the sample. Another interesting observation is that while the error along the x component of the Rodrigues space is quite uniform, the reconstruction along the y component presents one specific point where it becomes significantly bigger than the rest of the volume. This can be understood by inspecting figure 2.2, where the orientation space sampling of the reconstruction and the phantom orientation space distribution are plotted. We can immediately see that both the z -axis and the y-axis were not well covered. Even if especially from figures 2.2b and 2.2c, we see that the z -axis is getting the biggest angular deviations from the sampled region, it is not so obvious from figures 2.2a and 2.2c, that the y-axis is the one that suffers the most in terms of number of voxels that are not well sampled The difference between the actual orientation of the voxels and the coverage of the orientation space by the sampling grid is reflected in the form of higher error in the plots of figure 2.1, and it explains why for the z component we can see two high error regions, corresponding to two corners of the sample volume, which also correspond to the regions that didn't get correctly sampled. However, that explanation only applies to the biggest deviations in figure 2.1d, while the previous considerations represent the best explanation for the generally high reconstruction error throughout the whole volume.

To understand why the chosen sampling of the orientation-space was so bad, compared to the orientation distribution of the phantom, it should be considered that this was a preliminary study, where it was not clear if the correct determination of the bounding box in orientation space was possible, and the technique explained in section 1.2.6 was not developed, yet. The authors decided then to produce a sampling that had the size of the known orientation-spread, and to center it around the average orientation of the synthetic grain, which is the 64 CHAPTER 2. RECONSTRUCTION DISCARDING ω INFORMATION only piece of information that was already available at the indexing stage, at the time.

Introduction

Diffraction Contrast Tomography (DCT) is a non destructive technique used to characterize 3D grain microstructures and access their inner structure in a wide range of polycrystalline materials. Its standard approach is oriented to undeformed materials with negligible intragranular orientation spread. The technique is a truly three-dimensional tomographic imaging approach, sharing a common experimental setup with conventional X-ray microtomography. One of the most advanced and interesting ongoing developments of DCT is in the reconstruction of the local orientation inside of grains, which provides a full 65 characterization of the orientation microstructure of the sample. The aim of this work is to give a framework for determining which orientation is present in every single volume element of a grain.

DCT Experimental Setup

The sample is placed on a rotation stage and irradiated by a parallel monochromatic X-ray beam that is perpendicular to the rotation axis of the sample and whose dimensions are determined by slits.

As the polycrystalline sample rotates, the Bragg condition is met by the different grains at specific angular positions, giving rise to diffraction "spots". The spots correspond to two-dimensional projections of the three-dimensional grain volumes on the detector. The physics behind DCT is the same as for standard crystallography and diffraction measurements, so we will now briefly introduce the basic concepts, before explaining the peculiar properties of DCT. For more information we refer to [9] and [14].

We will now present the basics of the diffraction geometry, but for further details we refer to [12].

Crystal Domains Let us assume first that we describe a three-dimensional crystal lattice with the basis vectors a, b and c. The unit cell is the minimal space spanned by this basis.

The discrete periodic lattice in real space corresponds to a discrete periodic lattice in Fourier space, the so-called reciprocal space. The basis vectors in the reciprocal space will be a * , b * and c * . Given an incoming wave-vector k in and an outgoing wave-vector k out , X-ray diffraction is observed when the vector G = k ink out is close to a vector of the reciprocal lattice. A convenient way of representing G in the reciprocal lattice is by using the Miller indexes, so that G hkl = (h, k, l) T . The link between the direct space Cartesian system and the reciprocal space, indexed by the Miller indexes (h, k, l), is given by the matrix B. The B matrix is an upper triangular matrix that depends on both the direct and reciprocal basis vectors, and the angles between them. The result is a vector G c = BG hkl Diffraction Geometry In this paragraph we will now establish the connection between the diffraction directions, for a specific crystal lattice in a generic reference system, and the points G hkl in reciprocal space.

Let us now restrict ourselves to an individual grain rotated by an arbitrary angle along an arbitrary axis, and then fix one rotation as sample reference. For this grain, there will be a 3 x 3 orthogonal rotation matrix g -1 which aligns the crystal lattice to the sample reference. If we take a specific plane normal, which will correspond to a vector G c = BG hkl , the vector G c will become in the sample reference system:

G s = g -1 G c .
A given laboratory Cartesian reference system will be identified by the real space vectors x l , y l and z l . Assuming the sample reference system to have origin in the origin of the laboratory coordinate system, and to have rotation axis along the vector z, the two systems will be related by the unitary rotation matrix Ω ω which depends on the rotation angle ω. The sample reference system will coincide with the laboratory reference system when ω is 0. So we can finally express the relationship between a scattering vector in the laboratory system h, which is the direction of the diffracted beam, and the Bragg node of indexes (h, k, l) in the reciprocal space, at the given rotation angle ω:

h lab = Ω ω g -1 B   h k l   (1)
Projections of the Crystallographic Domains Assuming that each grain is a perfect lattice, we will observe a projection of the grain volume on the detector for those scattering vectors h lab , which give rise to a different beam intersecting the detector (fig. 1). In the simplest case of orientation unifor-Figure 1: Diffraction of one grain in a DCT experiment mity inside the grain, DCT is able to retrieve the three-dimensional shape of the grains in a polycrystal, by using oblique angle tomography reconstructions, which rely on Algebraic Reconstruction Techniques like the SIRT algorithm [7].

If we allow for regions in the grain to have different orientations from the average, the problem complicates heavily, because the projection geometry is no longer constant across the grain.

The two cases we will take as example in this article are: grains that exhibit a smooth variation of orientation, contained in a small bounded region in orientation space, and grains where sub-domains with a neat change in the orientation are formed, but there exists a notable symmetry operation between the orientation in the grains and the one in their sub-grains, called crystallographic "twins".

Model

The problem of determining the local orientation of a grain has already been heavily studied and many approaches have been proposed, both in terms of models and algorithms. Prominent examples are the forward modelling idea from [15], and GrainSweeper from S. Schmidt discussed in [13]. What they all have in common is that they try to assign orientations to the voxels, but each voxel is treated independently. Another very interesting approach is introduced in [1] and further discussed in [2]. This other approach is somehow closer to our model, but tries to solve the mathematical problems exposed in this section, using Discrete Tomography, with Gibbs priors tailored for Grain Maps.

Here we try to perform a global optimization, where all the voxels are processed simultaneously, using recent algorithms and techniques from signal theory. We also add the constrain of allowing only a discrete and restricted number of orientations. This results in a sampling of the 3D orientation space much similar to the common sampling of the 3D real space implicitly employed by traditional tomography, making this method a 6-dimensional method.

Basic Mathematics

The aim of this article is to show that within certain limits, for DCT measurements, it is possible to identify a discrete set of allowed orientations in the reconstruction.

The actual number of orientations to be taken into account is a compromise between two conflicting requirements: fidelity to the actual orientation distribution on one hand, and indeterminacy of the linear system to be solved in the reconstruction process. The reason lies in the fact that the measurements (diffraction spots) are usually highly convoluted, because in a three-dimensional volume there will be a discrete number of voxels projecting to the same pixel on the detector, and so, if we take too many orientations into account, the amount of information from the spots is not enough for the number of unknowns we allow for.

To explain better this statement we should think that from the reconstruction point of view, the orientation space is a three-dimensional space and the real space is another three-dimensional space. So the full characterization of a grain should be carried out in 6D [11].

If we model the forward-projection operator in tomography as a matrix A, the solution to a reconstruction problem can be seen as the solution to a linear system Ax = b, where we call x * a vector that makes the equality true. We know from linear algebra theory that for a guaranteed unique solution, the matrix A should be square and of full rank (otherwise uniqueness/existence depends on the particular vector b). This means that if we intend to sample the grain volume with a cubic grid of size n x n x n, and the orientation space with a cubic grid of size p x p x p, the number of unknowns will be of order O(n 3 p 3 ). If we now consider that a typical diffraction spot has O(n 2 ) pixels, we would then need O(np 3 ) diffraction spots to fulfill the requirements for a unique solution.

However the number of spots that we can actually successfully record on the detector, and use in the analysis is usually between 10 ∼ 100.

As the measurements will always be perturbed by noise and other experimental inconsistencies, the following formulation is preferred:

x * = min x ||Ax -b|| 2 (2)
where from the feasible solution space we select the vector x * which minimizes the l 2 norm of the residual. Unfortunately the solution to ( 2) is not unique if A does not have full rank. In fact every solution of Ax = b is also solution to (2). Thanks to the advances in the signal theory, it is possible to regularize this problem with the result of selecting solutions with specific properties. The regularization used in these examples is the Lasso formulation over the Haar transform [8] of the volumes:

x * = min x ||Ax -b|| 2 + λ||Hx|| 1 (3)
where H is the Haar transform and λ a weight. This formulation tries to assign a penalty to solutions that exhibit a large l 1 norm in the Haar domain, while the Haar transform is nothing else than a wavelet transform, on the Haar basis set.

The penalty is then weighted by the parameter λ which balances the data fidelity term ||Ax-b|| 2 and the regularization term ||Hx|| 1 . The parameter λ can hardly be related to some physical quantity, and remains for now a purely mathematical term. On the other hand it is also true that each value λ in the problem (3) corresponds to a specific value in the problem:

x * = min x ||Hx|| 1 (4a) subject to: ||Ax -b|| 2 ≤ (4b)
but there is no easy and reliable way of computing one, when the other is known. While ( 4) is more physical and more interesting to the scientific community, algorithms to solve (3) are generally simpler, more efficient and more effective.

In this work we slightly modify (3) to introduce a non negativity constraint for the reconstructed volumes:

x * = min x ||Ax -b|| 2 + λ||Hx|| 1 (5a) subject to: x ≥ 0 (5b)
69 which has both a physical meaning and a practical outcome. On the physical side, the voxels are supposed to indicate a scattering power for the given orientation at the fixed position in the sample, and this scattering power can never be negative. On the other side, imposing the non-negativity constraint, makes the reconstruction algorithm converge quicker and give better results. Finally we would like to briefly justify the choice of the l 1 -minimization over the Haar domain. It was demonstrated in [5] and [6] that l 1 -minimization can result in very accurate recovery of sparse signals, and it is known from the literature that in general it promotes the choice of sparse solutions, for problems where the solution is not unique [4].

So even if all the mathematical requirements for sparse recovery (see e.g. [5]) are difficult to be met for our physical case, where randomness of measurements is not an option, the use of sparsity promotion techniques can still be very useful if coupled with tools like the Haar transform. In fact, the choice of the Haar transform was driven by its well known property of being able to provide with a sparse or at least very compressible representation of transformed objects [8]. Moreover, white noise tends to never have a sparse representation in the Haar domain, because it always overlaps with the regions probed by the wavelets in frequency space. This means that the choice of (3) or (5) will also help in reducing the noise.

Algorithm

An implementation of the FISTA algorithm was used to reconstruct the examples [3]. FISTA is a strongly convergent first order algorithm, oriented to the solution of the Lasso problem. It is composed by three important pieces: a gradient iteration, a soft-thresholding phase and finally a memory step that tries to regularize the convergence of the algorithm.

To understand how these steps can be easily performed in our 6-dimensional problem, it is first important to describe the data-structures used by the algorithm.

We have two main data-structures that hold the information: the diffstack and the solution vector. The diffstack contains all the information recorded by the detector, while the solution vector contains the 6-dimensional volume in orientation and real space. These two data structures are linked by a third object, generically named geometry.

Solution vector

The inner structure of the solution vector is rather simple: it is a linear array of three-dimensional volumes, where each of these volumes has the same size and represents a specific point in orientation space (fig. 2). This means that for the same three-dimensional volume in real space, we allow only specific orientations in orientation space and each of these orientations is represented by a single volume in the solution vector.

Diffstack

The diffstack needs to be covered with more care than the other objects involved. Depending on the intrinsic nature of the grains we are trying If, instead, there is a smooth variation in orientation all along the grain, there will be a three-dimensional volume for each different event, and this volume in ω space is called Diffraction Blob.

While the case of two-dimensional parallel projections is rather simple to be dealt with, because it only requires to be careful enough to relate the correct spots to the correct orientations, the three-dimensional case can be handled in different ways. The so called "Omega spread" of the Diffraction Blob is surely a source of information, but if we don't provide enough sampling points in orientation space (different volumes in the solution vector ), we might have the risk of not covering completely the ω space. This would result in inconsistencies in the linear system. Moreover it is reasonable to think that for a given orientation, which was not sampled in the solution vector, the closest sampled points in orientation space will be approximating it, if the information about it (namely the intensity projected on the detector), can be reached by these nearby points in orientation space.

In case of limited orientation spread inside a grain, it can be beneficial to "squash" the Diffraction Blob into a Diffraction Spot, by summing the blob along the ω direction. The result is a regular diffraction spot, where the more detailed angular information is lost. While this is a really crude procedure, in some simple cases it was observed in simulations that it still provides diffraction Figure 3: On the left we can see a diffstack made of Diffraction Spots, while on the right, a diffstack made of Diffraction Blobs, where around each average ω we have multiple spots, with slightly different ω, that compose Diffraction Blobs spots which can be used to reconstruct relatively good approximations of the expected solution. In the simple case where only Diffraction Spots are used, the diffstack is nothing else than a three-dimensional volume where each spot is one slice of the volume. For more complicated cases, like the case where we want to consider the angular information contained within the blobs, the diffstack becomes a collection of Diffraction Blobs (fig. 3).

Geometry

The geometry is a collection of tables that fully describe the three dimensional scattering geometries for each of the sampling points in orientation space. In practice each line of each table will describe how the related volume will project on a specific spot into the diffstack.

The structure of the geometry table is fixed by the tomographic toolbox used to perform the projection and the back-projection of the volumes on the detector. The table presents 12 columns and as many rows as the number of slices in the diffstack, in fact each row determines the projection geometry for a given Diffraction Spot. The first three columns represent the scattering direction in the laboratory coordinates, while the second triplet of columns represents the position of the center of the detector, always in the laboratory reference system. The other two triplets of columns represent each a unit vector on the pixel grid that represents the detector. Doing so, it is possible to fully determine an arbitrary projection geometry that specifies how each of the volumes in the solution vector project on each of the slices of the diffstack.

Having described the data-structures we will now understand how these objects interact to perform the matrix operations described in the FISTA algorithm. Given the projection matrix A in (3), it is equivalent to the tomographic forward-projection of the volume on the detector in our model. On the other side, the back-projection is equivalent to the transpose A T .

Once we know this, in the gradient iteration of FISTA:

d k = y k - 1 L A T Ay k -b (6)
computing Ax k is nothing else than performing the forward-projection of each volume on the detector space (in a data-structure like the diffstack ), and summing all the contributions of each corresponding pixel, while computing A T of the detector residual will simply mean back-projecting the detector residual in each of the volumes that form the solution vector.

The soft-thresholding step happens in the wavelet domain, so performing the operation:

x k = H T T λt (Hd k ) (7)
is equivalent to computing the wavelet transform of each back-projected volume, performing the soft-thresholding (T λt ), and finally anti-transforming each volume again.

Finally the memory step:

y k+1 = x k + t k -1 t k+1 x k -x k-1 (8)
is simply a sequential operation over all the volumes in the solution vector, where t k is a purely mathematical parameter used to quench the memory step in the first few iterations, and defined in [3].

In our implementation the Haar transform is a simple function in matlab language, while the tomographic forward-projection and back-projection are performed using the ASTRA toolbox, an Open Source library that is able to provide high-performance tomographic primitives, thanks to the use of modern GPUs [10]. Despite possible slowdowns due to the memory transfers between the GPU and CPU memory domains, the generated overhead is almost negligible, compared to the cost of forward-projection and back-projection. Nevertheless, we are now considering to implement asynchronous transfers, in order to hide transfer latencies behind the most time consuming operations.

Results

We will now consider two important physical cases and analyse the results on synthetic data.

Twins

Twins are crystallographic sub-domains where the crystal undergoes a modification that can be represented as a rotation along a well defined direction of the crystal lattice. The result is that some of the scattering vectors will be the same as the ones from the parent domain, and some others will be changed.

If we assume these sub-regions to not be further distorted, the number of sampling points in orientation space to be able to reconstruct all these sub-domains at the same time, becomes exactly the number of these smaller regions inside of the grain.

We are going to assume the correct indexing of the spots was performed and show synthetic data reconstructions. For the realistic case where just one twin is formed, restricting the number of orientations to just two is perfectly physical and vastly reduces the number of needed measurements. Synthetic data was generated using a known grain shape, determined from phase contrast data [9]. The grain volume was then divided in two regions with different size and a different, but symmetry related, orientation was assigned to the regions. One of the two orientations, the one assigned to the smaller region, was chosen to be one of the possible twin variants of the other. The initial grain shape can be seen in figure (4a), where we can actually distinguish the two crystallographic domains by their color: the blue part is the parent, while the thin green region is the twinned part of the crystal. Since parent and twin usually share at least four reflections, given by the plane normal to the axis of rotation, in our example the diffstack generated was formed by 44 projections, 4 of the spots were in common between parent and twin, so that each of the domains was projecting to 24 spots.

As can be seen in figure (4b), the reconstruction almost completely recovers the shape of the two crystallographic domains. There are small conflicts at the interface (pink pixels) but they only represent very small deviations from the original solution.

Figure 5: Geometry of two-dimensional illumination

Smooth gradient deformation

In some materials, grains don't "twin" when subject to mechanical stress, but give rise to smoothly varying misoriented regions of the crystal lattice. This phenomenon affects the quality of the diffraction spots, giving rise to distorted projections. The diffraction geometry is no longer constant through the grain volume, and essentially every voxel will have one or more different orientations, which might slightly differ from the orientations of the close neighbours. This makes it not possible any more to restrict the number of orientations to very few possibilities, but still makes it possible to sample the orientation space around an average orientation, using a discrete number of variations.

For our case study we preferred to build a simplified model that allows us to use diffraction spots, instead of diffraction blobs. The experimental case simulated in our example corresponds to a two-dimensional slice illumination of the grain, that produces deformed spots on the detector. The angular range of misorientation in the example is confined between ±0.5 deg. Having a two-dimensional volume reduces the complexity of the data to the point that the spots are not too convoluted, and can be used by the algorithm without worrying about the ω spread around the average diffraction vector.

As already mentioned, the algorithm remains untouched, with every orientation projecting to every diffraction spot, and the spots built by a voxel-wise sum of the diffraction blobs along the ω direction.

In figure (6) we can see a comparison between the theoretical orientation domains in the sample (6a) and the reconstruction performed using a regular grid in orientation space with 341 sampling points (6b). The sampling grid can be seen as a body-centered lattice, or as two different cubic lattices, of sizes 6 x 6 x 6 and 5 x 5 x 5, where each occupies the centers of the cubes formed by the The image (6a) was produced assigning to each sampled orientation a color from the jet colormap, and then coloring each pixel with the color of the closest sampled orientation in orientation space. The image (6b) was instead generated assigning to each pixel the color corresponding to the orientation with highest intensity in the corresponding voxels, using the same jet colormap as in (6a). In figure (7) we can see the angular deviation in degrees between the theoretical phantom and the reconstruction, which showed a maximum value of ∼ 1.5×10 -3 deg on the top right corner. In this reconstructed volume, the average value between the sampled orientations was weighted over the reconstructed intensities.

Discussion

Even though the results presented in the article seem convincing on theoretical data, experimental data can be affected by deviations from the model that make the reconstructions look noisy or even fail. For twins, an important role is played by the alignment of the geometry, which needs to be carefully fine tuned, and by the re-normalization of the spots' intensity, which is necessary because of the different scattering powers of different (h, k, l) reflections, self-absorption and other effects which are hard to take into account. As a matter of fact, different scattering conditions have different scattered intensities, while tomography assumes that the total intensity of all the spots should be the same. The renormalization that takes place needs a reference which is arbitrary on a single grain basis, but the ratio between the single references is important for the twin case. Generally this reference depends on the total volume of the grains, but this information is not known in advance. In experimental data, we expect the preprocessing logic to be aware of this problem and find a renormalization technique that is able to estimate a good reference ratio between the parent and the twins.

For what concerns the two-dimensional orientation spread case, the approach followed in the article looks already enough for an approximate solution of the local orientation reconstruction problem. However, an approach that exploits diffraction blobs is expected to give better results.

Conclusions

From the reconstructions in section 3 we can confirm that the intuition of sampling the orientation space with a discrete number points can result in a good approximation of the local orientation in grains.

In the twin case, the needed orientations are determined a priori, and their identification is precise enough to make it possible to restrict them to the minimum. This results in a fast a quickly converging algorithm, which also seems to be very reliable on theoretical data.

In the other case, where the orientation spread was taken into account, the reconstruction technique was very simple and not completely adapted to this complex physical case. As already mentioned in 3.2, the results are showing good signs of progress and convergence towards a stable solution of the problem. Chapter 3

acknowledgments

KJB

Theoretical formulation of the six-dimensional model

This chapter will include the publication made on Journal of applied Crystallography (2014) [START_REF] Viganò | Reconstruction of local orientation in grains using a discrete representation of orientation space[END_REF], where the full six-dimensional model and algorithm were developed and discussed in great detail. Compared to the work presented in chapter 2, the work presented here underwent a considerable improvement, both in terms of the six-dimensional model formulation and of algorithmic implementation. Indeed, in this article it is possible to find a clear statement about the difference between the three-dimensional vector field representation, with its important non-linearity implications, and the six-dimensional scalar field representation discussed in section 1.3.3. The diffraction spot approach was completely abandoned, in favour of an approach that would use diffraction blobs instead, because in the meanwhile, the DCT analysis toolchain was able to produce the blob information, instead of just storing the spots into the grain data-structures, making this approach viable. Moreover, the minimisation func- , where the operator L is the identity matrix. This constraint doesn't enforce any really physical constraint, but it takes care of coping with the artificially increased number of unknowns by moving to a fully sampled six-dimensional space. In chapter 4, we will see a new article that surpasses this approach and introduces the Total Variation functional minimisation over the sum of all the orientations, to enforce the homogeneity of each single grain volume. This article represents a significant milestone for the innovation of the DCT technique, because it lays the ground for the work discussed in chapters 4 and 5. However, while the results in this publication were in line with the expectations, some other less important results had to be left out, in order to reduce the size of the manuscript. Hereafter we are going to analyse some additional results about the deformed two-dimensional slice with 1 degree smooth orientation spread case study from the previous chapter that was not included in [START_REF] Viganò | Reconstruction of local orientation in grains using a discrete representation of orientation space[END_REF], because the three-dimensional volume 1 degree smooth orientation spread was considered more relevant and even a bigger challenge to the algorithm. The only difference here with the example in chapter 2 is that a different randomly generated phantom was chosen. This new phantom has a slightly higher inter-voxel misorientation, and so it presents a more challenging distribution in orientation space. Also the coloring of the orientation of each voxel has been changed, in order to make the errors more obvious to spot to the human eye. In fact, in figure 3.1 it is possible to see a comparison of the results from the two reconstruction methods on a five-dimensional problem (slice illumination): in figure 3.1b the diffraction blobs are squashed into diffraction spots, omitting the ω information like for the case of chapter 2, while in figure 3.1c the result using the full information from the blobs. In figure 3.2 instead, the error of the reconstruction in figure 3.1 is presented. Indeed it can be seen from figure 3.3, in contrast to figure (5, d) in the attached paper, that for a given voxel in the reconstruction with the diffraction spots, the algorithm is not able to reconstruct the deformation along the z -axis, which is the axis of rotation of the sample. This last observation is coherent with the conclusions reached in the previous chapter. In figure 3.4, similarly to figure 2.2 in the previous chapter, it is possible to see the distribution in orientation space of the voxels in the phantom, but in this case overlaid with the distribution of the voxels in the reconstruction. Some deviations are indeed visible, but the general trend is in line with the expectations. Finally, in figure 3.5 it is possible to compare the forward projection to the detector of the blob reconstruction, with the actual data that was used for the reconstruction. The difference between figure 3.5a and figure 3.5b is that in the first case only a slice of one blob is used for the comparison, while in the second case one full blob has been squashed into a spot. We formulate the problem of orientation reconstruction in terms of an optimization over a 6D space X 6 = R 3 ⊗ O 3 , constructed from the outer product of real and orientation space, and we provide a strongly convergent first order algorithm that makes use of modern l1-minimization techniques, to cope with the increasing number of unknowns introduced by the six-dimensional formulation of the reconstruction problem. We then have a look at the performance of the new reconstruction algorithm on synthetic data, for varying degrees of deformation, both in a restricted line beam illumination, and in the more challenging full beam illumination. Finally we show algorithm's behavior when dealing with different kinds of noise. The proposed framework, along the reconstruction algorithm, looks promising for application to real experimental data from materials exhibiting intra-granular orientation spread of up to a few degrees.

Introduction

Over the past decade considerable efforts have been put into the development of three-dimensional X-ray diffraction techniques for structural characterization of polycrystalline materials. The ultimate goal of these grain mapping techniques is the non-destructive description of a materials three-dimensional micro-structure in terms of local phase and crystal orientation, the so-called micro-texture. For a review of the state of the art in this field, the reader is referred to the special issue on three-dimensional diffraction microscopy techniques of the current journal [4] and the book by Barabash and Ice [3]. In the general case, the crystalline micro-structure of a volume element of a polycrystalline material may have to be described in terms of a three-dimensional orientation distribution function, allowing for multiple orientations to be present in each volume element. However, depending on the size of the volume element and the deformation state of the material, simplified but still adequate representations of the micro-texture may be obtained by assigning an average orientation to each volume element and, for undeformed materials, even a single (average) orientation per grain may be sufficient. The general, six-dimensional framework for micro-texture analysis has been discussed by [15], who suggested that the use of algebraic reconstruction techniques may prove a viable route for microtexture analysis and related, lower dimensional sub-problems.

In this article we focus on near-field variants of the monochromatic beam rotation method like three dimensional X-ray diffraction microscopy (3DXRD) [17] and X-ray Diffraction Contrast Tomography (DCT) [19], well adapted for mapping 2D and 3D grain micro-structures in materials where the aforementioned simplifying micro-structure descriptions are applicable. Over the past years, a variety of solutions for sub-cases of the general problem of micro-texture analysis have been presented and remarkable progress has been made using algebraic reconstruction techniques as well as reconstruction strategies based on forward modeling and/or combinatorial optimization. For an overview of this work the reader is referred to [12], [17] and references therein. Restricting the illumination of the sample to a single slice through the volume, these methods have proved capable to produce orientation maps from metallic samples having undergone ten percent and more plastic deformation. From an experimental point of view, the restriction of the sample illumination to individual slices compromises the temporal resolution and may result in anisotropic voxel size in three-dimensional reconstructions obtained from stacking these layers. For this reason the development of algorithms allowing micro-structure reconstruction from three-dimensional diffraction data is highly relevant. DCT is an example of a truly three-dimensional tomographic imaging approach, sharing a common experimental setup with conventional X-ray micro-tomography. The algebraic reconstruction approach behind DCT may be considered as one of the sub-cases of the general six-dimensional framework, tailored to undeformed materials exhibiting limited (≤ 0.5 • ) intra-granular orientation spreads. In this case the orientation degrees of freedom inside each volume element are neglected and three-dimensional grain shapes are recon-structed assuming a single, constant orientation throughout the grain volume. In this work we introduce a six-dimensional extension of the three-dimensional tomographic reconstruction approach behind DCT, extending the applicability of this method to materials exhibiting intra-granular orientation spreads of up to several degrees. While a previous preliminary work in this direction already exists [24], this is a major extension, being the first extensive and complete treatment of the model and the reconstruction algorithm. In order to account for spatially varying orientations inside a grain and potentially also the presence of multiple orientations inside each of the individual volume elements, we describe a discrete, six-dimensional representation of the reconstruction problem in the form of a direct product of real space and orientation space [15]. The model is based on the assumption of kinematic diffraction, implying proportionality between crystal volume and integrated, diffracted intensity. We assume the grain average orientation has been determined by one of the existing polycrystal indexing approaches [11], [21], [19], [20], and that the experiment has been performed in such way that diffraction signal of different grains can be separated on the detector (negligible overlap with diffraction spots from other grains). In section 1.1 we summarize basics concepts of 3D diffraction imaging approaches like DCT and 3DXRD. We then describe in 1.2 how different properties of the analysed materials give rise to different types of images on the detector. In section 2 we introduce the proposed mathematical model for representing the problem. Section 3 is then dedicated to the application of the reconstruction algorithm to synthetic data generated from a grain with one degree of intra-granular orientation spread. Finally in section 4, we test the strength of the model and the algorithm, when applying them to more challenging cases, with higher degrees of orientation spread, and discontinuities (small angle boundaries) in orientation between neighbouring blocks. In the appendices we decided to add some extra considerations and details, like the discussion of some implementation details of the algorithm (Appendix A), and a robustness test of the algorithm when subject to extinction noise (Appendix B).

Near-field diffraction imaging setups

The sample is placed on a rotation stage and irradiated by a parallel monochromatic X-ray beam that is perpendicular to the rotation axis of the sample and whose dimensions are determined by slits. The data recorded during a nearfield diffraction imaging experiment will be over a range of 2π in the form of s = 2π/δω images, where δω is the angular range over which the signal is integrated on the detector, to form a single image. Typical values for δω are in the range between 0.05 and 0.2 degree. As the polycrystalline sample rotates, the Bragg condition is met by the differ- ent grains at specific angular positions, giving rise to diffraction "spots". For undeformed grains these spots correspond to two-dimensional projections of the three-dimensional grain volumes on the detector. As it can be seen in figure (1), a diffraction imaging experiment can be both performed with a full 2D beam or a restricted 1D line beam. The advantage of the second is that the beam dimension is very small in the direction parallel to the rotation axis which in turn reduces the complexity (convolution) of the reconstruction task (from 6D to 5D problem). The physics behind 3DXRD and DCT measurements has been outlined in previous work, for completeness we recall some basic equations for calculation of the diffraction geometry following the presentation in [16] and [23]. For more information on DCT we refer to [13] and [19]. We will consider the case of a face-centered cubic crystal and experimental parameters typically used in a near-field diffraction imaging experiment. The detector fully intercepts the first three diffraction cones, giving rise to about 60 diffraction spots per crystal. We recall the relationship between a scattering vector in the laboratory system G lab , which is the direction of the diffracting plane normal, and the Bragg node of indexes (h, k, l) in the reciprocal space, at the given rotation angle ω:

G lab = Ω ω g -1 B   h k l   (1)
where Ω ω is the rotation matrix associated with the angle ω, g is the orientation matrix, and B is the metric matrix. Bragg diffraction is observed for rotation angles ω fulfilling: where the incident wave vector is k i = kx [23]. Equations ( 1) and ( 2) fully determine the projection geometry for a volume element with a given orientation g, and can be used to derive the positions (u, v) and rotation angle (ω) onto which a given voxel will project on the detector.

G labx = - |G hkl | 2 2k ( 2

Projection of the Crystallographic Domains

Having introduced the diffraction geometry, we now describe the appearance of projection data for two different cases.

1. Undeformed grains. Assuming that each grain is a perfect lattice, we will observe a projection of the grain volume on the detector for those scattering vectors G lab , which give rise to a diffracted beam (in the direction k out = G lab + k in ), intersecting the detector (fig. 1 and fig. 2,a). In this case the shape of undeformed grains can be reconstructed by using oblique angle tomography reconstructions, which rely on Algebraic Reconstruction Techniques like the SIRT algorithm [10], [19]. Having the same orientation, all voxels will give rise to diffraction signal in the same direction, for the same angles of rotation ω n . In this idealized situation the diffracted beam will be parallel and each of the images on the detector will be a geometrical projection of the grain volume These projections will be called Diffraction Spots.

2. Deformed grains. A physically more relevant case corresponds to the situation when the whole grain volume diverges from the average grain orientation, but exhibits either smooth variations between nearby regions, or a combination of smooth variations and discrete jumps. This case will be investigated further in this article and we will provide various reconstruction examples with different values for the maximum degree of orientation spread.

As mentioned earlier, we model our grain as if it were made of small regions, each having a specific, average orientation, and so identifiable as a point in a three-dimensional representation of the orientation-space ("vector field"). Integrated over the grain volume, these points will form a dispersion in orientation-space and we define the maximum orientation spread of this distribution, as the diameter of the smallest ball in orientation-space that contains all of its points. We will restrict ourselves to the case of moderate orientation spread of up to a few degrees. As discussed later, this limit arises not only from memory restrictions by currently available hardware, but also by the increasing overlap of the diffraction blobs on the detector. Moreover, for higher levels of deformation, the vector-field description (one average orientation per voxel), may no longer be appropriate.

Let us now discuss more in detail how the data presents itself when dealing with grains affected by smoothly varying orientation. As it can be seen in figure (2, b), as opposed to figure (2, a), a grain doesn't produce a single detector image per reflection any more. Different regions of the grain will project on a range of neighboring images, forming a three-dimensional object per reflection in the uvω space, where u and v are the two-dimensional Cartesian coordinates of a detector image and ω is the rotation angle. Moreover, the directions of the diffracted beams will no longer be parallel, giving rise to distorted projection images of the grain volume. These distorted projection volumes are called Diffraction Blobs. Note that 3D formulations of the reconstruction problem like in [9] and [13] neglect this type of distortion, and use the spots resulting from the integration of the blobs along the ω-direction, treating them as parallel-beam projections of the grain volume, and applying a regular oblique-angles SIRT algorithm on them.

Model

The problem of determining the local orientation of a grain has already been heavily studied and many approaches have been proposed, both in terms of models and algorithms. Prominent examples are the forward modelling idea from [23], and GrainSweeper from S. Schmidt discussed in [17]. What they all have in common is that they try to assign orientations to the voxels, but each voxel is treated independently. Typically operating on less convoluted 1D line beam data, these algorithms have proven stable, even when neglecting local diffracted intensities for performance reasons.

Another interesting approach is introduced in [1] and further discussed in [2]. This other approach is somehow closer to our model because it operates a grain by grain optimization, but it tries to solve the mathematical problem in a 3D "vector field" representation, using discrete tomography, with Gibbs priors tailored for grain maps.

Here we try to perform a global optimization for the full 6D problem (2D full beam illumination), where all the grain voxels are processed simultaneously and matching of diffracted intensities is part of the optimization, using recent Figure 3: A three-dimensional volume in real-space, can be considered a single point in a three-dimensional orientation-space. By taking the corresponding real-space positions in each of these volumes, we would be able to construct a single-voxel ODF for each of the real-space positions.

algorithms and techniques from signal theory. We also add the constraint of allowing only a discrete and restricted number of orientations. This results in a sampling of the 3D orientation space similar to the common sampling of the 3D real space implicitly employed by traditional tomography, making this method a 6-dimensional method.

Discrete representation of the 6D reconstruction space

In order to add orientation degrees of freedom to the real space voxels in our grain volume, we introduce a discretized, six-dimensional reconstruction space X 6 = R 3 ⊗O 3 as the outer product of real space and orientation space. However, only a discrete set of orientations, sampled around the grain average orientation and confined to a small ball in orientation space are included in the analysis. We decided to represent the O 3 orientation-space as a Rodrigues space. Since we consider only small deviations from the known grain average orientation this orientation sub-space can be treated as Euclidean and both 3D spaces are tessellated using cube shaped voxels. From a material science point of view, this 6D space could be viewed as a collection of 3D voxellated orientation spaces, one for each of the volume elements, or in other words, with a discrete ODF per real-space volume element. For performance reasons discussed in appendix A, we organize the two subspaces in the inverse order, which means that we represent the 6D space as a collection of real-space volumes, one for each of the sampled orientations (fig. 3). We will use the vector x ∈ X 6 as a synthetic representation for any arbitrary volume in our 6-dimensional reconstruction space. If we assume that the grain volume was divided in a grid of size n×n×n, and the orientation space in a grid of size p × p × p, the length of the vector x is n 3 p 3 , and the scalar entries can be interpreted as "scattering powers" of the corresponding element. For a grain that can be described by a "3D vector field" (e.g. one 3D Rodrigues vector per real-space voxel), the majority of elements of this vector in X 6 will not contain any intensity.

If we model the forward-projection operator in tomography as a matrix A, the solution to a reconstruction problem can be seen as the solution to a linear system Ax = b, where we call x * a vector that makes the equality true, and b represents the images (pixel intensities) recorded by the detector. Note that the choice of having a fixed set of allowed orientations, decided a priori, makes the projection matrix a static object and the system of equations is linear.

To summarize, in the model proposed by this paper, the x vector is constructed as a collection of 3-dimensional real-space volumes, that correspond each to a specific point in orientation-space. This means that the projection matrix A will be a stack of projection matrices, one for each selected orientation, and that the linear system can be rewritten as the sum of each sampled orientation:

Ax = A 1 x 1 + A 2 x 2 + . . . + A N x N = N i A i x i = b (3)
where N is the number of sampled orientations.

Mathematical model

We will now first show why the reconstruction problem, exposed in section 2.1, is under-determined, and then introduce our model, that tries to cope with this situation.

We know from linear algebra theory that for a guaranteed unique solution, the matrix A should be square and of full rank (otherwise uniqueness/existence depends on the particular vector b). This means that if we intend to sample the grain volume with a cubic grid of size n x n x n, and the orientation space with a cubic grid of size p x p x p, the number of unknowns will be of order O(n 3 p 3 ). A typical diffraction blob has dimensions of order O(n 2 t) pixels, where t is the so called ω-spread, corresponding to the number of images onto which the blob extends. The requirement in terms of number of blobs, for a unique solution, is of order O( n t p 3 ). Note that grains with a smaller orientation spread will have smaller values of t. This also implies that we will need fewer sampling points in orientation space, and so a reduction of the number p. In fact we can conjecture that there will be a linear relationship between t and p, thus reducing the requirement of blobs to be of order O(np 2 ). Assuming 1deg maximum orientation spread, and 0.1deg angular integration steps (δω), the blobs will be ∼ 10 slices thick, and so t = 10. If then we also assume a 50 × 50 spot size and an orientation sampling that exhibit 0.1deg distance between the sampled points, we will have that n = 50 and p = 10. Thus the required number of blobs will be of order 50 × 10 × 10 = 5000. However the number of blobs that we can typically record on the detector in a single scan, and use in the analysis is usually between 20 ∼ 100, making the reconstruction problem heavily under-determined. As stated earlier, in our simulations we chose a fixed number of 60 blobs for all the test cases. As the measurements will always be perturbed by noise and other experimental inconsistencies, the following formulation is preferred:

x * = argmin x ||Ax -b|| 2 (4)
where from the feasible solution space we select the vector x * which minimizes the l 2 norm of the residual. Unfortunately the solution to ( 4) is not unique if A does not have full rank. In fact (4) implies that x * makes the equality in ( 5) true, and we see from ( 6) that ( 4) is equivalent to Ax = b.

∇ (||Ax * -b|| 2 ) = 0 (5)
So every solution of Ax = b is also solution to ( 4), but the minimization formulation better deals with noise.

∇ (||Ax -b|| 2 ) = ∇ (Ax -b) T (Ax -b) (6a) = 2A T Ax -2A T b (6b) = 2A T (Ax -b) (6c)
To regularize the problem (4), we can now impose some constrains that help in selecting specific solutions which exhibit the physical properties of real world samples. Under the assumption of a fine enough sampling in real-space we can expect that only few orientations in the six-dimensional space will be active for each real-space voxel, so that the tools for sparsity recovery will play an important role in selecting solutions which reflect the properties of real world samples.

Normally we would like to add a penalty term over the l 0 semi-norm, because this semi-norm would count the non-zero entries in the solution vector, and a minimization over it, would then yield a sparse solution. On the other side, using the l 0 semi-norm would require combinatorial search of the solution, through all the solution space. We can instead think of introducing a penalty term that uses the l 1 norm and for which simpler and better performing algorithms exist [7]. It was demonstrated in [6] and [8] that l 1 -minimization can result in very accurate recovery of sparse signals, and it is known from the literature that in general it promotes the choice of sparse solutions, for problems where the solution is not unique [5].

So even if all the mathematical requirements for sparse recovery (see e.g. [6]) are difficult to be met for our physical case, where randomness of measurements is not an option, the use of sparsity promotion techniques is a key element of the microtexture reconstruction approach proposed in the current work. We introduce this additional regularization term by exchanging the formulation (4) with a more powerful Lasso formulation, where the minimized l 1 norm would be in the full six-dimensional space:

x * = argmin x ||Ax -b|| 2 + λ||x|| 1 (7a) subject to: x ≥ 0 (7b)
where λ is a weight that balances the data fidelity term ||Ax -b|| 2 , from (4), and the regularization term ||x|| 1 . This formulation assigns a penalty to solutions that exhibit a large l 1 norm in the reconstruction domain.

The non-negativity constraint for the reconstructed volumes has both a physical meaning and a practical outcome. On the physical side, the voxels are supposed to indicate a scattering power for the given orientation at the fixed position in the sample, and this scattering power can not be negative under the assumption of kinematical scattering. On the other side, imposing the non-negativity constraint, makes the reconstruction algorithm converge quicker and give better results.

Algorithm and Data-structures

Having described the model that we wish to solve, we will now present an algorithm instance that can solve the same Lasso formulation. Recently, a new class of first-order primal-dual algorithms was proposed. This class of algorithms is called Chambolle-Pock [7], [22].

The CP algorithms can solve different types of optimization problems, and it is relatively easy to mathematically derive the algorithm tailored for the problem given in (7).

Algorithm 2.1: Synthetic description of CP Algorithm Input: x 0 , y, λ Output: Reconstruction in x k

1 p 0 ← 0 2 q 0 ← 0 3 θ ← 1 4 for k ← 1 to n do 5 p k ← p k-1 + σ 1 (Ax k-1 -y) 1 + σ 1 6 q k ← λ q k-1 + σ 2 x k-1 max(λ1, |q k-1 + σ 2 x k-1 |) 7 x k ← P 0 x k-1 -τ A T p k -τ q k 8 x k ← x k + θ (x k -x k-1 ) 9 end

94

The algorithm (2.1) is made of few important pieces: a projection in both the projection data space and the space associated to the regularization condition in the objective function, the back-projection, and a memory step.

To understand how all these steps in the algorithm can be performed in our 6-dimensional problem (eq. 7), we will now introduce the underlying datastructures. We have two main data-structures that hold the information: the diffstack (vector b) and the solution vector (vector x). The diffstack contains all the information recorded by the detector, while the solution vector contains the 6-dimensional volume in orientation and real space, introduced in section 2.1. These two data structures are linked by a third object, generically named geometry (matrix A). Each line of the matrix A contains the contribution of each element of x to one single pixel on the diffracted images, but using the ASTRA toolbox [14] it is possible to use a simpler description, which consists of a collection of tables that fully describe the projection of the 3D volumes onto the 2D detector images. In appendix A we will explain these concepts in more detail, and we will discuss the implementation details concerning the algorithm.

Reconstruction of slightly deformed materials

We will now show the result of reconstructions using the multiple orientations and diffraction blobs, demonstrating the improvement over the standard DCT approach in reconstructing a 3D cubic grain that exhibits smoothly varying deformation gradients in orientation-space, with 1deg of total orientation spread. Synthetic diffraction data was generated from a cube shaped phantom grain with 50 × 50 × 50 voxels using geometrical and experimental settings identical to those employed in a previous measurement. The diffraction images were simulated using a previously existing code [18], based on a (continuous) 3D vector field description of the orientation inside the grain. Before going into a full analysis of the orientation reconstruction, we can already show some results that concern the grain shape reconstruction (i.e. the signal integrated over the orientation space). As it can be seen in figure (4), the shape of the reconstructed grain is greatly improved by the use of the 6D algorithm. The gray-scale intensity in the said figure corresponds to reconstructed material density (or scattering power) for each real-space voxel. The features seen in the reconstruction of figure (4,c) are artifacts, probably caused by missing orientations, which are not well represented by the orientation sampling, and which result in a lower reconstructed intensity. The orientation sampling is a body centered lattice of 341, constructed from a

Visual Representation of Reconstructions

In order to ease analysis of the reconstruction results, we introduce three types of visualization schemes. Two of these visualization methods produce a colouring of all the voxels in a slice or set of slices, either based on orientation domains, or on the voxel-wise distance from the correct orientation. The third is a detailed visualization of the orientation-space distribution of intensities for a single voxel.

Domain Colouring

This method uses a sampling of orientation space different from the reconstruction sampling, and it assigns a colour to each of the points in the new sampled grid, called the colouring lattice. For the phantom, each voxel is then being coloured using the colour corresponding to its closest orientation in the colouring lattice. The same applies to the reconstructed volume, but the intensity weighted average orientation of the voxels, above a certain threshold, is taken as the reference orientation for determining the distance from the colouring lattice.

Distance Colouring This method takes the average reconstructed orientation of each single voxel, and computes the distance in degrees from this average point and the true orientation of the voxel in orientation-space.

Voxel Orientation Distribution This visualization technique concentrates on one single real-space voxel and it performs the three-dimensional plotting of the intensities assigned to each sampled orientation. In this representation, only the orientations that have intensity above a certain background are plotted, along with some few other special points. The true orientation of the voxel always takes the magenta color, the average orientation Finally the other points, colored using matlab's jet colormap for encoding the intensity, are the active orientations for the selected real-space voxel. The size of the points representing the active orientations is always slightly smaller than the size of the other special points.

Orientation reconstruction

For the reconstruction in figure (4), as it can be seen in the comparison between fig.

(5, a) and (5, b), the reconstruction algorithm is able to retrieve a good approximation of the true solution in all the three dimensions of the real-space.

Another interesting analysis of this test case can be found in appendix B, where we look at the reconstruction performance when the sample is affected by some 97 We conclude this section by looking at the error plot of this reconstruction (figure 5,c), which shows a relatively low angular voxel-wise error for the reconstruction. Moreover, the error is maximum on the borders of the grain. This could be expected in a non-perfect recovery, because that is the region at the border of the grain volume diffracted intensities transition to zero. As a consequence, even small errors in this region can become comparable in modulus to the total intensity of the voxel, and so resulting in higher errors in the reconstruction.

On the other hand, looking at figure (5, d), which shows which orientations are active in one real-space voxel, we see that the l 1 -minimization is effectively working.

Higher deformations

In many practical cases, the intra-granular orientation spread in a real material may be bigger than in the example of the previous section, and reach values of several degrees. In this section we will analyse how the algorithm behaves at higher levels of intra-granular orientation spread (5 degrees, including discontinuous jumps), and which factors influence the result.

In cases where we encountered memory limitations, we resort to using a less dense sampling of orientation space and the pseudo-slices approach, as described in Appendix A.

Acting on the number of orientations

We consider the case of 1D slice illumination of a phantom grain with 5 degree of orientation spread. An interesting result from the progression in figure (6), is that increasing the number of orientations, and so the resolution of sampling in orientation-space, does actually give a measurable improvement for the reconstruction quality. So, the algorithm is actually able to cope with the increasing number of unknowns, and the sparsifying techniques do seem to really help in What is clear from this picture is that, apart from some hot spots, the global error is normally around (if not lower) to the 10% of the maximum orientation spread.

Full beam illumination

Unfortunately, moving to finer sampling in orientation space is problematic when working with full-beam illumination and 3D real-space volumes. The size of the volumes and the blobs can be quite big and the increase in orientationspace resolution could result in the computer going out of memory. Nonetheless, as it can be seen in figure (8), the 6D reconstruction using 559 orientations is in-line with the equivalent reconstruction for the 2D real-space. While the number of orientations was quite limited for the given problem, due to memory constrains, the fact that the reconstructions presents the features of the phantom is an interesting result.

Sub-grain Boundaries

Real materials may display discontinuous jumps in the orientation, so-called subgrain boundaries of only few degrees or even fractions of degrees. This kind of grain sub-structure is usually observed along with the smooth orientation gradients considered previously and we decided to analyse the performance of the algorithm on a test case, having four sub-grains and smooth orientation in each sub-domain. In fact, even if the sub-grain discontinuities could lead to separable diffraction spots, we would observe overlap on the detector due to the smooth deformation.

For solving this problem, a higher number of orientations is needed and so we decided to only analyse the case of line beam illumination (5D reconstruction). Figure (9, a) is a plot of the inter-voxel misorientation, where both the smooth orientation changes and the sharp sub-grain boundaries can be seen. As it can be seen in figure (9), the introduction of the jumps doesn't completely break the reconstructions. This is another very interesting result, because even if the boundaries are regions of higher error especially when reaching the borders (fig. 9,d), the bulky part of the sub-domains seems to be reconstructed in a similar way to the previous 5deg orientation spread test case.

Discussion

The model and algorithm proposed in this work perform a global optimization over the six-dimensional inverse problem defined as the local orientation reconstruction of the crystal lattice in grains, using data from Diffraction Contrast Tomography or other variants of near-field diffraction imaging experiments using X-rays or neutrons. Having an algorithm that can perform such kind of optimization for 2D (extended beam) illumination, by taking local diffracted intensity explicitly into account, is quite critical for time-resolved experimental studies and a prerequisite for studies carried out on low brilliance sources like laboratory or neutron sources.

The model is, however, still based on a number of strong assumptions. The first is the assumption of kinematical diffraction, which in most practical cases is not fulfilled. Primary and secondary extinction as well as absorption may give rise to deviations from the idealised case of mathematical projections of the diffracting grain volumes. The impact of these perturbations can to some extend be reduced optimizing the experiment conditions, for instance via the choice of the X-ray beam energy. Moreover, as shown in appendix B, the algorithm can handle moderate levels of "extinction" noise.

Next, the model explicitly ignores elastic distortion of the crystal unit cell. Introducing elastic strain, which can be modeled with a 2nd rank tensor, in the problem, would add additional six dimensions. Unfolding and sampling all the twelve dimensions would be not feasible in terms of computational power and memory needs. Fortunately, many crystalline materials have elastic limits below one percent and the errors introduced by elastic distortion of the crystal lattice may be considered as small perturbations, compared to the misorientations of up to several degrees discussed in this article.

Another assumption, tightly liked to the choice of an l 1 minimization over the six-dimensional space, is the postulate of a limited number of orientations per voxel. Depending on the actual size of a voxel and the type and deformation state of the material, voxel-wise sparsity in orientation space may or may not be a good description of the materials' micro-structure. In other words, the physical relevance of the selected solution having minimum l 1 norm may depend on the experiment conditions. It would be interesting to evaluate the performance of the current framework in cases where each voxel actually contains an orientation distribution instead of a single orientation. Finally we mention that the current approach treats each grain separately, which in turn translates in the need for limited diffraction blobs overlap with the blobs coming from other grains. The number of grains which can be simultaneously illuminated during data collection is inversely proportional to the intra-granular orientation spread of the grains. During an experiment, this condition can be adjusted to some extent by appropriate selection of the sample dimensions and setting of the beam-defining slits. We expect that the current approach will overcome some of the common problems encountered in tomographic grain reconstruction from near field diffraction data. First of all, reconstructions based on this new approach are expected to outperform tomographic reconstructions treating integrated diffraction spots as parallel projections of grain volumes in terms of accuracy of reconstructed grain shapes. As shown in the previous sections, operating in a six dimensional reconstruction space results in much improved accuracy of shape reconstructions for grains containing some non-negligible orientation spread. This article was focused on testing and validating the model and the algorithm against a number of synthetic test cases. As a next step, the performance of the algorithm will be evaluated on experimental data, using EBSD and other 3D orientation mapping techniques for cross-validation. By adapting the projection and back-projection operators to the case of polychromatic and/or cone beam geometry, the same algorithmic framework could potentially be used to solve the inverse problem of orientation reconstruction from laboratory X-ray or neutron diffraction imaging experiments. For polychromatic variants of diffraction imaging orientation gradients inside a grain give rise to noticeable distortions of the diffraction spots both in azimuthal and radial direction, thereby increasing the need to take these kind of distortions into account.

Conclusions

We conclude this study of a 6D orientation imaging framework, stating that it looks promising for future developments and applications to real experimental data, where the intra-granular orientation spread is within a few degrees. While the current implementation, when scaling to higher levels of orientation spread, hits the computational and memory limits in the current desktop hardware, we showed that for the study of limited inter-granular orientation spread of up to a few degrees, the algorithm can be tuned to give encouraging results, by acting on the orientation resolution parameters.

The reformulation of the non-linear inverse problem of six-dimensional microtexture analysis into a linearized version of the same, and the use of globally convergent algorithms for the given model, give a solid theoretical ground for the future developments in the same direction. Diffstack It is represented in figure (10), and it is the data-structure that holds the diffraction data recorded by the detector. As discussed in section 1.2, data can be of different types: a collection of two-dimensional images or threedimensional volumes (fig. 10). Note that different orientations will typically project to different slices in the blobs, as depicted in figure (11).

Geometry Mathematically speaking, the geometry describes how the matrix A acts on the solution vector, to produce the diffstack.

As a physical consideration, the matrix A doesn't include a full scattering model, and so the Lorentz factors and atomic form factors of the blob will be ignored. This makes it impossible to use the absolute intensities of the blobs, but an easy way to avoid this problem is to re-normalise all the blobs to the same integrated intensity, because we assume that the relative scattering power between the voxels does not change as a function of the reflections. The structure of the geometry table is fixed by the tomographic toolbox used to perform the projection and the back-projection of the volumes on the detector. The table presents 12 columns and as many rows as the number of slices in the diffstack, in fact each row determines the projection geometry for a given Diffraction Spot. The first three columns represent the scattering direction in the laboratory coordinates, while the second triplet of columns represents the position of the center of the detector, always in the laboratory reference system. The other two triplets of columns represent each a unit vector on the pixel grid that represents the detector. Doing so, it is possible to fully determine an arbitrary projection geometry that specifies how each of the volumes in the solution vector project on each of the slices of the diffstack. The actual expansion of the geometry tables into the projection matrix A used in the reconstruction process is handled internally in the toolbox.

Figure 11: On the left we can see a diffstack made of Diffraction Blobs, and on the right the full volume associated to each sampled orientation. The lines connect each full volume associated to a sampled orientation, and the slice in the blobs that the orientation reaches.

Having described the data-structures we will now provide a description of how these objects interact to perform the matrix operations described in the CP algorithm. Given the projection matrix A in (7), it is equivalent to the tomographic forward-projection of the volume on the detector in our model. On the other side, the back-projection is equivalent to the transpose A T . In a simple case like the algorithm that only deals with the integrated diffraction spots, computing Ax k is nothing else than performing the forward-projection of each volume on the detector space (in a data-structure like the diffstack ), and summing all the contributions of each corresponding pixel, while computing A T of the detector residual will simply mean back-projecting the detector residual in each of the volumes that form the solution vector. For more complex problems like the one dealing with the diffraction blobs, one intermediate step needs to be added for practical reasons. It is indeed easier to let each orientation project the associated volume to an intermediate projection data-structure, containing N discrete ω slices with N the number of observed reflections. The ω slices of the blob-diffstack receive the summed contribution from each volume that projects to that slice (fig. 12). This approach, shown in figure (12), has the advantage that since the projection and back-projection are computed on GPU devices, it doesn't require the full blob-diffstack to be loaded on the GPU memory at the same time. An important parameter in the model is the distance of sampling points in orientation space which in turn governs the distance between real and reconstructed orientations in the solution vector, as well as the distance between the discrete ω slices, reached by these sampling points. Having only a fixed number of orientations, for large deformations, they might not be enough to project to all the ωs in each blob, giving rise to inconsistencies in the linear system. Figure 13: Grouping of slices: (a) most simple approach to let the orientations reach all the slices, (b) pseudo-slices approach, where we perform a linear interpolation of the original slices in the blobs.

A simple way to avoid this problem is to subdivide the slices in a blob into groups, from which pseudo-slices (or sub-spots) are produced. This technique reduces the angular resolution of the scans, but helps in relieving the problems coming from undersampling. There can be different strategies on how the pseudo-slices can be produced. The simplest strategy consists in summing all the images in a group, to produce a diffraction blob with an reduced angular resolution (fig. 13,a). While this approach can be considered quite crude, it has indeed some sense, and for groups made of a small number of slices, probably little difference can be observed from other methods. On the other hand, this method doesn't help in letting close orientations talk to each others. What was in fact observed in the previous article [24], is that if a voxel had an average orientation which was not sampled but between other sampled orientations, the voxel would gain signal from those close orientations, giving an average orientation, close to the original orientation.

To try to better allow this behavior, we developed another approach, showed in figure (13, b), which builds the pseudo-slices from the linear interpolation of the original slices in the blobs. This is done by selecting a level of interpolation n, picking the geometry of one blob slice every n, and building the corresponding pseudo-slice from the 2(n -1) blob slices around it. The assigned coefficients, 105 In terms of implementation, the memory step of the algorithm

x k = x k + θ (x k -x k-1 ) (8)
is simply a sequential operation over all the volumes in the solution vector, where θ is a purely mathematical parameter, introduced in [7].

In our implementation the tomographic forward-projection and back-projection are performed using the ASTRA toolbox, an Open-Source library that is able to provide high-performance tomographic primitives, thanks to the use of modern GPUs [14]. Despite possible slowdowns due to the memory transfers between the GPU and CPU memory domains, the generated overhead is almost negligible, compared to the cost of forward-projection and back-projection. Nevertheless, we are now working on the implementation of asynchronous transfers and asynchronous computations on both CPU and GPU, in order to hide transfer latencies and blobs↔sinograms transformations behind the most time consuming operations.

B Stability to noise

In this section we will analyse the robustness of the algorithm to the presence of a common type of noise in a diffraction experiment, when crystals show little to none deformation. So we performed this test for the 1deg deformation case only.

Other sources of inconsistencies in the projection data, like shot noise in the pixel signals were not considered. Typically in DCT experiments, only the brightest diffraction spots/blobs are segmented and then selected for the reconstruction of the single grain volumes. Given that the shot noise for a pixel grows with the square root of its intensity, and the fact that DCT uses only the brightest spots, this type of noise can be considered negligible.

Figure 15: Illustration of the origin and apprearence of noise related to secondary extinction in a polycrystalline sample.

B.1 Extinction and attenuation effects

The source of noise affecting the measurements of diffracted intensities on the detector, that we are now going to consider, is related to extinction and attenuation of the X-ray beam travelling through the polycrystalline material. Whereas attenuation can in principle be predicted from the known distribution of the attenuation coefficient (measured simultaneously with the diffraction signal in a DCT experiment). In many cases the influence of attenuation can be reduced by selecting a X-ray energy for which photoelectric absorption in the sample becomes negligible.

On the other hand, the possibilities to predict or reduce the influence of extinction effects is more problematic. In general two sources of extinction have to be considered: primary extinction (multiple diffraction) within the grain of interest leads to deviation of diffracted intensity from our model assumption of a linear relationship between crystal thickness and diffracted intensity. Fortunately, deformed materials are known to develop a mosaic structure which limits the influence of this effect encountered in rather perfect crystals. In polycrystalline samples another important source of noise is secondary extinction, referring to local changes of the intensity of the diffracted beams, caused by attenuation of the incoming and / or diffracted beams by other grains which happen to fulfill the Bragg condition on the trajectory through the sample volume (figure 15). Knowing the 3D grain structure of the material, this phenomenon could in principle be predicted, as well. However, the number of grains and diffraction events to be taken into account has prevented so far any attempts to explicitly take secondary extinction into account. Due to the randomized nature of this phenomenon, the behaviour of the algorithm in respect to this kind of inconsistencies can be analysed by adding synthetic extinction noise to the projection data.

We decided to model extinction noise as a reduction of intensity in spherical regions of the blobs. The position and the radius of these spheres has also been

B.2 Reconstructions of noisy blobs

We now show the results of reconstructions performed using blobs affected by the extinction noise. In the reconstruction shown in figure (16, b), all the blobs were affected by a random value of extinction noise, like in figure (14, c). Looking to figure (16), compared to the same reconstruction without noise (fig. 5), we observe only a moderate increase of the reconstruction error, which seems to indicate that the lower limits of the error in this test case are rather dominated by the lack in resolution of the orientation sampling and other deficiencies of the model.

Chapter 4

Reconstruction of grain clusters

This chapter will be about a recently accepted article on Scientific reports, that in our opinion is going to significantly change the future of DCT, and it is among the most important result in this thesis. The concepts discussed in chapter 3 have been further refined and, most importantly, applied to experimental data, and showing very convincing results, especially because they were compared to measurements performed using EBSD.

The notable change is that the six-dimensional framework is not only applied to the single grain reconstruction as in traditional DCT fashion, but it has also been applied to a cluster of grains. We call a certain region of the sample a "cluster of grains", when within this region all the adjacent grains have orientations that are confined in a small region of orientation space, and they could be viewed as sub-grains of a bigger crystallographic domain, resulting for instance from the process recovery of a somewhat deformed grain. Clustered regions represent one of the biggest challenges for DCT reconstructions because the diffraction spots from the grains in these regions are likely to cluster on the detector as well, giving rise to overlaps for some of the reflections and being detached for some other reflections. Overlaps are very detrimental for the reconstruction of the single grains, because if it is not possible to decide whether a spot is the result of an overlap or not, it is not even possible to correctly determine its center of mass, and all the subsequent steps in the DCT analysis toolchain get broken.

The method used to reconstruct the clustered regions was conceptually simple, yet very effective. Given the correct indexation of some sub-grains with similar orientation, in a region that presents some holes, and so unindexed grains, it is very likely that the missing regions will belong to other sub-grains that have similar orientations to the surroundings. The simplest idea here is to take a bounding box in real-space that accommodate for the whole region in sample, a bounding box in orientation-space that accommodates all the indexed grains, CHAPTER 4. RECONSTRUCTION OF GRAIN CLUSTERS oversize them slightly, and define the combined six-dimensional bounding box as the region to reconstruct. At this point, a total of 64 vertexes in the combined six-dimensional space can be calculated, given by the outer product between the 8 vertexes of the real-space bounding box and the 8 vertexes of the orientationspace bounding box. It is possible to predict on the detector the (θ, η, ω) bounding box for each reflection that hits the detector. The reconstruction of the clustered regions in this article took place from the raw images, without any kind of mask, to not remove unidentified spots, or to not alter the projection data in any way. It is interesting to note in article's figure (2), that even if the raw images were used, where we can see overlaps from reflections of grains from other regions of the sample, the algorithm was able to discriminate between the signal coming from the to-be-reconstructed region and the noise. However the most interesting figure is probably the render of isosurfaces in the orientation-space from article's figure (4, a) where some previously unindexed regions were definitely reconstructed by the six-dimensional algorithm and the corresponding regions in real-space were now filled, as it can be seen from the comparison of figure (1, d) with figure (1, a). This result also inspired the six-dimensional bounding box extension approach that was not discussed in the paper, to not deviate from the main message. In fact, it can be seen from the comparison of article's figures (1, a) and (1, d), that the orange region in the very lower part of the sample surface, underwent a noticeable improvement, but that region could hardly be called a cluster, or work as a cluster like region. Figure4.1 shows that the region now under consideration was indeed a grain, exhibiting a sub-grain boundary that would trick first the segmentation phase into thinking that some of spots of the grain were split into different spots, and then the indexing and forward simulation would not be able to assign the correct spots to the grain. As a consequence, especially from the comparison of figures 4.1a and 4.1b, it becomes clear that applying an extension of the predicted orientation-space bounding box, and using the raw images from the detector, would allow to catch all the regions of the considered grain. This "extended grain" approach is very similar to the "cluster of grains" approach and in fact there is no difference in the algorithm for the analysis of the two types of data. The consequences of these findings are outlined in the attached article, and clearly change the possible future of DCT. By being able to circumvent indexing's shortcomings when dealing with deformed materials, this would open the possibility to attack new types of problems with this technique. Moreover, DCT is based on a fast acquisition scheme, and using the six-dimensional model, we could now think of applying the enhanced 6D-DCT to new scenarios where progressive deformation of the sample is applied.

February 5, 2016

Abstract A previously introduced mathematical framework for full-field X-ray orientation microscopy is for the first time applied to experimental nearfield diffraction data acquired from a polycrystalline sample. Grain by grain tomographic reconstructions using convex optimization and prior knowledge are carried out in a six-dimensional representation of positionorientation space, used for modelling the inverse problem of X-ray orientation imaging. From the 6D reconstruction output we derive 3D orientation maps, which are then assembled into a common sample volume. The obtained 3D orientation map is compared to an EBSD surface map and local misorientations, as well as remaining discrepancies in grain boundary positions are quantified. The new approach replaces the single orientation reconstruction scheme behind X-ray diffraction contrast tomography and extends the applicability of this diffraction imaging technique to material micro-structures exhibiting sub-grains and/or intra-granular orientation spreads of up to a few degrees. As demonstrated on textured sub-regions of the sample, the new framework can be extended to operate on experimental raw data, thereby bypassing the concept of orientation indexation based on diffraction spot peak positions. This new method enables fast, three-dimensional characterization with isotropic spatial resolution, suitable for time-lapse observations of grain microstructures evolving as a function of applied strain or temperature.

Introduction

Their ability to map crystallographic phase and orientation in three dimensions make X-ray diffraction imaging techniques a highly complementary tool to established 2D electron microscopy techniques like electron backscatter diffraction (EBSD). In this article we present first experimental results obtained with a new sixdimensional data analysis framework [22], applicable to monochromatic beam, near-field X-ray diffraction data, aiming at the reconstruction of spatially resolved, three-dimensional orientation maps. As opposed to 3D X-Ray Diffraction (3DXRD) variants based on pencil beam (1D) or line beam (2D) illumination which involve 3D or 2D scanning procedures limiting the ultimate acquisition speed, the full-field variant described in this paper requires only one rotational scan and can be executed on standard X-ray tomography instruments, offering simultaneous X-ray absorption and phase contrast inspection of the illuminated sample volume. Depending on the experiment settings and material characteristics, the reconstruction task related to X-ray orientation microscopy [15] can be sub-divided into one of the following three categories: (1) negligible intra-granular orientation spreads (2) presence of intra-granular orientation distributions which can be described as a 3D vector field and (3) presence of intra-granular orientation distributions which require a description as 3D orientation distributions per sampled volume element. In the first case the task of 3D grain shape reconstruction reduces to the classical problem of image reconstruction from parallel projections and can be solved using algorithms developed in the field of medical imaging ( [5], [6], [15], [13], [12], [16]). Common to this class of micro-structure reconstruction techniques is a two step process based on poly-crystal indexing followed by grain by grain reconstruction. In the indexing step grain orientations and positions are identified based on a systematic search through the list of scattering vectors derived from the measured diffraction spot peak positions. In the second step, a projection stack and the corresponding projection geometry are constructed for each of the grains, which are reconstructed individually and assembled into the common sample volume. However, it has to be noted that the concept of indexing grain orientations from diffraction spot peak positions applies to materials and acquisition conditions giving rise to limited diffraction spot overlap on the detector. This concept is known to gradually fail with increasing complexity of the microstructure and macroscopic plastic deformation (giving rise to increasing values of intragranular orientation spread and lattice strain) and eventually results in non-spacefilling grain maps (see Supplemental Information for a more detailed discussion of the interplay and the typical values of material and acquisition parameters determining the applicability of polycrystal indexing approaches).

Alternative strategies for X-ray orientation microscopy replace the above-mentioned two step process with a forward modeling approach performing a voxel-wise optimization in order to maximize the overlap between experimental data and the current description of micro-structure. These strategies are well suited for the second class of reconstruction task (i.e. description of the material microstructure in terms of a 3D orientation field) and have proven to be applicable to samples which have undergone significant plastic deformation [10]. Prominent examples are the algorithms entitled 'Grain Sweeper' [17,2] and IceNine [20,11], established in the frame of of the High Energy Diffraction Microscopy (HEDM) project at the Advanced Photon Source. Both algorithms share the concept of local (voxelwise) optimization and operate on a binarized version of the diffraction data. In this article we introduce a framework which deals with the third, and most complicated case, accounting for a 3D orientation distribution for each sampled volume element. The formulation of the reconstruction task in six-dimensional position-orientation space was first proposed in earlier work by Poulsen [14], but an actual implementation became only possible recently, thanks to progress in computing hardware and emergence of new concepts and algorithms in the field of mathematical optimization. Introducing a regular sampling over the sub-volume of 3D orientation space occupied by a grain, we have shown that the inverse and underdetermined problem of X-ray orientation microscopy can be solved with the help of first-order, convex optimization algorithms, using physical meaningful regularization terms [22]. Whereas this previous work was based on synthetic diffraction data we present here a first application to experimental data acquired from a polycrystalline halite (rock-salt) sample containing more than 300 indexed grain orientations in the illuminated sample volume. We demonstrate that the quality and applicability of fast, full-field acquisition techniques can be significantly enhanced when switching to such a six-dimensional formulation of the reconstruction problem. Three different types of grain maps, all reconstructed from the same set of experimental data are presented and compared to each other in this work. The first two maps are based on the two step process of indexing grain orientations, followed by algebraic reconstruction using the conventional 3D (single orientation) reconstruction algorithm [12,16] and the new six-dimensional framework [22], respectively. The material under study presents a pronounced grain substructure and some of these sub-grains could not be identified using conventional indexing schemes based on diffraction spot center of mass positions. We then present an extension of the six-dimensional approach, operating on experimental raw data and circumventing the steps of diffraction spot segmentation and orientation indexing. Seeded with approximate information concerning the real space and orientation space sub-volumes of the region to be reconstructed, this extended approach correctly identifies the missing sub-grain orientations and results in a space-filling grain map. For visualization and evaluation of the reconstruction results we used the reduced, 3D vector field representation of the orientation field, derived from the 6D reconstruction output. The orientation of a real space voxel in this repre-sentation is calculated as the center of mass of the local orientation distribution associated to this voxel. The resulting 3D orientation maps have been compared to EBSD measurements acquired on one of the lateral sample surfaces and local disorientation and discrepancies in the grain boundary positions have been analysed.

Method

In a monochromatic beam, near-field diffraction experiment the sample is positioned on a rotation stage, and as it rotates by angle ω, it gives rise to diffracted beams each time the Bragg condition is met for one of the grains. A limited number of the diffraction spots will intersect the high resolution imaging detector positioned a few millimeters downstream of the sample. Each of the sub-orientations present in a grain is associated to a slightly different projection geometry and the diffraction signal associated to a given Bragg reflection is observed as a three-dimensional diffraction "blob" volume, parametrized by two spatial coordinates u and v (detector pixel coordinates) and a rotation angle ω (image number). Whereas in the previous implementations of DCT the reconstruction units were grains, described by an average orientation and a corresponding projection geometry applied to the integrated diffraction spots (2D), we now introduce an explicit, discrete sampling of the local orientation space centered around the grain average orientation and exploit the additional information encoded in the intensity variation of the 3D diffraction blob volumes as a function of the rotation angle. We neglect the possible (in metals typically ≤ 1%) elastic distortion of the crystal unit cell and introduce a 6-dimensional reconstruction space X 6 = R 3 ⊗ O 3 obtained by the outer product of cartesian real-space and 3-dimensional orientation space. Assuming kinematic diffraction and neglecting photoelectric absorption and extinction effects, the process of diffraction image formation (forward projection) can be formulated as a linear operator:

Ax = b (1)
where x is a vector containing N P elements, representing the scalar "scattering power" for each of the sampled volume elements in the six-dimensional position -orientation space (N P = n 3 × p 3 for the case of a regular sampling over cubeshaped sub-volumes with grid length n and p in position and orientation space, respectively). Each line of the projection matrix A S×N P contains the contribution of the 6D volume elements to a given detector pixel and the vector b S holds a list of measured pixel intensities, specified by their (u, v, ω) coordinates in the 3D diffraction image stack. S corresponds to the total number of detector pixels in the 3D image stack reached by the M projections (hkl reflections) of the grain volume(s). If the elements of x are arranged as a succession of 3D real space volumes, each representing one of the sampled orientations, and A as an array of M × P projection matrices, each one describing the projection geometry for one of the M hkl reflections intercepted by the detector, the vector b will be composed of M blocks corresponding to contiguous 3D subvolumes (i.e. the aforementioned diffraction blobs), spread throughout the entire stack of detector images:

   A 11 • • • A 1P . . . . . . . . . A M 1 • • • A M P       x 1 . . . x P    =    b 1 . . . b M    (2)
(Removing lines of b i with elements corresponding to ω positions not reached by current orientation j, the resulting sub-matrices Ãij x j = bi correspond to parallel beam projections of single orientation volumes onto 2D diffraction spot images bi , as used in previous 3D implementations of DCT [12].)

Formulation as a 6D optimization problem and underlying assumptions

While equation ( 2) is a concise formulation of the diffraction image formation (forward projection), it is not suited for the actual reconstruction. Whereas in the case of 3D reconstructions based on single (grain average) orientation one had to solve for one unknown (local scattering power) per real space volume element, we now have to solve for a much larger number of unknowns (i.e the scalar scattering power for each of the discretized orientations (∼ p 3 -typically several hundreds up to a few thousands of orientations per real space voxel). At the same time the number of measurements typically only increases by a factor of p [22]. Moreover measurements from a polycrystalline sample are affected by diffraction spot overlaps from other grains and various sources of noise. As a consequence, equation ( 2) is an ill-posed, heavily under-determined problem and has to be re-written as a minimisation problem, in which case additional terms enforcing physical meaningful constraints and prior knowledge about the solution can be incorporated in its optimization functional. While in previous work on synthetic diffraction data [22] we chose to use an l 1 minimization over the full six-dimensional space, because of the a priori knowledge that the phantom was based on a (sparse) 3D vector field representation of the grain volume, in the current work we used a slightly different formulation (equations 3 and 4) which promotes positive solutions minimising spatial variations of the integrated scattering power assigned to the real space voxels, as to be expected from a purely kinematic diffraction model.

x * = argmin x ||Ax -b|| 2 2 + λ|| (|∇Sx|) || 1 (3)
subject to:

x ≥ 0 ( 4
)
where S is the operator that sums all the orientation components for each realspace voxel, and the l 1 -norm over the absolute value of the gradient is the total variation operator ([3], [18]).

This formulation is well adapted to the reconstruction of materials showing a pronounced sub-grain structure, since cross-talk between the sampled orientations tends to assign intensity to the grain boundary regions and results in sharper boundaries and better homogeneity of the reconstructions compared to the previous l 1 -norm formulation. Like in the previous work [22], a specific class of algorithms called Chambolle-Pock [4] is used for the minimization of the functional in equation 3. As shown in the supplementary material this functional can be rewritten in a different form, leading to a very similar algorithm, but now with a fixed weight λ = 1, more suitable for practical work on large dataset containing a large number of individual reconstruction problems.

For more details on the calculation of the projection matrix, on the 6-dimensional space and on the optimization algorithm, we refer to the supplementary material and [22] [1] [15].

Results

A parallelepiped sample with 0.9 × 0.9 × 2.2 mm dimension was prepared from a two-phase materials consisting of large (100-400 µm) NaCl grains and a fine dispersion of micrometer sized Cu particles. (The fine dispersion of Cu serves as a contrast agent for digital volume correlation and will be used for determination of the 3D displacement fields introduced by incremental compressive loading of the specimen, foreseen in a follow-up experiment.) Halite is a viscoplastic geomaterial deforming under the action of grain boundary sliding and dislocation mediated plasticity [19,24]. Here we report on the characterization of the 3D orientation field prior to plastic deformation of the specimen. Closer inspection of the EBSD surface mapping (figure 1a) reveals that part of the bigger grains consist of smaller sub-grains with typical dimensions of order of 100 µm and which are delineated by small angle boundaries with a few degrees of misorientation. A 560 µm high sub-volume of the specimen containing more than 300 grains was scanned on a conventional X-ray imaging setup available at the bending magnet beamline BM05 of the European Synchrotron. After the synchrotron experiment one of the lateral sample surfaces was characterized by electron backscatter diffraction for comparison and cross-validation of the X-ray orientation maps calculated from the X-ray diffraction data. More details on the sample preparation and experiment conditions used for the DCT and EBSD characterization are given in the supplementary material.

DCT reconstructions

Three different types of grain reconstructions were performed in order to illustrate the improvements of the new reconstruction framework compared to the previous (single orientation) approach: (1) a standard 3D DCT reconstruction of the grains that have been identified in the polycrystal indexing step. Integrated 2D diffraction spot images were used for the reconstruction and the resulting 3D grain map underwent a two voxels dilation step as described in [12] (2) a 6D DCT reconstruction of the same set of grains, but now reconstructed from the 3D diffraction blob volumes using a variant to of the recently introduced 6-dimensional algorithm [22] and (3) an extension of the latter where sub-regions of the 3D orientation map corresponding to some extended regions in real space and orientation space have been reconstructed from unsegmented experimental raw data.

As mentioned earlier, the 6D reconstruction output has been reduced to the conventional 3D vector map representation, commonly used in orientation imaging microscopy. In all three cases the final shape of the 3D grain volumes is determined by an automated thresholding operation at the time of assembling the individual grain reconstructions into the common sample volume. Note that no additional dilational postprocessing was applied to the grain orientation maps obtained by the six-dimensional reconstruction approaches.

Grainmaps comparison

Figure 1 shows a comparison of the EBSD map acquired close to one of the sample surfaces with the corresponding section through the reconstructed grain volume obtained with the three approaches described in the previous section. The grain maps are color coded according to the inverse pole figure mapping along the surface normal and overlaid on the X-ray attenuation map (as explained in Supplementary Information the DCT volume and orientation matrices have been rotated by about 2.4 degrees in order to coincide with the reference frame of the backscatter electron image. The IPF colour key is provided in figure S1.) The absorption image has been reconstructed from the simultaneously acquired X-ray projection images and it is intrincically aligned with the DCT grain map. Due to their higher attenuation coefficent the Cu particles show up as bright particles, decorating part of the high angle grain boundaries in the NaCl matrix material. The construction of the depicted X-ray orientation maps involves the steps of segmentation and assembly of individually reconstructed grain volumes into the common sample volume. This process may lead to non-space-filling orientation maps in case of grain orientations not been identified in the previous indexing step. Inspection of the EBSD map indicates that the "holes" in proximity of the orange and rose grains at the bottom and right side of the depicted grain maps figures (1b,1c) correspond to sub-grains separated by low-angle boundaries with less than 3 • misorientation. Since the connectivity of diffraction blobs originating from grains with a pronounced sub-structure changes as a function of the (hkl) reflection, the center of mass and shape based indexing procedures outlined in [12,16] may fail to identify part of the existing sub-grains in this situation. However, as long as one of the sub-grains has been identified a straightforward extension of the six-dimensional reconstruction approach discussed in the next section can be used to find the others and eventually results in a space-filling grain map (figure 1d) with a much reduced number of ambiguously or unassigned voxels.

As the attenuation map is simultaneously measured with the 3D orientation map, with the same detector, in the same reference system, it is possible to complement the grain map with the information about the position of the Cu particles. This type of information is also interesting because the Cu particles are only observed at the grain boundary positions, and their presence is both a direct confirmation of a correct indexing and a measure of the spatial accuracy of the reconstructed grain maps.

Reconstruction of sub-grain clusters and intra-granular orientation

Having identified one of the sub-grain orientations in a textured and/or badly reconstructed region of the sample a straight forward extension of the current framework allows to identify and reconstruct neighboring sub-grains without the need to explicitly identify and isolate the corresponding diffraction blobs on the detector. It is enough to extend the real-space bounding box to include the missing region, and to (iteratively) extend the orientation-space bounding box until the missing sub-grains are fully included. A new stack of diffraction blobs (difstack as defined in [22]), that would include the region on the detector covered by the forward projection of the extended real-space and orientation-space volumes is assembled and then reconstructed with the very same 6D algorithm used for the indexed grains. There is, however, one important difference: unlike the difstack of indexed grains which is constructed from a subset of pre-selected and segmented diffraction blob volumes, the difstack corresponding to a clustered region is directly assembled from the background corrected experimental raw images. This in turn may lead to considerable degree of overlap between the reflections of the grain under analysis and spurious reflections from other grains, as depicted in figure (2a). In fact, these spurious overlaps are mathematically inconsistent to each other, from one blob to another one, and the reconstruction algorithm is able to greatly reduce their impact on the final reconstruction. Indeed, figure 3c), and the overlaps, affecting about 50 % of the projections, are completely filtered out. It should be noted that even with indexed spots overlaps may occur, so this tolerance is important.

Figure (3a) is a zoom on the clustered region in the right mid-bottom part of figure(1d). It shows that the intra-granular misorientation is characterized by a maximum value of 2.5 degrees misorientation from the chosen reference orientation and illustrates the local orientation characterization capability of the new six-dimensional reconstruction framework. Each of the real space voxels carries a local orientation determined from the weighted average of the intensity assigned to the 3D orientation space voxels related to this specific position. The calculation of the kernel average misorientation in figure (3b) clearly reveals the presence of small-angle grain boundaries in this region. Note that these boundaries are not sharp, which can be partly attributed to the fact that the orientation-space resolution for this calculation was limited to a grid of 11×11× 11 orientations over a bounding box of 3.1 • × 2.5 • × 3.9 • , because of memory constrains. 123

122

Interesting insight about the structure of the solution is also gained from inspection of the local 3D orientation space associated to this sub-structured region of the sample volume. Figure (3d) confirms that the algorithm was able to identify a few new sub-grain orientations in addition to the five indexed sub-grain orientations detected by the conventional indexing procedure. The application of this extended reconstruction framework to other regions which were incompletely reconstructed in figure (1c) results in a significant improvement of the reconstruction which now is close to space-filling (figure 1d). Finally, figure (4), shows a quantitative comparison between the EBSD map in figure (1a) and the 6D-DCT map in figure (1d). More precisely, figure (4a) shows the voxel-by-voxel misorientation between the two maps. Regions corresponding to Cu particles and voxels close the grain boundaries, assigned to different grains are shown in white and have been excluded from this analysis. The histogram of the misorientation distribution (figure (4c) shows a pronounced peak at 0.2 degrees and some low intensity tails up about 1 degree misorientation. In figure (4b) instead, an overlay of the grain boundaries computed with the two techniques is presented. The local distance between the two maps has been determined by intersecting the (binary) DCT boundary map with the distance transform calculated from the EBSD boundary map. The corresponding histogram with an average distance of 1.06 voxels (corresponding to 3.7 µm) is presented in figure (4d).

Discussion

With the transition from a 3D to a 6D reconstruction framework, the previous restriction of diffraction contrast tomography to materials with negligible values of intragranular orientation spread disappears. By assigning an average orientation to each of the reconstructed voxels, the output becomes identical to well established 2D and 3D grain mapping techniques based on scanning electron or X-ray diffraction techniques. However, an important difference remains: unlike scanning techniques or reconstruction schemes based on voxel by voxel optimization, the global optimization approach as presented in this article was seeded with some approximate information concerning grain position and orientation, obtained from a previous polycrystal indexing step. This implies a remaining restriction to moderately deformed microstructures (typically ≤ 5% plastic deformation), which can be still be described as a collection of grains, representing 3D crystal domains with limited orientation spread and well-defined boundaries (see Supplemental Information for a discussion of the requirements for successful orientation indexation from diffraction spot peak positions). As expected from a previous feasibility study on synthetic diffraction data, the transition from the previous single-orientation (3D) to a 6D reconstruction framework results in a visible improvement in the accuracy of the grain shape reconstructions. The extended model can account for the non-parallel projection geometry of deformed grains and the iterative reconstruction process produces consistent grain shapes which can be assembled into the 3D sample 124 volume with much reduced tendency for overlaps (i.e. voxels simulateneously claimed by adjacent grains) and reduced gaps between grains, corresponding to unassigned voxels. Indeed figures confirm that the grain boundaries obtained with the 6D-DCT approach (using the "grain cluster" generalization for the clustered regions and no dilational postprocessing of the grain map), are in good agreement with the EBSD measurements. The accuracy of the grain boundary positions determined with the 6D method appears to approch the voxel size of the reconstruction (3.75µm in this study). It can be expected that values down to 1 µm resolution can be reached when working with high resolution configurations of the detector system, as typically employed in studies of materials with grain sizes in the range of 10-50 µm.

For what concerns the reconstruction of local orientations, figures (4a) and (4c) suggest that also in this case the algorithm is able to retrieve an average orientation for each voxel that is close (i.e. within the accuracy of the EBSD measurement) to the results obtained from EBSD. The full width at half maximum of the misorientation distribution (0.25 degree) is comparable to the expected accuracy of the EBSD measurement. As explained in the Supplementary Information, the overall shift of about 0.2 degree in misorientation can be explained by remaining uncertainties in the alignment of the real-space and orientation space reference systems of both measurement modalities. We furthermore recall that the results presented in figures (1c) and (1d) are a projection of the full six-dimensional reconstruction output into the 3D vector field representation that describes the local average orientations in each voxel. This means that by using this representation we have lost the information concerning the possible presence of sub-voxel orientation domains. Further tests are required in order to evaluate the accuracy of the local (per voxel) 3D orientation distribution output by our optimization routine. Note that the 2D EBSD measurements presented in this study are not sufficient for giving a ground truth for the six-dimensional reconstruction, since the interaction volume and information depth (of order of 100 nm) of EBSD is much smaller than the voxel size (3.75µm) used in this study. Three-dimensional, high angular resolution 3D EBSD on dual-beam FIB-SEM instruments [23] or Differential Aperture X-ray Microscopy (DAXM [9]) measurements shall be considered for this task. The current model is based on an idealized (purely kinematic) description of the diffraction process and calculations are performed on diffraction blob volumes which have been renormalized to constant intensity, which is a pragmatic but rather crude way to account for intensity variations due to differences in structure and Lorenz factors, spatial and temporal inhomogeneity of the incoming beam profile and attenuation of the incoming and diffracted beam (the latter also known as self-absorption) due to photoelectric absorption and extinction effects. The model could in principle be extended to account for these effects: structure and Lorenz factors can be calculated and the incoming beam profile is periodically updated by taking images without the object during the tomographic acquisition procedure. Moreover, the 3D attenuation coefficient distribution is reconstructed from the transmitted beam and after a first 3D reconstruction one could in principle check if grains on the trajectory of the incoming and diffracted beam simultaneously fulfill the diffraction condition, giving rise to local intensity variations. An obvious extension of the approach applied to the sub-grain clusters would consist in enlarging the sampling of the orientation space to the entire fundamental zone. Covering the fundamental zone for high symmetry space groups like the cubic with an orientation resolution of 0.2 degrees would require of order of 10 8 degrees of freedom. With currently available computing hardware featuring up to 512 GB of memory, the model could possibly handle high-resolution diffraction data acquired in slice beam illumination mode.

Another promising development route would consist in simultaneous acquisition of full field (extended beam) diffraction data on 3D detector systems featuring two or more screens with different pixel size and positioned at different distances. The reduced number of spatial degrees of freedom on the outermost low spatial resolution (diffraction) detector would allow the identification of orientations and approximate positions [17,7] which in turn can be reconstructed at higher spatial resolution from the near-field diffraction data, using the extended approach presented in the current paper.

The design and implementation of the reconstruction algorithm presented here is also independent from the type of tomographic projector used and could be easily modularized to allow for polychromatic projectors [21] to be used in both polychromatic X-ray or neutron imaging applications [8]. The advantage of using such framework over the previous approaches would be once again in the possibility to explicitly take the distortion of diffraction images due to the local variations of the diffraction angles into account. Due to additional variation of the Bragg angle, these distortions are known to be more severe in polychromatic imaging [8] and severely restrict the choice of materials which can be analysed in the single orientation framework. As a final remark, it should be stated that the current work is based on a 6dimensional kinematic diffraction model which does not take into account the possible distortion of the unit cell as result of elastic strains, which would require additional six degrees of freedom. Due to the exponential growth of memory requirements with the number of sampled dimensions the discrete sampling approach used in this work will not be an appropriate framework for addressing the full 12-dimensional problem.

Conclusions

We have demonstrated the feasibility of three-dimensional full-field X-ray orientation microscopy from extended beam near-field X-ray diffraction data. A new six-dimensional reconstruction framework yields spatially resolved 3D attenuation and orientation maps, substituting and outperforming the previously introduced three-dimensional (single orientation) reconstruction algorithm behind X-ray diffraction contrast tomography.

The results obtained on a two phase material made from NaCl and containing a fine dispersion of Cu particles have been cross-validated against EBSD measurements on the sample surface and indicate that our approach yields an orientation resolution comparable to EBSD and a 3D spatial resolution consistent with the pixel size of high resolution X-ray imaging detectors.

As demonstrated on a material displaying a pronounced grain sub-structure, the introduced six-dimensional frame can be extended to operate on non-segmented raw data corresponding to sub-volumes of the six-dimensional position -orientation space. This finding indicates possible future extensions of the framework, such as replacing the two-step process of orientation indexing and reconstruction by a global optimization procedure by including extended regions of orientation space into the reconstruction process.

Finally, as a full-field approach, compatible with simultaneous absorption and phase contrast imaging, the proposed methodology is currently about one order of magnitude faster than competing techniques based on slice or pencil beam illumination. This specific feature opens interesting new possibilities for timelapse observations of processes like plastic deformation, coarsening or phase transformations in bulk polycrystalline structural materials.

Three-dimensional full-field X-ray orientation microscopy Nicola Viganò, Alexandre Tanguy, Simon Hallais, Alexandre Dimanov Michel Bornert Kees Joost Batenburg, Wolfgang Ludwig Supplementary Information 1 Experimental Details 1.1 X-ray diffraction contrast tomography data acquisition X-ray diffraction contrast tomography measurements were performed at the bending magnet beamline BM05 of the European Synchrotron Radiation Facility (ESRF). The central part of a parallelepiped sample was illuminated by a monochromatic X-ray beam (35 keV; 1.5 (h) × 0.56 (v) mm cross-section), delivered by a Si 111 double crystal monochromator. The transmitted and diffracted beams were captured on a detector system positioned 15 mm downstream of the sample position. The detector system consisted of a back-illuminated Frelon CCD camera, coupled via visible light optics to a 100 µm thick transparent luminiscent screen made of LuAG [5] and provided an effective pixel size of 3.75 µm and a field of view of 7.7 × 7.7 mm. 3600 images with an exposure time of 8 sec were recorded during a 360 • continuous rotation movement of the sample. With these settings the five inner-most Debey-Scherrer rings were intercepted by the detector and on average 45 out of the 60 diffraction blobs could be unambiguously assigned to one of the 345 indexed grains in the illuminated sample volume.

EBSD acquisition

In order to compare the 6D-DCT reconstructions to a classical EBSD surface measurement, the surface of the sample was further prepared by broad argon ion beam sputtering, using the most recent ion-beam polishing system from Gatan [PECS-II], which combines cross sectioning and planar polishing capabilities. In order to prepare a large mm-sized area, we used the planar polishing technique, applying an acceleration voltage of 2 kV and 5 • beam incidence, for a duration of about 1 hour.

The EBSD data were acquired on an environmental Quanta 600 FEG-SEM 1 132 instrument, operated at 20 Hz acquisition rate, an acceleration voltage of 20 kV , a working distance of 17.6 mm and a residual air pressure of 100 Pa in order to avoid surface charging. A grid of 5 × 4 partially overlapping EBSD maps, of 300×270 µm size each, allowed for covering an entire lateral sample surface with a step size of 1.3 µm. The acquired diffraction patterns were analyzed using the AZtec HKL system from Oxford Instruments. The standard procedure with a minimum of 12 indexed Kikuchi lines allowed for 99.6% indexation success. A surface layer of approximately 20 µm in depth was removed in five successive ion-sputtering steps before a full correspondence between the micro-structures observed by EBSD and DCT was obtained.

IPF coloring

The coloring key used for figure (1) in the text of the article is reported in figure (1).

Registration of datasets and remaining sources of error

The two-dimensional EBSD mapping and the three-dimensional tomographic reconstruction have been carefully aligned with respect to each other, using an optimization routine based on the position of a set of Cu particles identified in both imaging modalities. Cu particles are visible in the 3D absorption image and in the backscatter electron image acquired in conjunction with the EBSD characterization. The optimization provides a best fit between the X-ray absorption map and the electron backsacttered electron image of the polished sample surface. Nevertheless, the histograms in figure 4 show some remaining discrepancies.

133

Whereas part of these discrepancies may be attributed to actual errors introduced by the reconstruction and segmentation procedure (due to uncertainties in the experiment geometry, non-idealities of the detector system and various other sources of noise), a series of other factors complicate a one to one comparison between both imaging modalities:

1. The voxel size of the DCT map (3.75 µm in our case) is significantly bigger than the typical interaction volume of EBSD measurements (of order of 100 nm). Even though the EBSD maps have been sampled with a smaller step size of 1.3µm, the volumes from which the average orientations in both maps are derived from are not identical.

2. The comparison is carried out between a not ideally flat surface (there is some remaining curvature and surface relief created by the ion beam polishing procedure) and a geometric plane inside a 3D volume.

3. There may be some remaining spatial distortion of the (stitched) EBSD map after the affine transformation (70 degree tilt) back into the coordinate system aligned with the electron beam direction.

4. The alignment of the orientation-space reference system between the two datasets does not perfectly coincide with alignment of the real-space reference systems: while the real-space alignment suggested a rotation of 2.6 degrees, the orientation-space alignment suggested a rotation of only 2.4 degrees. This is however is in-line with the observation of a shift in the peak of misorientation histogram shown in figure 4d (main text). We attribute this to a remaining curvature of the sample surface and small errors in the calibration between the reference frames of the EBSD map and the back-scattered electron image.

5. The distance transform calculated in 2D can only provide a conservative estimate of the true distance in three-dimensions (boundaries cutting the observation plane at a small angle, can give rise to large errors).

Sample preparation

This work is part of a research project aiming at the characterization of plastic deformation mechanism in Halite (rock-salt, NaCl) by combined use of digital volume correlation (DVC) [2,4] and 3D grain mapping techniques. Abandoned rock-salt mines are currently considered as natural reservoirs for compressed air energy storage power plants, and there are questions concerning creep deformation and damage modes of salt as a geo-material [1].

In order to enable the characterization of 3D displacement fields in the bulk of the polycrystalline microstructure as a function of applied plastic deformation a two-phase material providing internal (absorption) contrast was prepared for this study. Pure synthetic NaCl powder and a fine grained dispersion of Cu particles (3 volume %, average size 3 µm) were mechanically mixed at room temperature. A nearly dense cylindrical block of 7 cm in diameter, with about 3 % of porosity and widely spread NaCl grain size distribution (20 -300 µm) was obtained after a first stage of hot isostatic pressing (HIP) at a temperature of 150 • C and 100 MPa for 168 h. In order i) to promote NaCl grain growth to a size compatible with DCT characterization of a mm-sized specimen, and ii) to reduce the residual porosity for the purpose of the mechanical testing, additional high temperature HIPing of cm-sized cylinders was performed at 700 • C and 2.5 MPa over 24 hours. The latter procedure resulted in a homogeneous two-phase micro-structure, consisting of equiaxed NaCl grains with narrow size distribution (100 -400 µm) and an average size of about 250 µm (as determined by the linear intercept method, considering grain boundaries with ≥ 10 degrees misorientation ) and a residual porosity of less than 1 vol%. Interfacial dragging of the copper particles during NaCl grain growth did result in a micro-structure presenting individual NaCl grains and clusters of a few NaCl grains delineated by a more or less continuous 3D network of segregated Cu particles. A parallelepiped shaped sample was extracted from the center of the hot pressed material using a low speed diamond saw, and was carefully polished down to grid 4000 SiC paper to final dimensions of 0.9 × 0.9 × 2.2 mm. Even if the material was tailored for this type of combined study, the microstructure of the material shows very similar characteristics to the microstructures encountered in stuctural materials like metals and their alloys. In both cases the material can deform through dislocation mediated plasticity, and dislocations may organize in sub-structures leading to the formation of cells and small-angle grain boundaries.

Sample requirements

Although the full-field diffraction approach presented in this work can deliver comparable output to other well established orientation imaging techniques, the higher convolution of the diffraction signal (illumination of the full 3D sample volume) imposes more stringent conditions on the type of microstructures which can be successfully analysed, as compared to 2D or 3D scanning techniques.

Diffraction spot overlap One of the principal limitations is related to the concept of orientation indexing based on a systematic search through scattering vectors, derived from diffraction spot peak positions. Acquisition conditions resulting in diffraction spot overlap will bring this concept to fail, because the erroneous determination of the centers of mass will result in erroneous determination of the associated scattering vectors and poor indexation results. The probability of diffraction spot overlap increases with the number of grains in the illuminated sample volume, the intragranular orientation spread, and the strength of texture. While it is difficult to provide absolute numbers for the individual conditions, due to the interplay of these different factors, it is possible to alleviate the probability of overlap by reducing the number of grains (by reducing the 4 135 sample diameter and/or reducing the dimensions of the illuminated sample volume by closing the gap of the beam-defining slits), and optimizing the distance between sample and detector.

Number of grains

In the absence of deformation and texture, peak search based indexing techniques can work with up to a few thousand grains, while for a 5% deformed samples, a more realistic number is of order hundred.

Sample size

The sample dimensions in the direction perpendicular to the rotation axis should in general not exceed 20 grain diameters. Deformed and textures materials might require further reduction to values below 10 grain diameters.

Grain size The minimum grain size limit is closely linked to the pixel size of the high resolution X-ray detector system employed. In order to spatially resolve 3D grain shapes, one should aim at about 10 voxels across the typical dimensions of the smallest grain size to be reconstructed. Given the physical limitation of high resolution X-ray detector systems to about 0.5 µm resolution in the best case, the minimum grain size which can be handeled with type of approach is of order of 5-10 µm.

Scaling Provided the material has a mono-modal grain size distribution, the pixel size can typically be adjusted over a wide range (0.5 µm (in U O 2 ) up to 30 µm (in ice) have been demonstrated).

4 Mathematical Details

Construction of the 6D reconstruction space

To work with high symmetry space-groups and limited values of intra-granular orientation spread, a regular sampling of Rodriguez parametrization of orientation space can be used (see [8] and suggestions therein for a computationally efficient representation applicable to all space groups, based on a compound description based on quaternions and "local" Rodriguez vectors). The voxel size in the real-space volumes will be defined by the acquisition resolution and the volume size can be estimated from the convex hull of the back-projected diffraction spots.

A natural choice of the orientation space resolution would be the angular step size used in the acquisition procedure. However, since an efficient implementation of the iterative reconstruction algorithm requires the full 6D volume to be loaded in computer memory, the sampling interval may have to be decreased to multiples of the acquisition step size in order to avoid saturation of the memory. Note that the orientation space resolution is linked to the total number of orientations and hence memory requirements by the third power of the inverse of the orientation step size. The use of adequate mathematical priors and interpolation techniques discussed in [10] help to mitigate this problem and allow 5 136

for the discretized orientations to represent the signal on the detector coming from the missing orientations. A different approach could be the reduction of the real-space resolution, to compensate for a finer orientation-space sampling grid, which could be used for reconstruction of bigger portions of the sample.

A conservative estimate of the bounding-box of the orientation distribution function of each grain can be obtained from the extreme values of the diffraction blob volumes in real space. The direction of the diffracted beam can be parametrized by two angles: 2θ being the (Bragg) angle between the diffracted beam and the incoming beam, and η being the angle between the projection of the rotation axis and the line connecting the intercepts of the direct beam (u 0 , v 0 ) and the diffracted beam (u, v) on the detector (azimuth angle on the virtual Debye Scherrer ring associated to the current grain position). While 2θ is not supposed to change in this 6-dimensional framework, for a given reflection, the deformation will be observed as a spread in both the ω and η directions. As described in [7], a given set of (θ, η, ω) angles will define a line in the Rodrigues representation of orientation-space. For a given reflection the lines produced by the extremal values of η and ω will not perfectly lie in the same plane, but the deviation will be small enough to use them to approximately define a limiting plane in orientation space for the Orientation Distribution Function of the grain being analysed. By collecting each of the four planes defined by the four extreme deviation values of η and ω for each reflection, it is possible to build a circumscribing polyhedron to the grain ODF.

Convex optimization algorithm

The Chambolle-Pock algorithm (CP, [3]) is a first-order primal-dual algorithm for convex optimization problems, which is among the most popular algorithms for solving total variation minimization problems. The CP algorithms can solve different types of optimization problems, and it is relatively easy to mathematically derive the algorithm 1, tailored for the formulation of the optimization problem in (eq. 3 of the article) [9]. A practical solution to avoid tuning of the free parameter λ in equation ( 3) of the article consists in rewriting the functional as

x * = argmin x || (|∇Sx|) || 1 (1)
subject to:

x ≥ 0 and ||Ax -b|| 2 ≤ (2)
where is now an estimation of the total noise of the projection data. The formulation in equation 1 still contains a free parameter. However, performing a pre-selection of diffraction blobs based on intensity criteria one can reduce the noise level in the projection data. Considering the transition → 0 leads to:

x * = argmin x || (|∇Sx|) || 1 (3)
subject to: x ≥ 0 and Ax = b (4) 6 137

The CP algorithm instance tailored for this new form of the functional is given below in algorithm 2.

The parameter σ 1 is typically chosen as an approximation to the inverse of the l 2 -norm of the matrix A, while the parameter σ 2 is typically chosen as 1/2 which is the inverse of the l 2 -norm of the matrix representation of the gradient.

ALGORITHM 1: CP Instance for eq. 3 of the article Require: x 0 , y, λ Ensure:

Reconstruction in x k p 0 ← 0 q 0 ← 0 θ ← 1 for k ← 1, . . . , n do p k ← p k-1 + σ 1 (Ax k-1 -y) 1 + σ 1 q k ← λ q k-1 + σ 2 ∇Sx k-1 max(λ1, |q k-1 + σ 2 ∇Sx k-1 |) x k ← P 0 x k-1 -τ A T p k -τ S T div q k x k ← x k + θ (x k -x k-1 )
end for ALGORITHM 2: CP Instance for eq. 3

Require: x 0 , y Ensure: Reconstruction in x k p 0 ← 0 q 0 ← 0 θ ← 1 for k ← 1, . . . , n do p k ← p k-1 + σ 1 (Ax k-1 -y) q k ← q k-1 + σ 2 ∇Sx k-1 max(1, |q k-1 + σ 2 ∇Sx k-1 |) x k ← P 0 x k-1 -τ A T p k -τ S T div q k x k ← x k + θ (x k -x k-1 ) end for
If we consider that in our applications, the norm of the matrix A can be large (and therefore the parameter σ 1 is small), we see that 1 + σ 1 ≈ 1, such that Alg. 2 is almost identical to Alg. 1 with the choise λ = 1.

138

The algorithms (1) and ( 2) are made of few important pieces: a projection in both the projection data space and the space associated to the regularization condition in the objective function, the back-projection, and a memory step.

To understand how all these steps in the algorithm can be performed in our 6-dimensional problem (eq. 3 of the article) and (eq. 3), we will now introduce the underlying data-structures. We have two main data-structures that hold the information: the collection of all the used diffraction blobs, also known as diffstack (vector b) and the collection of all the volumes related to the sampled orientations, also known as solution vector (vector x). This means that the diffstack contains all the information recorded by the detector, while the solution vector contains the 6dimensional volume in orientation and real space, introduced in the "Method" section of the article. These two data structures are linked by a third object, generically named geometry (matrix A). Each line of the matrix A contains the contribution of each element of x to one single pixel on the diffracted images, but using the ASTRA toolbox [6] it is possible to use a simpler description, which consists of a collection of tables that fully describe the projection of the 3D volumes onto the 2D detector images.

5 Hardware, software implementation and computational times

The reconstructions were performed on the ESRF computing cluster and distributed over 10 machines with double Tesla K20 GPUs and 128GB of RAM.

The software was implemented in Matlab 1 and C++, using the ASTRA Toolbox (https://github.com/astra-toolbox ) for the projection and back-projection of the volumes. The computational times scale linearly with the number of sampled orientations, size of the real-space volumes and projection data, which means that for the smaller grains, reconstructions can take up to a few minutes, while for the biggest cluster reconstruction twelve hours are needed on one of the above-mentioned machines. The computation times could still be greatly reduced using the newer generations of graphics chipsets, because the biggest cost center in the reconstruction resides in the forward-projection and back-projection of the real-space volumes.

Chapter 5

Twin reconstructions

This chapter concludes the series of chapters dedicated to the inclusion of published articles or articles that are undergoing submission to scientific journals as of the writing of this thesis. This final article is a further extension and refinement of the framework introduced in chapter 3 and 4, again demonstrated on experimental data related to one of the more challenging types of microstructures, encountered in practice: the presence of thin bands of deformation twins.

In some materials, as a result of the applied strain or stress, the crystal structure can deform in sub-regions of some selected grains, to release such stress.

The new structures will typically be in the form of plate-like grains, called twins inside bigger grains called parents. Twins are intrinsically difficult to deal with, because they appear in deformed materials and so the possible presence of elastic strain is not captured in the six-dimensional model, but also because they tend to be very small and this is the source of many effects that make them difficult to index and process. The fact that twins are very thin, and so they only have a small volume compared to the other grains, results in their reflections being usually much less intense and so more sensitive to noise. In the case of the more intense "edge-on" reflections, instead, the point spread function (PSF) of the detector blurs their shapes out.

It is however possible to predict the new orientation of the twins among few possibilities, once the orientation of the parent is known. Indeed the transformation known as twinning acts as a rotation of the crystal along specific symmetry axes of the crystal, by well defined associated angles.

The interesting fact about the six-dimensional reconstruction used in this article for the twins is that it has to deal with two or more disjoint bounding boxes in orientation-space at the same time: the parent's bounding box and the twins' bounding boxes. This makes the orientation-space description in the six-dimensional framework non-contiguous any more. On top of that, other complications arise: (a) to incorporate the reflections of the parent and the twins in one single reconstruction, it is necessary to correctly estimate the hklfamilies' structure factor, to be able to correctly renormalize those reflections to 142 CHAPTER 5. TWIN RECONSTRUCTIONS not create inconsistencies in the data given to the algorithm, (b) the PSF of the detector has to be introduced to account for the smearing of the twin projection shape, but also to prevent the complementary phenomenon, the smearing of the dip caused by the twin in the reflections of the parent, from happening (c) in the case of non-indexed twins, the prediction of the twin variants, to be able to reconstruct them anyway. Once these complications have been taken care of, the reconstruction problem, turns out to be easy and straight forward in the six-dimensional framework developed in this thesis, because it is just enough to sample both regions in orientation-space, compute where they will project on the detector and use the compute projection geometries on the associated reflections.

What is also very interesting is that for some reflections, since the twinning happens along specific symmetry elements of the crystal, both parent and any of the twins will project to the detector at the same position. This can be viewed as a sort of crosstalk between the two bounding boxes in orientation-space. Moreover, if in the reconstruction functional, some term like the TV (•) operator is applied to the sum of the real-space parts of the two six-dimensional bounding boxes of the parent and the twin, which happen to be coincident, there will be another sort of crosstalk between the two grains. This is actually a wanted effect because the twins, as said before, spawn from the parent's volume.

Another very interesting topic mentioned in this article, that we couldn't cover enough because it was not the main message we wanted to pass, is that it can be possible to roughly estimate the scattering intensities of each hkl-family from near-field data. One can think of different approaches to this estimation and depending on the experimental conditions, one could work better than another one. First of all we have to observe that each reflection from a given hkl-family will have a fixed θ that is almost identical for each grain, if strain is negligible. So, Friedel pairs define diffraction spots' θ, and this means that we can identify on the detector all the spots from all the grains that belong to a certain family.

If the spot segmentation is done correctly, and the Debye Scherrer rings associated to a given set of hkl families are complete, it would be enough to compute a relative scattering intensity for each of the hkl-families by using the following formula:

I hkl = 1 N hkl n∈hkl I n | sin (η n ) | (5.1)
where the sum goes over all the spots identified as belonging to the family hkl, N hkl is the degeneracy multiplicity of the given family, I n is the total intensity of each spot, and η n is the associated η, to account for the Lorentz factor. Polarization effects, incoming beam profile and beam attenuation could also be taken into account, but they would just be small corrections and except for the first of them, the others would cancel themselves out during the integration. The obtained intensity I hkl is the scattering intensity of the whole volume for the given family hkl.

143

An alternative way of computing an equivalent estimation of the scattering intensities would be to have first reconstructed enough grains without any assumption on the experiment scattering intensities, by just using 1 as renormalization factor for the spots, and computing the average spot intensity for a given hkl family, renormalized by the Lorentz factor and the total number of voxels in the segmented volume of the grain they belong to. Repeated for all the grains, this would result in the scattering intensity of a single voxel for the given hkl family.

In the attached article, the scattering intensities were computed using the first method and they showed a much improved match with the experimental data, compared to the theoretically computed ones from pure titanium, since the alloy used in the study case was a mixture of titanium with a 4% of aluminium. Moreover, as it was earlier said, the correct estimation of the scattering intensities for the larger θ intercepted by the detector might be underestimated due to the full ring not falling on the detector, but also to some sort of vignetting on the detector, that is not corrected by the earlier stages of DCT. Finally, to further discuss the distortions introduced by the detector and the optics used, at section 2.6 in the attached article, it is possible to find a description of the point spread function phenomenon and its implications. Even if the modelled PSF in the article was quite crude and definitely under-estimated in terms of width and richness of the modelling function, it made however possible to reconstruct thin plate-like features like the twins inside larger structures as the parent grains. The main conclusion from that result is that a better measurement and modelling of the PSF of the used optics and detector will be needed to reach higher resolutions in the DCT reconstructions and so to allow even better accuracy in the determination of the grain boundaries and in the reconstruction of smaller features like sub-micron twins.

A feasibility study of full-field X-ray orientation microscopy at the onset of deformation twinning Nicola Viganò a,b,c , Laura Nervo Abstract 3D X-ray orientation microscopy based on X-ray full field imaging techniques such as diffraction contrast tomography is a challenging task when it comes to materials displaying non-negligible intra-granular orientation spread and/or intricate grain microstructures as a result of plastic deformation and deformation twinning. As shown in this article, the optimization of the experiment conditions and a number of modifications of the data analysis routines enable detection and 3D reconstruction twin lamellae down to micrometer thickness as well as more accurate 3D reconstruction of grains displaying intra-granular orientation spreads of up to a few degrees. The reconstruction of spatially resolved orientation maps becomes possible through the use of a recently introduced six-dimensional reconstruction framework which has been further extended in order to enable simultaneous reconstruction of parent and twin orientations and to account for the finite impulse response of the X-ray imaging detector. The simultaneous reconstruction of disjoint orientations domains requires appropriate scaling of the scattering intensities based on structure and Lorentz factors and yields 3D reconstructions with comparable density values for all the grains. This in turn enables the use of a global, intensity guided assembling procedure and avoids problems related to the single grain thresholding procedure used previously. Last but not least, carrying out a systematic search over the list of known twin variants (forward modeling) for each of the indexed parent grains, it is possible to identify additional twins which have been left undetected at the previous stage of grain indexing based on diffraction spot peak positions. The enhanced procedure has been tested on a 1% deformed specimen made from Ti 4% Al alloy and the result has been cross-validated against a 2D EBSD orientation map acquired on one of the lateral sample surfaces.

Introduction

Four-dimensional (time-lapse) multi-modal observation of polycrystalline materials microstructures as they evolve as a function of strain and/or temperature can provide valuable benchmark data for adjustment and validation of image based, numerical simulations of polycrystalline aggregates. The process of deformation twinning encountered in structural materials with hexagonal crystal structure like Mg, Ti and Zr is an example of industrial relevance (forming of critical components in aerospace and nuclear applications) where the predictive capabilities of theoretical models for microstructure evolution are not yet satisfactory and where refinement and extension of existing models could potentially be facilitated by direct comparison of model predictions with experimental observations. However, time-lapse observations of deformed and/or twinned microstructures by means of three-dimensional synchrotron X-ray diffraction techniques turn out to be a challenging task. Established methods suitable for orientation mapping in deformed microstructures like Differential Aperture X-ray Microscopy (DAXM) [7], X-ray diffraction tomography [1,4] or 3DXRD and HEDM [14,8] use focused X-ray beams and are based on two or three dimensional acquisition procedures. Due to limitations in instrumentation and detector technology these scanning procedures currently suffer from time overheads which limit the ultimate acquisition rate and lead to compromises in terms of sample volume and/or time steps which can be recorded in the course of a typical synchrotron beamtime. Tomographic, full-field microstructure mapping techniques like Xray diffraction contrast tomography (DCT) [9,15] require only a single rotational scan. Given the longer exposure time for individual images (as compared to 1D or 2D focused beams) they are less affected by above-mentioned technical limitations. As a consequence, relatively large (10 6 -10 7 voxel) sample volumes can be characterized in a fraction of the time currently reported for 2D and 3D scanning approaches1 However, until recently the reconstruction approach behind DCT has been lacking the ability to address materials with non-negligible intra-granular orientation gradients. With the introduction of a six-dimensional reconstruction framework [20,21] this technique can now be extended to the case of moderately deformed materials displaying intra-granular orientation spreads of up to a few degrees. In this paper we present an extension of this full-field X-ray orientation imaging approach tailored for the specific challenges related to reconstruction of materials having undergone moderate levels (up to a few %) of plastic deformation, partly accommodated by the formation of deformation twins.

Challenges related to the reconstruction of deformation twins

Deformation twins are typically observed as disk-like crystal domains with thickness and size distributions spanning from nanometer to micrometer length scales [3]. Their crystallographic orientation is related to the parent orientation by a well defined pair of rotation axis and rotation angles around one of the symmetry axes of the crystal (see appendix A for a brief recall of the calculation of twin variant orientations given the orientation of the parent grain and the twin axis and angle pair of the twinning mode under consideration). This means that all Bragg reflections associated to scattering vectors co-linear to the said twin rotation axes will be shared between parent and twin and project to the same position and diffraction spot on the detector. As a result, while most of the reflections recorded on the detector will be independent projections of the two domains, a subset of them will be shared. If these shared reflections shall be included in a joint reconstruction procedure, where parent and twin domains are reconstructed simultaneously (as envisaged in this work), the diffracted intensities have to be scaled properly based on structure and Lorentz factors 2 . The high aspect ratio and small scattering volume of deformation twins result in diffraction signals which can be orders of magnitude weaker than the one of the corresponding parent grain. Moreover, due to the high aspect ratio, the local intensity received on the high resolution screen varies as a function of projection direction and systematically drops to the noise level for projections close to the direction of the twin plate normal. These systematic absences are in turn a great challenge for the indexation and the tomographic reconstruction algorithm, suffering from the poor angular coverage. In addition, particularly projections in directions perpendicular to the plate normal ("edge-on" views, observed as intense lines) suffer from blurring due to the limited impulse response of the high resolution imaging detectors employed in this type of near-field diffraction imaging experiments. This together with the partial-volume effect (drop of reconstructed intensity for objects with thickness below the pixel size of the reconstruction) results generally in substantially broadened twin reconstructions and lower reconstructed intensities compared to the parent grains. Last but not least, deformation twins are triggered by and accompanied with other types of plastic deformation (dislocation glide), giving rise to increasing levels of disorder and intragranular orientation spread with increasing deformation. Even moderate plastic deformations in the range of one percent can already result in the formation of grain-substructures and orientation spreads of order of one degree. The resulting three-dimensional diffraction blob volumes are deformed grain projections which violate the parallel projection assumption of conventional (single orientation) DCT, which can no longer produce accurate grain shape reconstructions in this case.

The two obvious experimental solutions to alleviate some of the above-mentioned problems are improving the signal over noise ratio by reducing the background on the detector and increasing the illumination on the sample, as described in section 3.1.

In terms of data analysis strategy the approach pursued in this work was to create the infrastructure for simultaneous reconstruction of the parent and twin, using an upgraded version of the recently introduced six-dimensional reconstruction framework [20,21], as described in the next section. Furthermore, we show that a basic search over all possible twin variants (derived from the known parent grain orientations in the sample volume) can help to reveal additional twins, left undetected by the indexing procedure.

Method

The usual setup of a near-field diffraction experiment, when using a monochromatic beam, consists of a rotation stage where the sample is positioned, and as the sample rotates by angle ω, it gives rise to diffracted beams each time the Bragg condition is met for one of the grains. Some of those diffracted beams will intersect the high resolution imaging detector positioned a few millimeters downstream of the sample, and give rise to diffraction spots, which, in the absence of intragranular orientation spread, correspond to 2D projections of the 3D grain volume. After diffraction spot segmentation and indexation based on Friedel pairs (hkl and hkl reflections of the same grain observed for ω 0 and ω 0 + 180 • ) the 3D grain structure can be reconstructed by means of iterative tomographic reconstruction techniques. The reader interested in details concerning the setup, acquisition procedures and initial processing steps like segmentation, Friedel pair matching and indexation of near-field X-ray diffraction data is referred to [9] and [15].

In the presence of non-negligible intra-granular orientation spread, the parallel projection assumption used in conventional (three-dimensional) DCT gets violated. Each of the sub-orientations present in a grain is associated to a slightly different projection geometry and the diffraction signal associated to a given Bragg reflection is observed as a distorted, three-dimensional diffraction volume, which then takes the name of diffraction "blob". It is parametrized by two spatial coordinates u and v (detector pixel coordinates) and a rotation angle ω (image number) 3 .

It has been recently shown that this type of diffraction data can be inverted by switching to a six-dimensional description, in which case maps of local (intragranular) orientation can be produced using suitable optimization and regularization techniques [20,21]. Unlike the case of 3D-DCT, where only one single (grain average) orientation is considered per grain, the six dimensional formulation assigns a discrete 3D orientation distribution to each of the real space voxels inside a grain. In order to compare this 6D scalar field reconstruction output to conventional 3D vector field representations of the orientation field, we transform the 6D output into a 3D orientation map, by assigning the average orientation (center of mass of the local orientation sub-space) to each real space voxel. Since this six-dimensional framework intrinsically accounts for the distortion of the diffraction blob volumes it results in improved shape reconstructions compared to the single orientation (3D) reconstruction framework.

In order to adapt the 6D framework to the case of deformation twins we further extend the previously introduced concept of "cluster" reconstructions [21] to the case of disjoint orientation domains (applicable to other cases of grain sub-structures with a priori known orientation relationships, or grain neighborhoods), thereby enabling simultaneous reconstruction of parent and twin domains from experimental raw data. As mentioned earlier, for this to become possible the previous concept of diffraction spot intensity re-normalization has been abandoned and replaced by an appropriate scaling of projection matrix, taking intensity variations between different hkl reflections due to variations in structure factor and Lorentz factors explicitly into account. Apart from alleviating problems with segmentation of parent and twin reconstructions, the simultaneous reconstruction has the added benefit that additional (shared) reflections between parent and twin domain can be included in the reconstruction process. Moreover our model of diffraction image formation has been enhanced and now includes blurring effects caused by the finite impulse response of the X-ray imaging detectors employed in this type of experiments.

Mathematical formulation

As in [20], the reconstruction space, if elastic strain is supposed to be negligible, is a 6-dimensional space X 6 = R 3 ⊗ O 3 obtained by the outer product of cartesian real-space and 3-dimensional orientation space. The representation of this 6-dimensional space is discretized by choosing a sampling that can be viewed in two equivalent ways. The first and most intuitive sampling would consist of creating 3-dimensional orientation sampling for each of the real-space sampled positions. This would give rise to a collection of Orientation Distributions Functions for each of the real-space voxels. Although this representation may be intuitive, it is not very computationally efficient, and we preferred to use a collection of identical real-space volumes, each of which associated to a different sampled point in the orientation-space.

beam and the diffracted beam is called 2θ, and the angle between the projection of the rotation axis and the line connecting the intercepts of the direct beam (u0, v0) and the diffracted beam (u, v) on the detector is called η.)

We would also like to point out that while the parent grain and its twins are reconstructed within the same bounding box in real-space, their bounding boxes are instead separate in orientation-space, and that, within limits of memory occupation, it is possible to include several sub-volumes (e.g. the different twin variants) in orientation-space.

If now we assume that the image formation can be described by kinematic diffraction and that we can neglect photoelectric absorption and extinction effects, the forward-projection of the grain volumes can be represented as a linear operator:

Ax = b (1)
where x is a vector containing N P elements, representing the scalar "scattering power" for each of the sampled volume elements in the six-dimensional position -orientation space (N P = n 3 × kp 3 for the case of a regular sampling over a cube-shaped sub-volume with grid length n in real space and k disjoint subvolumes in orientation space, each sampled with grid length p). Each line of the projection matrix A S×N P contains the contribution of the 6D volume elements to a given detector pixel and the vector b S holds a list of measured pixel intensities, specified by their (u, v, ω) coordinates in the 3D diffraction image stack. S corresponds to the total number of detector pixels in the M = kM hklc diffraction blob volumes corresponding to the joint set of parent and twin reflections recorded during the scan. Here M hkl represents the total number of reflections intercepting the detector, which is usually of order of a few tens up to one hundred, and c represents the small number of reflections shared between the parent and twin orientations. If the elements of x are arranged as a succession of 3D real space volumes, each representing one of the sampled orientations, and A as an array of M × P projection matrices, each one describing the projection geometry for one of the M reflections intercepted by the detector, the vector b will be composed of M blocks corresponding to contiguous 3D subvolumes (i.e. the aforementioned diffraction blobs), spread throughout the entire stack of detector images:

   A 11 • • • A 1P . . . . . . . . . A M 1 • • • A M P       x 1 . . . x P    =    b 1 . . . b M    (2)
The beauty of equation ( 2) is however limited by the fact that it is not suited for the actual reconstruction. The number of diffraction blobs that fall on the detector (M 30 -100) is limited and these spots are affected by noise and overlaps, so equation ( 1) is a heavily underdetermined and perturbed equation system which needs to be re-written as a minimization problem. Some additional constraint and prior should be added to the minimization functional in order to promote solutions which are compliant with the known physical properties of the microstructure under investigation. In the case of a moderately deformed metal sample, imaged at 1.5 µm resolution, as presented in this study, two physically justified prior assumptions are positiveness and sparsity of the 6D solution vector. Indeed, we only expect a limited number of "active" orientations inside a given real space volume element and this behaviour could be enforced by adding a l 1 penalty term to the optimization functional, as proposed in [20]. Here we have used an alternative formulation, minimizing the spatial variability of the reconstructed intensity4 as proposed recently in [21]:

x * = argmin x ||Ax -b|| 2 + λ|| (|∇Sx|) || 1 (3a) subject to: x ≥ 0 ( 3b
)
where S is the operator that sums all the orientation components from both the parent and the twins orientation-spaces for each real-space voxel, and the l 1 -norm over the absolute value of the gradient is the total variation operator [2,16].

In this work however, while we formally keep the same functional because the general idea is still completely valid, we had to adapt many details under the cover of equation ( 3), to be able to successfully deal with the intrinsic nature of the twin reconstructions.

From physics we know that the parent and its twins occupy exclusive but also contiguous regions in the real-space reconstruction domain. This means that they do not overlap but also that there is no empty spacing among the two domains. To force an homogeneous real-space volume with sharp boundaries at the grain edges and transitions at the parent-twin boundaries that are also homogeneous, the most natural choice is a total variation (TV) minimization over the sum of all the orientation components of the real-space voxels, which means that this time, the sum is over two disjoint regions in orientation-space.

As in [20] and [21], the algorithm that solves this optimization problem was derived from a specific class of algorithms called Chambolle-Pock, slightly modified to cope with the new type of functional, as described in [16].

In the next sections we will also analyse other peculiar details of equation ( 3), that mark the progress achieved compared to the work in [21].

Scaling of the projection matrix

From the mathematical point of view, the renormalization comes from the decomposition of equation ( 1), as the following:

Ax = C Ãx = C N i Ãi x = b (4)
where the matrix C is a diagonal matrix, which predicts the scattering intensities for each of the reflections, and the matrix à only contains the geometric projection coefficients to the detector. The re-normalization procedure simply pre-multiplies both sides of the equation 4 like:

Ãx = N i Ãi x = C -1 b = b (5)
where now b are the re-normalized diffraction blobs. The advantage of this new formulation is a better conditioning of the problem, being the diagonal of the matrix à more homogeneous. Finally, while most of the reflections are not shared between a parent and its twins, and so the projection to the diffraction blobs happens separately, for what concerns the shared reflections, the forward-projection and back-projection should simultaneously happen for both the parent and the twins to the same diffraction blobs.

Point spread function and blurred reconstructions

Near-field diffraction imaging experiments are typically performed with detector systems based on CCD cameras, coupled via visible light optics to a transparent luminiscent screen that converts the incoming X-Rays into visible light. Due to the large interaction volume of X-rays and the diffraction limit of the visible light microscope optics common designs of these systems are known to be limited in spatial resolution to values of typically a few times the effective detector pixel size. The performance of such an optical system can be described by its impulse response (point spread function, PSF) and/or the modulation transfer function. For X-ray imaging systems, both depend on details of the experiment conditions (i.e. energy of X-ray beam, type and thickness of scintillator crystal, visible light optics, ...). In the present study the point spread function was approximated by a symmetric, translational invariant Gaussian profile with a full width of half maximum of about 2 pixels. (As further discussed in section 5 a more realistic description of the PSF would depend on the detector region and on the incidence direction of the photons on the scintillator screen.) If we now take equation (5), where the scattering intensity has been translated into a renormalization of the diffraction blobs, and so the forward model is expressed in a pure geometrical form, it can be modified to include the PSF blurring like the following:

F Ãx = F N i Ãi x = b (6)
where the matrix F is the operator responsible for applying the PSF to the forward-projected images. Practically speaking this translates in a simple convolution between the forward projected diffraction blobs and the PSF specific to each of them. When trying to invert the forward operator F Ã, by first pre-multiplying it by its transpose, we have to take into account that the PSF is generally speaking not symmetric, and so the adjoint of F should be taken, resulting in:

ÃT F † F Ãx = ÃT F † b (7)
where the superscript † stands for adjoint (complex conjugate transpose).

Finally, these considerations about the PSF were always to be considered negligible in the previous works using the 6-dimensional approach, being the lack of orientations-space resolution the major limiting factor in the reconstruction quality. However, in the case of twin reconstructions, where the twins can be sometimes less than 1µm thick, and the detector pixel resolution is above 1 µm, even moderate effects like the PSF blurring of very few pixels, could heavily impact the quality of reconstruction of these thin regions of the grains.

Finding twins from the parent orientation

For the reasons explained in section 1.1, deformation twins are intrinsically difficult to index from near-field diffraction data using conventional indexing schemes based on diffraction spot peak positions. Nevertheless, the detection ratio can be improved by a systematic search for intensity (forward modeling) in the background corrected raw images, starting from the known orientation of the successfully indexed parent grains. In the following we restrict our discussion on one specific twinning mode, predominant for the material studied in this work, but the suggested strategy can be applied to any other known twinning modes.

The Ti-4 wt.%Al alloy studied in this work has a hexagonal crystal structure with a c a ratio smaller than the one for optimum close packing. When deformed in compression in the direction parallel to the initial rolling direction, the most frequently observed twinning mode for this type of material is one known as tensile twin of type 1 (TT1) [3]). This twinning mode corresponds to a rotation of about 85 • around [10 10] and the hexagonal crystal symmetry allows for six different variants of the same type. Given the parent orientation, the orientation matrix of the twin variants can be precomputed (see appendix A) and from this the approximate twin diffraction spot peak position ( (u, v, ω) coordinates in the detector-space) can be determined. Since twins are confined within the volume of the parent grain we can restrict the search of the twin reflections to diffraction blobs that do not exceed the size of the parent grain. At the same time we allow for blobs that have a high variability in aspect ratios, as expected from the disk-like morphology of deformation twins. However, the actual twin orientation might be slightly different from the theoretically predicted orientation (derived from the average orientation of the parent grain), and the position of the twin inside the parent grain is not known a priori. So it is convenient to express the detector region for the search of the twin reflections in terms of the (θ, η, ω) coordinates, allowing for small deviations from the nominal values, in order to account for these uncertainties.

Experiment and material

A sample of 2.5 mm length and a diameter of 0.5 mm was extracted by electrodischarge machining parallel to the cylinder axis and rolling direction of a larger compression sample made from binary Ti-4 wt.% Al alloy. The larger specimen had been deformed in compression to about 1 % plastic strain in the direction parallel to the cylinder axis, i.e. parallel to the former rolling direction. The preferential orientation of the < c > axis in the direction perpendicular to the cylinder axis and loading direction promotes tensile twinning during compression loading as the < c > axis is strained in tension [5]. The grain size of the material determined by the linear intercept method was 70 µm. In a parallel study on the same specimen an analysis of twin nucleation criteria based on a 3D reconstructed sample volume containing 400 grains and about 60 twins revealed a correlation of the propensity for twinning with grain size and a clustering of twinned grains in neighborhoods with favourable prismatic slip transmission conditions [12].

DCT data acquisition

The experiment was performed on the 3D XRD instrument of the materials science beamline ID11 of the European Synchrotron Radiation Facility (ESRF). The incident beam has been monochromatized to 60 keV and focused to dimensions of approximately 0.2 mm times 0.5 mm by means of a set of compound refractive lenses (CRL). Under experiment conditions described above, the focusing results in about a factor of 20 increase in flux density at the sample position compared to the configuration without lenses. This increase turns out to be a crucial requirement for successful detection of deformation twins which are know to nucleate as disk-like tapered lamellae of (initially sub-) micrometer thickness. The cylindrical sample was mounted with its cylinder axis parallel to the rotation axis of the diffractometer and the focused incident beam has been confined by slits to a section of 0.15 mm in the vertical and 0.6 mm in the horizontal direction. The sample to detector distance was set to 8 mm and the diffraction images were recorded on a CCD based and light optic based high resolution detector system featuring a 50 µm thick, transparent luminiscent screen prepared from a Eu doped Gallium Gadolinium Garnet single crystal (GGG) [11] and a long-working distance visible light microscope objective (10×, NA= 0.23). The resulting effective pixel size was 1.4 µm. As mentioned earlier, increasing the flux density by at least one order of magnitude compared to acquisitions in materials with comparable grain size turns out to be a prerequisite for successful indexation and reconstruction of deformation twins down to micrometer thickness.

The introduction of the absorber significantly reduces the background signal on the scintillation screen and allows to raise the exposure time and hence diffraction spot intensities by the inverse of its X-ray transmission (≈ 20 times in the case of 0.5 mm of Ta at 60 keV).

For the DCT scan a total of 3600 images with 2.5 s exposure time were recorded during a synchronized, continuous 360 • motion of the air-bearing diffractometer stage. A second scan with 0.2 s exposure time and 360 projections but without the additional absorber was recorded right after the DCT scan, maintaining the relative position of sample and detector. These data were used for tomographic reconstruction of the X-ray attenuation coefficient, serving as a mask in the process of assembling the individual grain volumes into the common sample volume.

EBSD analysis

The EBSD characterization was performed on a cross-section parallel to the axis and loading direction of the cylindrical sample. After embedding the sample into epoxy resin a layer of approximately 200µm was removed using 4000 grit SiC grinding paper. For best Kikuchi pattern acquisition, the sample was polished with Oxide polishing suspension (80% colloidal silica + 20% H 2 O 2 by volume) for ≈ 2 h and etched repeatedly with Krolls reagent. The surface finish was further improved by a final vibratory polishing with Oxide polishing suspension in a Buehler VibroMet 2 machine. To overlap the 3D region scanned in DCT, a 600 x 90 µm 2 area was selected and scanned using a step size of 0.4 µm in a field emission gun scanning electron microscope (FEI Quanta 650) equipped with a NordLYS EBSD detector. The EBSD mapping was performed at an operating voltage of 20 kV with a working distance of 13 mm. The AZTEC software was used for EBSD data acquisition. Post analysis of the EBSD data like inverse pole figure map generation, axis angle misorientation pair for determination of active twin mode was performed using HKL Channel 5 software.

Results

From inspection of the EBSD map figure (1) acquired with a step size of 0.4µm, it appears that most of the twins intersecting the sample surface correspond to tensile twins (85 • rotation around [10 10])). In the common surface area shared between the EBSD and DCT acquisition two grains have developed (sub-)micrometer thickness twins of this type. One of them (fig. 2,a) has been successfully identified during the indexing step of the DCT data reduction process [9,15] whereas two other domains could only be detected during a systematic search, using the prediction (forward simulation) of twin orientation variants calculated from the orientation of the parent grain (fig. 2,b), as it can be seen in appendix A. In figure (3) shows a slice though the reconstructed 3D intensity volume (i.e. the sum of all orientation space components of the 6D reconstruction output) of the grain rendered in figure (2, a), illustrating the difficulty of correctly reconstructing the precise shape and volume of these sub-micrometer thickness twins domains. In fact, while the thickness of twin lamellas should be smaller than a single voxel, in the reconstruction they extend over several voxels, and while they should present a reconstruction intensity that is comparable In addition, an even thinner twinned region is visible in figure (5, a). Its thickness appears to be smaller than the acquisition step size of the EBSD map (0.4µm), and it was not possible to identify this second twin with the DCT acquisition conditions described in this paper.

In figure (5) we take a closer look to the smooth deformation gradients in one of the grains of the analysed surface region, and compare the results between the EBSD map and the 6D-DCT reconstruction. The analysed grain was unfortunately damaged by few scratches during the polishing phase, but it is also the one that shows the strongest and most easily recognizable orientation gradient among the reconstructed surface grains. Indeed, if we exclude the obvious deformation due to the scratches in figure (5, a), and remember that the same limitations introduced before still apply, a similar deformation gradient can be observed in both figures (5, a) and (5, b). What is not completely obvious from figure (5), is the presence of slip bands inside the grain. This is instead noticeable both from the band contrast map obtained for the EBSD measurements ment, the local misorientation in vicinity of these slip bands were too small to be clearly distinguished. Finally, in figure (7), we can see a comparison of the assigned average orientation per each real-space voxel, and the total scattering intensity assigned by the 6-dimensional reconstruction algorithm. The similar average intensities among all the grains (the variations are around few percents) are due to the renormalization of the diffraction blobs discussed in the previous section 2.2, where instead of renormalizing all the spots belonging to the same grain to a common, arbitrary value, we decided to re-normalize them based on theoretically predicted scattering intensities per each crystallographic family and different Lorentz factors. The main outcome of this procedure is that these intensities are now more physically meaningful and can be compared among different grains, in order to solve grain conflicts. This in turn results in better accuracy of the grain boundary and triple junction positions, and yields grain maps which show an higher coverage of the sample volume, even before any type of dilation is applied. Some remaining deviations in the average intensity can be observed for some of the grains in figure (7). The physical origin of these variations remains to be clarified.

Discussion

The high aspect ratio of deformation twins gives rise to large (up to two orders of magnitude and higher) variations of the signal received on individual pixels of the high resolution imaging detector as a function of rotation angle. In order to raise the signal of unfavourable reflections (i.e. those with direction of the diffracted beam close to the direction of the surface normal of a twin lamella) above the background, and hence assure uniform angular coverage of the projections from a twin domain, it is mandatory to increase the photon flux density at the sample position far beyond the level usually employed for an equiaxed grain structure. With the flux provided by state of the art undulator insertion devices and the use of tunable X-ray focusing optics [19] the total scanning time for the 150 × 500 µm sample sub-volume was 2.5 hours. In order to extend the detection limit towards even smaller (thinner) twins, further increase of the photon flux density is required and will become available with the one to two order of magnitude increase in brilliance predicted for the next generation of low emittance storage rings. Alternatively, one may consider switching to 1D focused line beam illumination mode, compromising ultimate time resolution but enabling alternative reconstruction strategies (forward modeling [18,8] or global 5D optimization approaches for simultaneous reconstruction of entire sample cross section) which may outperform the current 6D framework based on grain by grain reconstruction at higher levels of deformation.

In the present work we have incorporated a rather crude estimate of the point spread function modeled as a rotational symmetric Gaussian profile. Further improvement can be expected from a measurement of the actual profile for different angles of incidence on the detector and introduction of an orientation dependent PSF, reproducing the directional asymmetry caused by the parallax of the diffracted beam impinging with different angles (2θ Bragg , η) on the scintillator screen with finite thickness. The oversimplified model of the detector impulse response together with the partial volume effect known from conventional absorption tomography are supposed to be the principal reasons for the lower reconstructed densities of voxels corresponding to twinned regions. Figure 4 b shows some discrepancies in the grain shapes as determined by EBSD and our full-field reconstruction approach. These discrepancies are partly related to the fact that the X-ray orientation map has been extracted from a virtual surface situated approximately two pixels (3 µm) below the actual sample surface, as determined from the absorption volume.

As mentionned earlier, we notice some gradient in the reconstructed intensity values over the first few pixels of a grain volume. The presence of this gradient renders gray scale segmentation of these volumes problematic and the small offset of the surface from which the orientation values have been extracted from helps to overcome this segmentation problem but starts to introduce some visible shifts, especially for boundaries intersecting the sample surface at a small angle. Moreover, the sample surface is slightly curved (artefact from mechanical polishing procedure). We have determined the actual shape of the surface from the phase contrast reconstruction of the sample and we do have extracted the orientations (figure 4, b) from the intersection of this curved surface with the (2 pixel shifted) 3D orientation volume. Within the limits of computer memory available, the formalism used in this paper allows for simultaneous reconstruction of arbitrary sets of disjoint orientations occupying a common volume in real-space. By extending the real-space volume to a cluster of neighboring grains, it is therefore possible to use the same framework for simultaneous reconstruction of 3D grain neighborhoods. This procedure would circumvent the step of grain shape determination by means of gray-level segmentation of the intensity volumes of the individual grains, and it can be expected that such combined reconstructions will result in an overall improvement of the 3D grain shape reconstruction accuracy, similar to the one observed for deformation twins presented in this study.

The current reconstruction approach is still essentially a two-step process of grain indexing based on scattering vectors determined from diffraction spot peak positions. Due to increasing levels of spot overlap, the concepts behind grain indexing are known to fail gradually for increasing levels of plastic deformation.

On the other hand, by repeated characterization of a sample volume at increasing levels of plastic deformation it should be possible to seed the reconstruction of the next step with the information from orientation maps reconstructed at previous step(s), thereby extending the range of plastic deformation which can be addressed with the current approach beyond the level of the feasibility study presented in this work. It can be expected that additional diffraction data acquired on a far-field detector with a pixel size comparable or larger than the grain size and higher quantum efficiency would allow for for identification of twins at early stages during deformation since the entire twin volume would typically project on individual pixels, thereby removing the systematic intensity variations as a function of the angle between the twin lamella and the diffracted beam direction, observed on the near-field detectors. Using the center of mass position and orientation determined from far-field data, it would be possible to use these twin "seeds" in the joint 6D reconstruction of the corresponding parent grains. The example shown in figure 6 confirms the visibility of topographic orientation contrast related to the presence of slip bands in grains which are favourably oriented for prismatic slip. Unfortunately, the sensitivity of the 6D-DCT reconstruction is not high enough to reveal the faint lattice rotations giving rise to topographic orientation contrast observed in diffraction spots with favourably oriented scattering vectors. We anticipate that by using a combined characterization approach with DCT as the initial characterization step and topotomography [10], or dark-field microscopy [17] as subsequent "zoom" on individual grains and grain neighborhoods, it will be become possible to resolve the 3D configuration of these structures. This in turn will allow for studying slip transfer and slip -twin interactions as the one observed at the intersection at the grain boundary between grains A and C (figure 1 a and b) in the bulk of the material.

Last but not least it should be mentioned that part of the remaining inaccuracies might be due to non-negligible elastic distortions of the crystal lattices in the vicinity of twins. The six-dimensional framework used in this work does not account for these additional degrees of freedom and additional work is required in order to estimate in how far this and other previously mentioned simplifications of our model for diffraction image formation affect the accuracy of the orientation maps.

Conclusions

The previously introduced six-dimensional framework for reconstruction of spatially resolved grain orientation maps from extended beam, near-field diffraction data has been extended and adapted to the characterization of structural materials with hexagonal crystal structure displaying intra-granular orientation spreads of up to a few degrees and micrometer thickness twin lamellae as a result of plastic deformation. The optimization of acquisition conditions and an explicit search over possible twin variants enables detection of deformation twins down to sub-micrometer initial thickness. Introducing an appropriate scaling of the diffracted intensities by structure and Lorentz factor, accounting for the finite impulse response of the high resolution detector and a global optimization carried out over a 6D volume covering the parent and twin orientation(s) result in a clear improvement of the accuracy of the orientation maps when compared to previously used single orientation (3D) reconstruction framework of X-ray diffraction contrast tomography. Cross-validation with EBSD data acquired on the sample surface show good overall agreement but also some remaining degree of blurring and thickening of the twin domains in the X-ray orientation maps.

The non-destructive nature, straightforward combination with other tomographic imaging modalities (phase contrast, topotomography, dark-field microscopy) and improved time resolution compared to other 3D orientation mapping techniques makes this approach a promising candidate for time-lapse observations at initial stages of plastic deformation in structural materials.

where the vector r (a i , α) = a i tan 1 2 α.

B Alignment of EBSD and DCT reference systems

The DCT reconstruction and the EBSD map, naturally come with different reference systems, and so, their orientations are not directly comparable. Indeed, in the orientation-space, is the underlying crystal symmetry is not not considered, the relationship that connects the same orientation expressed in two different reference systems (namely, the DCT and the EBSD reference systems) can be represented like:

g DCT = g EBSD g DCT →EBSD ( 10
)
where g DCT is the orientation matrix of the chosen orientation in the DCT reference system, while g EBSD is the orientation matrix of the chosen orientation in the EBSD reference system, and g DCT →EBSD is the orientation matrix of the EBSD reference system expressed in the DCT reference system. However, if the underlying crystal symmetry is taken into account, the following equation is then valid:

F (g DCT ) = F (g EBSD g DCT →EBSD ) (11)
where F : O 3 ⊆ R 3 → O 3 ⊆ R 3 is the mapping of the orientations to the fundamental zone of the crystal orientation space. Unfortunately, while equation ( 10) is always invertible, equation (11) is not, unless the magnitude of the orientation-space Rodrigues vector associated with g DCT →EBSD is small enough to make equation (11) become equivalent with equation (10). This means that for small calibrations of two reference systems, it is enough to invert equation (10) like the following:

g -1 EBSD g DCT = g DCT →EBSD (12)
But for bigger misorientations among the two coordinate systems, it is important to know the sample geometry and how it is represented in those reference systems, so that an approximate relationship among the two can be found.

Chapter 6

Conclusions and perspectives

The structure of this thesis, except for the introduction, intentionally introduces the work performed during the author's Ph.D in temporal order, in the form of articles. By following the evolution of the proposed method, from a simple description of a discrete representation for the six-dimensional space, in chapter 2 [START_REF] Viganò | Discrete Representation of Local Orientation in Grains using Diffraction Contrast Tomography[END_REF], to the more mature framework used to model the deformation in polycrystals, in different scenarios like the ones of chapter 4 and 5, it is possible to understand how some of the limitations imposed by traditional DCT have been tackled and solved, to address more interesting scientific cases.

The algorithm and the associated six-dimensional representation were initially designed for the case of smoothly deformed single grain reconstructions, while it was not completely clear from tests performed in [START_REF] Viganò | Reconstruction of local orientation in grains using a discrete representation of orientation space[END_REF] (the article included in chapter 3), whether this work could scale to bigger and more challenging cases like the ones in the following chapters. An interesting fact is that between the article included in chapter 3, and the analyses in the following chapters, the available computation power to the end users grew to the point that what was considered a limit in terms of memory in [START_REF] Viganò | Reconstruction of local orientation in grains using a discrete representation of orientation space[END_REF], became quickly addressable soon after. The further development of the combined model and algorithms will be fundamental for the success of the proposed framework, but the continuous and fast paced improvement of the current computational technologies, and the use of hybrid CPU-GPU computations, will contribute to make even more challenging experimental conditions possible in the future.

It is indeed possible that once larger amounts of memory will be available on common workstations (512GB or more), larger portions of the analysed sample volumes could be reconstructed in once instance, following the principles introduced in this thesis. On the other hand, having a fine sampling of both the real-space and the complete orientation-space, which would make it possible to completely discard the indexing phase, is extremely demanding in terms of CHAPTER 6. CONCLUSIONS AND PERSPECTIVES memory and computational power, so it will not be likely to happen anywhere in the near future.

A more affordable intermediate solution could be to perform a reconstruction restricted to the five-dimensional space, where the real-space is limited by a one-dimensional beam illumination to only be two-dimensional. This would slow down the acquisition to the same speed as other approaches like the one discussed in section 1.4.1, but on the other hand it would have the advantage to still perform a global optimization, as opposed to the local voxel optimization of the model in equation 1.49, and possibly, thanks to the use of the just mentioned hybrid CPU-GPU types of computations, it could avoid the use of super-computers or large computational clusters.

Another possible solution to the same problem could be to not restrict the space dimensionality, but to use a multi-resolution approach instead, where the reconstruction is first performed on a very coarse real-space grid, that gets successively refined, as the orientation-space gets instead restricted to the relevant regions.

A third and probably easier route is based on the fact that some commonly used materials exhibit a heavily textured structure, that tends to cluster the local orientations of large portions of the samples in reduced regions of the orientation-space, making it possible to scale the principles of chapter 4 to the size of the said regions, without hitting the mentioned computational barrier. Moreover, heavily textured samples are known to break the principles assumed by traditional DCT, because the diffraction blobs tend to be extremely large, due to bigger regions of the samples having similar orientations that heavily overlap on the detector. The properties of the said blobs, become then extremely unreliable, and this is known to break the indexing algorithm used in DCT. Skipping indexing entirely, and letting the six-dimensional reconstruction algorithm reconstruct those regions could be a viable solution, but it might not just work out of the box, because if no sort of indexing is possible, it is not possible to estimate any of the needed parameters to sample the orientation-space. An interesting solution to this problem could come from the extension of DCT, which traditionally relies on a single-scan near-field acquisition, to the use of scanning procedures with simultaneous acquisition on both, a near-field and a far-field detector (similarly to the concept of a three-dimensional detector proposed by H.F. Poulsen and co-workers). Not only, this would allow other poly-crystal indexing techniques, like "GrainSpotter" [START_REF] Sø Ren Schmidt | GrainSpotter : a fast and robust polycrystalline indexing algorithm[END_REF] (discussed in section 1.4.2), to become compatible with the six-dimensional framework here introduced, but also to enormously extend the capabilities of DCT itself. Far-field data indeed cannot report highly accurate spatial information about the sample, but this means that having a bigger focus on just the orientationspace morphology and distribution analysis, it can be of great help if incorporated in the DCT technique. In the simplest scenario, this would already allow to gather information about the ODF of individual grains inside the sample, and make better decisions when sampling the orientation-space for the single grain reconstructions. But especially for highly textured samples, it is possible to determine, to a certain degree, the ODF of a certain spatially confined region, and then use the current framework exposed in chapter 4, together with the six-dimensional space bounding box determination from the far-field data and reconstruct the said region from the near-field raw images.

Another extremely valuable use of the far-field data would be to use it as an additional constraint for the six-dimensional reconstructions. If we can model a forward projection operator for the far-field problem, which is geometrically very similar to the near-field operator, and call it W , we could think of modifying any of the tomographic minimization problems to incorporate such constraint, like:

x * = argmin where y are the far-field images.

As mentioned earlier, the far-field data provides more accurate information in the orientation-space, while the near-field data provides more accurate spatial information, so the two combined approaches could significantly improve the six-dimensional reconstructions, especially with moderately to highly deformed materials.

As important as adding far-field information to the DCT reconstructions, could be to extend the acquisition geometry with other near-field measurements for different distances of the detector or even positions that don't intercept the direct beam any more. In fact, the acquisition geometry has a rather limited field of view in real-space. This can be easily understood by looking at the acquisition geometry represented in both figures 1.6 and 1.7, where the measurements are, in the sample coordinate system, distributed on a ring around the sample, defined by the rotation of the detector surface around the sample. From this picture it is immediately obvious that directions more parallel to the rotation axis are not sampled in this type of geometry. Unfortunately, even in the comparison carried out in [START_REF] Reischig | Advances in X-ray diffraction contrast tomography: flexibility in the setup geometry and application to multiphase materials[END_REF] between the "inline" geometry, and the geometry with the detector at π/2¸the detector still lied in the said ring around the sample, and thus, even if it gave access to other hkl-families that are not usually intercepted by the inline detector, it didn't bringing more spatial information to the reconstruction in extended-beam mode. An example of alternative approach could be to position an extra detector outside of such ring, possibly at π/2 with the direct beam direction and at minimum π/6 from the sample rotation axis direction, to reduce the ω-broadening given by the Lorentz factor. Such geometry would not however be very easy to optimize, due to the horizontal beam polarization in synchrotrons, forcing some 172 CHAPTER 6. CONCLUSIONS AND PERSPECTIVES additional constraints in the choice of the detector position. Moreover, the absence of the direct beam in the extra detector, would make the calibration and fine tuning of its tilts and position very difficult. A similar kind of geometry is used in topotomography experiments [22] (recently extended by the technique known as darkfield microscopy [START_REF] Simons | Dark-field X-ray microscopy for multiscale structural characterization[END_REF]), where for a specific plain normal of a chosen grain, the diffractometer stages are tilted and optimized to always have such plane normal in diffraction condition across the whole 2π range of ωs. In topotomography, jointly with the rotation around the sample rotation axis, acquisitions for different base tilts are performed for each ω. This extra degree of freedom in the acquisition, takes into account the local misorientations in the crystal that would reach diffraction condition for slightly different base tilts. This technique is currently integrating these tilts into single images, reconstructing the grain shapes by using regular tomography tools like filtered backprojection algorithms. It would then be extremely interesting to use the six-dimensional framework developed in this thesis, to be able to extract orientation information and possibly obtain better reconstructions for deformed grain structures in this configuration, because the appearance of slip bands and the increase of deformation in the recorded images are indeed observed in increasingly deformed grains. The reconstructed orientation information, however, will probably miss the component along the plane normal of the atomic plane aligned for diffraction, in a similar fashion to the DCT orientation reconstructions where the ω information was discarded by converting diffraction blobs into diffraction spots. This problem could be alleviated by acquiring additional projections of reflections perpendicular to the aligned hkl-vector (possibly at different energies to increase number of projections), or, if possible, to mount the sample 90 degrees rotated in order to acquire a second topotomo of an hkl-vector perpendicular to the first one. Another kind of experiments that would benefit from the use of the six-dimensional framework developed in this thesis, is the type of DCT experiments performed with laboratory or neutron sources [17]. These experiments are usually characterized by a poly-chromatic beam, and especially in the first case, by non negligible beam divergence. The development of forward and back projection functions for laboratory experiments [START_REF] Wim Van Aarle | An accurate projection model for diffraction image formation and inversion using a polychromatic cone beam[END_REF], is completely compatible with the algorithms used in the synchrotron based DCT, making the port to the laboratory DCT reconstruction code easy and straight forward, while another interesting solution could be to use energy-resolving detectors. Finally, it is worth mentioning some limitations remain for this six-dimensional framework: (a) the assumption of kinematic diffraction, (b) no correction for self-absorption or incoming beam absorption effects, (c) no correction for inhomogeneities in the beam profile. While for the cases (b) and (c), these corrections could be introduced, at least in an approximate form, without much effort, the case (a) results more problematic, because it is difficult to tackle the dynamic diffraction effects and to model the extinction inside the proposed model. Also the total variation functional used in chapters 4 and 5, might still not be the optimum solution, that exploits all the a priori knowledge available about the grain reconstruction problem. On the other hand, it proved to be fairly good choice so far, and it was able to produce encouraging results. Diffraction contrast tomography (DCT) is a technique that is mainly used in synchrotrons, with a monochromatic and parallel high energy X-ray beam, for the analysis of undeformed or nearly undeformed polycrystalline materials, described by "grains" with a certain average orientation. Here we propose a six-dimensional extension of this tomographic method, enabling the reconstruction of local orientation distributions in polycrystalline materials.

An orientation can be parametrized by the used of three numbers, so the local orientation in the grains is modelled by a six-dimensional space X6 = R3 ⊗ O3, that is the outer product between a three-dimensional real-space and three-dimensional orientation-space. The reconstruction problem is then formulated as a global optimization problem, where the reconstruction is the solution that minimizes a functional. The different functionals used in this thesis include a data fidelity term which ensures that the reconstruction is consistent with the measured diffraction data, but depending on the type of expected features in the reconstruction, additional terms like the l1-norm minimization of the solution vector, or a Total Variation operator over the real-space part of the six-dimensional voxels, are used for regularization of the heavily under-determined linear problem. While the first three chapters of this thesis introduce and explain in depth the model just described, the following two chapters study the application of this six-dimensional framework to the analysis of experimental data, hard to reconstruct within the previous, three-dimensional framework. In both cases the results are compared to electron backscatter diffraction (EBSD) surface measurements with which they turn out to be in excellent agreement with.

Chapter 4, more specifically, analyses the case of a NaCl sample, where the grains have the tendency to form subgrains. The results show an excellent agreement of the local orientations with the EBSD map, and some of the reconstructions were even performed from raw-data, bypassing the steps of segmentation and indexing of the diffraction spots. It was also demonstrated that the algorithm can reconstruct clusters of grains, by defining a bounding box in the six-dimensional space and simply using the raw-data associated to it. This means that we can surpass single grain reconstructions, by reconstructing larger regions inside a volume, as needed for highly textured samples, but more importantly that indexing is not a limiting factor any more for the reconstruction of deformed samples.

1. 1 . ORIENTATIONS 35 Figure 1 . 1 :

 13511 Figure 1.1: Coordinate systems

Figure 1 . 2 :Figure 1 . 3 :

 1213 Figure 1.2: Bragg's law

Figure 1 . 4 :

 14 Figure 1.4: Two-dimensional tomographic geometry

1. 3 .Figure 1 . 5 :

 315 Figure 1.5: Comparison of how different choices for the norms can affect the final solution.

Figure 1 . 6 :

 16 Figure 1.6: DCT projection geometry for an undeformed grain, from the sample coordinate perspective. The projection geometry is defined by the sixcoordinates (x, r), where x is the position of the center of mass of the grain, and r is the orientation.

Figure 1 . 7 :

 17 Figure 1.7:DCT projection geometry for a deformed grain, from the sample coordinate perspective. The projection geometry is defined in each voxel by the six components (x i , r i ), where x i is the position of the voxel i, and r i is its orientation.

1. 3 .Figure 1 . 8 :

 318 Figure 1.8: Comparison between the two different interpretations of a sixdimensional space that combines orientation and position information.

Figure 1 . 9 :

 19 Figure 1.9: Orientation space super sampling, as shown by the top right front voxel of the sampled grid. While the usual sampling would pick the center of the voxels associated to the volumes like in figure 1.8b, in this case we instead pick the centers of the sub-voxels and associate all of them with the same real-space volume.

  62 CHAPTER 2. RECONSTRUCTION DISCARDING ω INFORMATION(a) Global distance (b) Distance over the x -axis (c) Distance over the y-axis (d) Distance over the z -axis

Figure 2 . 1 :

 21 Figure 2.1: Distance map in degrees between the reconstruction and the original solution.

  (a) Sampling in the XY plane (b) Sampling in the XZ plane (c) Sampling in the YZ plane

Figure 2 . 2 :

 22 Figure 2.2: Orientation space sampling (red points) and phantom local orientations (blue points).

Figure 2 :

 2 Figure 2: A three-dimensional volume in real space, can be considered a single point in a three-dimensional orientation space

Figure 4 :

 4 Figure 4: Twinned Grain

Figure 6 :

 6 Figure 6: Two-dimensional orientation spread

Figure 7 :

 7 Figure 7: Reconstruction Error in degrees

Figure 3 . 1 :Figure 3 . 2 :

 3132 Figure 3.1: Comparison of the reconstruction of a five-dimensional phantom (slice illumination) between the spot reconstruction and the blob reconstruction.

Figure 3 . 3 :

 33 Figure 3.3: Orientation-space visualization for one real-space voxel in the spot reconstruction.

84CHAPTER 3 .

 3 THEORETICAL FORMULATION OF THE SIX-DIMENSIONAL MODEL (a) Orientation distribution in the XY plane (b) Orientation distribution in the XZ plane (c) Orientation distribution in the YZ plane

Figure 3 . 4 :

 34 Figure 3.4: Orientation space distribution of the voxels in the phantom (red points) and voxels of the reconstruction (blue points).

  (a) Comparison of one slice in a blob (b) Comparison of the spot generated from a blob

Figure 3 . 5 :

 35 Figure 3.5: Screenshots that show a comparison the synthetic data (left), and the same data computed from the forward projection of the reconstruction (right).

Figure 1 :

 1 Figure 1: Diffraction of one grain in a near-field diffraction imaging experiment: (a) full beam illumination, (b) 1D line beam illumination.

Figure 2 :

 2 Figure 2: Projections of: (a) uniform grain, (b) a 'deformed' grain that will give rise to distorted diffraction blobs, spreading over a range of rotation angles.

Figure 4 :

 4 Figure 4: Comparison of one slice in the 3D reconstruction of grain shapes, where the red mark identifies the original grain shape (50 × 50 voxels), and the gray-scale intensities identify the (integrated) "scattering power" coming from the voxels: (a) Phantom, (b) Single orientation SIRT, (c) Reconstruction with 341 orientations, integrated over all sampled orientations.

Figure 5 :

 5 Figure 5: Comparison of phantom and reconstruction of a 1deg orientationspread sample, using 60 blobs: (a) Phantom, (b) Reconstruction, with 341 orientations, (c) Error plot (in deg) of the reconstruction, (d) orientation-space visualization for one real-space voxel in the reconstruction

Figure 6 :

 6 Figure 6: Results for different reconstruction parameters, using 60 blobs: (a) Phantom, (b) Orientations 341, (c) Orientations 559, (d) Orientations 1241

Figure 7 :Figure 8 :

 78 Figure 7: Error plot (in deg) of reconstruction in figure (6, d)

  99

Figure 9 :

 9 Figure 9: Reconstruction of a 5deg orientation-spread sample with sub-grain boundaries, using 60 blobs: (a) Phantom, (b) Reconstruction, with 1241 orientations, (c) Inter-voxel orientation distance in deg, (d) Error plot (in deg) of reconstruction

Figure 10 :

 10 Figure 10: On the left we can see a diffstack made of Diffraction Spots, while on the right, a diffstack made of Diffraction Blobs, where around each average ω we have multiple spots, with slightly different ω, that compose Diffraction Blobs.

Figure 12 :

 12 Figure 12: Volumes project first to an intermediate projection data-structure, and then each blob slice becomes the pixel-wise sum of all the contributions from the corresponding slices in the intermediate projection data.

Figure 14 :

 14 Figure 14: (a) One slice of the original blobs, (b) Same slice with noise, (c) Extinction noise

Figure 16 :

 16 Figure 16: Comparison of phantom and reconstruction of a 1deg orientationspread sample with extinction noise, using 60 blobs: (a) Phantom, (b) Reconstruction with 341 orientations, (c) error in deg for each voxel

Figure 4 . 1 :

 41 Figure 4.1: Reconstruction of grain 63, which exhibits a sub-grain structure, that could not be reconstructed properly, through the standard DCT reconstruction route, but that was completely reconstructed by performing an extension of the six-dimensional reconstruction bounding box, and using raw images.

Figure 1 :

 1 Figure 1: Comparison of EBSD surface mapping with the different reconstruction approaches for full-field X-ray orientation microscopy discussed in this work.

  (a) Experimental spot image (b) Theoretically fwd-projected image

Figure 2 :

 2 Figure 2: Comparison of the experimental images for an [2 2 2] reflection at θ = 6.21 degrees, η = 112.112 degrees, with a ∆ω = 6.7 degrees (67 images), with the same forward-projected spot from the result of the reconstruction.

Figure 3 :

 3 Figure 3: Reconstruction of a grain cluster using the extended 6D approach. (a, b, c) same slice through the 3D reconstruction showing: (a) Intra-granular Misorientation, (b) Kernel Average Misorientation (c) inverse pole figure colour coding scheme revealing the presence of sub-grains and small angle boundaries from skeletonization of (a) (gray: ≥ 0.5 • , black: ≥ 2 • ), (d) iso-surface of the orientation sub-space reconstructed for the clustered region. Red points indicate sub-grain orientations which had been successfully identified using conventional indexing procedures, along with their corresponding grain ID.

  (a) Voxel-wise misorientation (b) Grain-boundary differences (c) Histogram of misorietations (d) Histogram of boundary distances

Figure 4 :

 4 Figure 4: Quantitative comparison between EBSD map and the selected surface slice in the 6D-DCT volume, where we can find: a) voxel-wise distance in degrees between the orientations found by EBSD and 6D-DCT, b) Overlay of grain-boundaries from EBSD (red) and 6D-DCT (yellow), c) Histogram of misorientations in a), d) Histogram of pixel distances in b).

Figure 1 :

 1 Figure 1: IPF Color key used for the orientation coloring in figure (1).

Figure 1 :Figure 2 :Figure 3 :Figure 4 :

 1234 Figure 1: (a) Orientation map of a longitudinal sample cross section parallel to the rolling direction, recorded by EBSD (IPF colour code). The dashed lines indicate the sample sub-volume characterized in the DCT scan, (b) Quality map (Band Contrast) of the map in (a) on the byte scale. Deformation twins are clearly discernible on the EBSD map whereas slip traces are more easily discerned on the quality map (e.g. inside the grain labeled with the letter C).

Figure 5 :Figure 6 :

 56 Figure 5: Intragranular misorientation of the grain labelled as C on the surface shown in figure (1). EBSD and DCT maps have voxel sizes of 0.4 µm and 1.4 µm, therefore the number of voxels shown on the figure axes are different, while the scale bars are similar. The pictures refer to: (a) EBSD, (b) 6D-DCT.

Figure 7 :

 7 Figure 7: Slice in the full reconstructed volume: (a) Inverse Pole Figure coloring of the grain orientations, using the color coding from the legend in figure (1, a), (b) intensity assigned by the 6D-DCT reconstruction in each real-space voxel.

x ||Ax -b|| 2 2 ( 6

 26 .1a) subject to: x 0, W x = y (6.1b) and:x * = argmin x ||Ax -b|| 2 2 + λTV (x)(6.2a)subject to: x 0, W x = y (6.2b)

  -field X-ray orientation imaging using convex optimization and a discrete representation of sixdimensional position -orientation space NATURE : Doctorat Numéro d'ordre : 2015ISAL0095 Ecole doctorale : Matériaux Lyon RESUME :

Chapter 5 ,

 5 instead, analyses the case of a 1% plastically deformed microstructure in Ti 4% Al alloy, giving rise to the formation of deformation twins. Twins are crystallographic domains that are very difficult to deal with, because they tend to be thin, plate-like regions, which produce weak signals on the detector, and are more sensitive to factors like the point spread function (PSF) of the detector. The six-dimensional framework was able to improve the reconstruction of twinned grains, compared to the previous three-dimensional implementation.MOTS-CLÉS : six-dimensions, diffraction contrast tomography, DCT, orientation imaging, three-dimensional X-ray diffraction microscopy, convex optimization, total-variation minimization, l1-norm minimization Laboratoire (s) de recherche: iMinds-VisionLab Antwerp; MATEIS, INSA-Lyon UMR CNRS 5510 ; ESRF Grenoble Directeur de thèse : Wolfgang LUDWIG et K. Joost BATENBURG et Jan SIJBERS Président de jury : Olivier CASTELNAU Composition du jury : H. F. Poulsen (Rapporteur) A. Borbely (Rapporteur) C. Le Bourlot (Examinateur) O. Castelnau (Examinateur) K. J. Batenburg (Directeur de thèse) J. Sijbers (Directeur de thèse) W. Ludwig (Directeur de thèse)

  Comparison of the reconstruction error in degrees between the spot reconstruction and the blob reconstruction . . . . . . . . . . 82 3.3 Orientation-space visualization for one real-space voxel in the spot reconstruction. . . . . . . . . . . . . . . . . . . . . . . . . . 83 3.4 Orientation space distribution of the voxels in the phantom (red points) and voxels of the reconstruction (blue points). . . . . . . 84 3.5 Screenshots that show a comparison the synthetic data (left), and the same data computed from the forward projection of the reconstruction (right). . . . . . . . . . . . . . . . . . . . . . . . . 84 4.1 Reconstruction of grain 63, which exhibits a sub-grain structure, that could not be reconstructed properly, through the standard DCT reconstruction route, but that was completely reconstructed by performing an extension of the six-dimensional reconstruction bounding box, and using raw images. . . . . . . . . . . . . . . . . 113

	LIST OF FIGURES
	3.2 Chapter 1
	Introduction
	1.1 Orientations
	. .
	2.2 Orientation space sampling (red points) and phantom local ori-
	entations (blue points). . . . . . . . . . . . . . . . . . . . . . . .
	3.1 Comparison of the reconstruction of a five-dimensional phantom
	(slice illumination) between the spot reconstruction and the blob

1.1 Coordinate systems . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Bragg's law . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Diffraction geometry of a DCT experiment, where, for simplicity, the detector is supposed to be on the plane perpendicular to the incoming beam direction . . . . . . . . . . . . . . . . . . . . . . . 1.4 Two-dimensional tomographic geometry . . . . . . . . . . . . . . 1.5 Comparison of how different choices for the norms can affect the final solution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.6 DCT projection geometry for an undeformed grain, from the sample coordinate perspective. The projection geometry is defined by the six-coordinates (x, r), where x is the position of the center of mass of the grain, and r is the orientation. . . . . . . . . . . . 1.7 DCT projection geometry for a deformed grain, from the sample coordinate perspective. The projection geometry is defined in each voxel by the six components (x i , r i ), where x i is the position of the voxel i, and r i is its orientation. . . . . . . . . . . . . . . . 1.8 Comparison between the two different interpretations of a sixdimensional space that combines orientation and position information. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.9 Orientation space super sampling, as shown by the top right front voxel of the sampled grid. While the usual sampling would pick the center of the voxels associated to the volumes like in figure 1.8b, in this case we instead pick the centers of the sub-voxels and associate all of them with the same real-space volume. . . . . 2.1 Distance map in degrees between the reconstruction and the original solution. . . . . . . . . . . . . . . . . . . . . . . . . . . . reconstruction. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  acknowledges financial support of the Netherlands Organization for Scientific Research (NWO), project nr. 639.072.005. NV and WL acknowledge financial support of the French National Research Agency (ANR), project: ANR 2010 BLAN 0935. The authors acknowledge COST Action MP1207 for networking support.

  b,f , Lorenzo Valzania b , Gaurav Singh f , Michael Preuss f , Kees Joost Batenburg d,c,e , and Wolfgang Ludwig a,b a MATEIS, INSA Lyon, Université de Lyon, UMR5510 CNRS, France b ESRF Grenoble, France c iMinds-Vision Lab, Univ. Antwerp, Belgium d CWI, Amsterdam, Netherlands e Universiteit Leiden, Netherlands f University of Manchester, United Kingdom

	February 5, 2016

Cette thèse est accessible à l'adresse : http://theses.insa-lyon.fr/publication/2015ISAL0095/these.pdf © [N. R. Viganò], [2015], INSA Lyon, tous droits réservés

Cette thèse est accessible à l'adresse : http://theses.insa-lyon.fr/publication/2015ISAL0095/these.pdf © [N. R. Viganò], [2015], INSA Lyon, tous droits réservés 1.2. DIFFRACTION CONTRAST TOMOGRAPHY

× 6 × 6 cubic lattice, with another 5 × 5 × 5 cubic lattice occupying the centers of the cells of the bigger lattice. The total volume covered by this BCC lattice in orientation-space is big enough to envelop all the orientations in the sample, but not too big, in order to not loose resolution (each edge is ∼ 1.1deg).

A typical value at ID11, ESRF would be 400 3 voxels within

h, corresponding to a gain of a factor 10 and 100 in acquisition time as compared to currently reported acquisition times for slice beam and point beam scanning techniques, respectively. With progress in instrumentation, the slicing and scanning approaches may become competitive in time resolution. Since these scanning approaches offer less convoluted and hence easier to analyse diffraction data, adaptations of the global (6D) optimization framework to these acquisition modes shall be envisaged.

Since the relative volume fraction of parent and twin are not known beforehand, a joint reconstruction of parent and twin can no longer rely on the re-normalization of the diffraction blob intensities used in previous implementations of the reconstruction algorithm where all the identified and segmented reflections would be renormalized to the same intensity.

(Another parametrization of the diffraction "blob" positions on the detector, instead of (u, v), is given in chapter 3 of[13], where the aperture (Bragg) angle between the incoming

the intensity value obtained from projecting the 6D space onto the 3D real space grid can be interpreted as the local "scattering power" of the material

Acknowledgments 7 Acknowledgments

KJB acknowledges financial support of the Netherlands Organization for Scientific Research (NWO), project nr. 639.072.005. NV and WL acknowledge financial support of the French National Research Agency (ANR), project: ANR 2010 BLAN 0935. NV and WL acknowledge Péter Reischig for having created a tool for generating synthetic test datasets, which was used in the creation of the tests cases for this article. The authors acknowledge COST Action MP1207 for networking support.

Acknowledgments

The authors acknowledge the ESRF for providing access to beamtime and computing infrastructure as well as Eva Héripré, Ababacar Gaye, Mathieu Bourcier and Fanny Jambon for their help during the experiment campaigns. NV, WL, MB and AD acknowledge financial support of the French National Research Agency (ANR), project: ANR 2010 BLAN 0935 and the German Research Foundation (DFG), project SPP1466. KJB acknowledges financial support of the Netherlands Organization for Scientific Research (NWO), project nr. 639.072.005. The authors acknowledge COST Action MP1207 for networking support.

Acknowledgments

The authors acknowledge ESRF for providing access to inhouse research beamtime and computing infrastructure. We would further like to thank Arnas Fitzner for providing the material. NV and WL acknowledge financial support of the French National Research Agency (ANR), project: ANR 2010 BLAN 0935. MP and LN acknowledge financial support from EPSRC, UK (EP/F020910/1). KJB acknowledges financial support of the Netherlands Organization for Scientific Research (NWO), project nr. 639.072.005. The authors acknowledge COST Action MP1207 for networking support.

A Implementation Details

Let us now have a look at the deta-structures first introduced in section 2.3.

Solution vector

It is the data-structure that is responsible for representing the 6-dimensional space in the computer memory. Its inner structure reflects the schematic representation in figure 3: we decided indeed to use a linear array of three-dimensional volumes, where each of these volumes has the same size and represents a specific point in orientation space. The number of volumes is equal to the number of sampled orientations.

Author Contributions

NV carried out the X-ray data analysis and wrote the manuscript jointly with WL. AD, MB and WL performed the X-ray experiment and SH, AD and AT prepared the sample and performed the EBSD analysis. WL and JB supervised the project. All authors have reviewed the manuscript.

Additional Information

The authors declare no competing financial interests.

A Theoretical determination of twin variants

Twin orientations are related to the parent orientation by means specific pairs of rotation axis and associated rotation angles. Given the set of N symmetry operators S = {S i } with i = 1, . . . , N , in the form of unitary rotation matrices, and the rotation axis a which defines the type of axis for a given twinning mode, we obtain all the possible twinning rotation axes for the said mode, by simply taking the unique results from the set A = {S i a}, that will be another set called A . If then we take the orientation matrix g p associated to the Rodrigues vector r p which is the orientation of the parent grain, the orientation matrices of the twin variants will then be found by the simple operation:

for each element i in A , where R (a i , α) is nothing else than the rotation matrix defined by the axis a i and the angle α.

The same operation in equation ( 8) could be carried out in Rodrigues space, by using the multiplication formula from [6] for successive rotations:

Appendix A

Prediction of scattering ωs

We will now show the solution to equation 1.25. The first step is the decomposition of the rotation matrix Ω ω in its three components from the Rodrigues formula, where these three matrices only depend on the rotation axis, while the final rotation matrix is given by:

As a result, substituting the right side of equation A.1 into equation 1.25, leads to:

If we then call A = bT Ω cos g -1 h, B = bT Ω sin g -1 h, and C = bT Ω const g -1 h, equation A.2 can then be rewritten as a set of two equations:

These two equations can be then recast as second order equations and they would give rise to two solutions each, giving a total of four possible omegas.

There are many different ways to solve them, but here we will present the tangent based solution. We first rewrite equations A.4 into one single equation by creating C 1,2 = C ∓ sin θ, and remembering the equivalence formulas: