Elisa Ruhland

Caroline Bund

Hassiba Outilaft

Martial Piotto

Izzie-Jacques Namer
email: izzie.jacques.namer@chru-strasbourg.fr

A metabolic database for biomedical studies of biopsy specimens by high-resolution magic angle spinning nuclear MR: A qualitative and quantitative tool

Keywords: biopsies, database, HRMAS NMR, metabolomics

Purpose:

The aim of this study is to generate a metabolic database for biomedical studies of biopsy specimens by high-resolution magic angle spinning (HRMAS) nuclear MR (NMR). Methods: Seventy-six metabolites, classically found in human biopsy samples, were prepared in aqueous solution at a known concentration and analyzed by HRMAS NMR. The spectra were recorded under the same conditions as the ones used for the analysis of biopsy specimens routinely performed in our hospital. Results: For each metabolite, a complete set of NMR spectra (1D 1 H, 1D 1 H-CPMG, 2D J-Resolved, 2D TOCSY, and 2D 1 H-13 C HSQC) was recorded at 500 MHz and 277 K. All spectra were manually assigned using the information contained in the different spectra and existing databases. Experiments to measure the T 1 and the T 2 of the different protons present in the 76 metabolites were also recorded.

Conclusion:

This new HRMAS metabolic database is a useful tool for all scientists working on human biopsy specimens, particularly in the field of oncology. It will make the identification of metabolites in biopsy specimens faster and more reliable.

Additionally, the knowledge of the T 1 and T 2 values will allow to obtain a more accurate quantification of the metabolites present in biopsy specimens.

| 63

RUHLAND et AL.

is particularly useful to characterize the biology of the tumor, monitor treatments, and predict toxicity. [START_REF] Claudino | Metabolomics: available results, current research projects in breast cancer, and future applications[END_REF] Nuclear MR (NMR) metabolomics of tissues has been used to study several cancers [START_REF] Bathen | Magnetic resonance metabolomics of intact tissue: a biotechnological tool in cancer diagnostics and treatment evaluation[END_REF] and has the potential to enable personalized treatment.

High-resolution magic angle spinning (HRMAS) NMR spectroscopy is particularly adapted to study small samples of unprocessed tissue specimens for medical indications such as tumor analysis [START_REF] Erb | Toward improved grading of malignancy in oligodendrogliomas using metabolomics[END_REF] and transplantation. [START_REF] Faitot | Impact of real-time metabolomics in liver transplantation: graft evaluation and donor-recipient matching[END_REF] The preparation of biopsy samples is straightforward and fast, [START_REF] Battini | Metabolomics approaches in pancreatic adenocarcinoma: tumor metabolism profiling predicts clinical outcome of patients[END_REF] obviating the need for lengthy chemical extraction procedures. HRMAS NMR spectroscopy can identify and quantify several metabolites from 1D [START_REF] Griffin | Metabolic profiles of cancer cells[END_REF] H NMR spectra with excellent resolution and signal-to-noise ratio. In clinical routine, this technique could eventually be used to characterize biopsies and find predictive and prognostic factors during surgery. In our hospital, samples are transported rapidly by a pneumatic tube from the operating room to the room where the NMR spectrometer is located. Acquisition of HRMAS data takes only around 20 minutes. [START_REF] Erb | Toward improved grading of malignancy in oligodendrogliomas using metabolomics[END_REF][START_REF] Battini | Metabolomics approaches in pancreatic adenocarcinoma: tumor metabolism profiling predicts clinical outcome of patients[END_REF][START_REF] Battini | High-resolution magic angle spinning (1)H nuclear magnetic resonance spectroscopy metabolomics of hyperfunctioning parathyroid glands[END_REF][START_REF] Imperiale | In vivo detection of catecholamines by magnetic resonance spectroscopy: a potential specific biomarker for the diagnosis of pheochromocytoma[END_REF] HRMAS NMR has the potential to provide important medical benefits; however, metabolite identification and quantification needs to be fast and reliable. Identification is usually based on the literature and existing databases that have been developed mainly for urine, plasma and CSF [START_REF] Wishart | HMDB 3.0-The Human Metabolome Database in 2013[END_REF] (HMDB: http:// www.hmdb.ca/).10,11 A specific HRMAS database for the study of human biological tissues is still lacking and needs to be developed. Moreover, quantification is essential for network analysis and discovery of new biomarkers. The aim of this study is to generate a metabolic HRMAS database for the 76 most common metabolites found in human tissues.

| METHODS

2.1 | Preparation of metabolite solutions 2.1.1 | Pure metabolite standards A total of 76 metabolite solutions were prepared to build the metabolic database. For each metabolite, a solution containing approximately 3 mM of the pure metabolite was prepared in D 2 O containing 9.4 mM of phosphate-buffered saline (PBS) buffer (pH 7.4) and 0.125% trisodium phosphate (TSP). Twenty-four microliters of this solution were placed in a 25-µL disposable insert in KeLF. The number of nmoles actually present in the inserts of each standard solution was determined exactly using the Eretic 12 method (Table 1

| Biopsies preparation

The protocol is the one classically used in our laboratory to prepare biopsy specimens: All tissue samples were collected during surgery immediately after removal of the tumor. Specimens were either frozen in liquid nitrogen in an appropriate Dewar in the operating room or transferred from the block to the spectrometer room by a pneumatic system. Samples were then prepared at a temperature of -20°C. Amount of tissue used ranged from 15 to 20 mg. Each tissue sample was placed in a 25 μL disposable insert using a biopsy punch. A volume of 10 µL of deuterium oxide was then added to each biopsy insert. After a centrifugation step, deuterium oxide was added again to fill the insert. Finally, inserts were stored at -80°C. until the HRMAS analysis was performed. The insert was then placed in a 4-mm ZrO 2 rotor and inserted into the NMR spectrometer. All assays were performed at 4°C. After the analysis, inserts were stored back at -80°C.

| HRMAS NMR spectroscopic analysis

All experiments were conducted on a Bruker Avance III 500 MHz NMR spectrometer (Bruker BioSpin, Billerica, MA) installed at Strasbourg University Hospitals (www. mnms-platform.com). The spectrometer is equipped with a triple-resonance ( 1 H, [START_REF] Aue | Homonuclear broad-band decoupling and 2-dimensional J-resolved NMR spectroscopy[END_REF] C, and 31 P) HRMAS probe. The temperature was set to 277 K for all acquisitions. Spinning speed was set to 3502 Hz in order to keep the rotation sidebands out of the spectrum window and to minimize sample degradation. For each standard metabolite solution, a complete set of NMR experiments was recorded under exactly the same experimental conditions as the ones used in our laboratory when analyzing human tissues. [START_REF] Battini | Metabolomics approaches in pancreatic adenocarcinoma: tumor metabolism profiling predicts clinical outcome of patients[END_REF][START_REF] Imperiale | In vivo detection of catecholamines by magnetic resonance spectroscopy: a potential specific biomarker for the diagnosis of pheochromocytoma[END_REF] A 1-dimensional 1 H CPMG spectrum was first recorded and followed by a 1-pulse [START_REF] Faitot | Impact of real-time metabolomics in liver transplantation: graft evaluation and donor-recipient matching[END_REF][START_REF] Battini | Metabolomics approaches in pancreatic adenocarcinoma: tumor metabolism profiling predicts clinical outcome of patients[END_REF] experiment in order to check the stability of the metabolites. For all metabolites, the following 2D experiments were then recorded: J-Resolved, [START_REF] Aue | Homonuclear broad-band decoupling and 2-dimensional J-resolved NMR spectroscopy[END_REF] TOCSY, [START_REF] Braunschweiler | Coherence transfer by isotropic mixing: application to proton correlation spectroscopy[END_REF][START_REF] Kupce | Adiabatic TOCSY MAS in liquids[END_REF] T A B L E 1 • Group has been defined as cluster 0 in Chenomx. and 1 H- [START_REF] Aue | Homonuclear broad-band decoupling and 2-dimensional J-resolved NMR spectroscopy[END_REF] C HSQC. [START_REF] Davis | Experiments for recording pure-absorption heteronuclear correlation spectra using pulsed field gradients[END_REF] For metabolites containing phosphorus, additional experiments were also recorded: 1D 31 P with 1 H decoupling and 2D 1 H-31 P HMBC experiment (see Supporting Information).

| Relaxation times measurements were performed after the two 1D sequences

Transverse relaxation times were measured using a pseudo 2D 1 H Carr-Purcell-Meiboom-Gill (CPMG) pulse sequence with presaturation. This experiment consists in the acquisition of several 1D spectra with an increasing length of CPMG times. The experiment was performed with the following parameters: spectral width 14.0 ppm, number of points 32 k, relaxation delay 2 seconds, acquisition time 2.33 seconds, and 64 scans. Total acquisition time was of 1 hour 26 minutes.

Longitudinal relaxation times were measured using a pseudo 2D 1 H inversion recovery pulse sequence with water suppression, which consists in the acquisition of several 1D spectra with an increasing relaxation delay after the 180° pulse was used. A CPMG module of 93 ms was added just before the acquisition time in order to reduce broad spectral features, which arise from lipids and macromolecules. The experiment was performed with the following parameters: spectral width 14.0 ppm, number of points 32 k, relaxation delay 8 seconds, acquisition time 2.33 seconds, and 64 scans. Total acquisition time was 2 hours and 6 minutes.

| HRMAS NMR data processing

All spectra of the metabolic database were processed using Topspin 3.5 (Bruker BioSpin) and calibrated with respect to the TSP signal at 0.00 ppm in 1 H and at -0.25 ppm in [START_REF] Aue | Homonuclear broad-band decoupling and 2-dimensional J-resolved NMR spectroscopy[END_REF] C. Phosphorus spectra were referenced to the 31 P signal of the PBS (H 2 PO 4 -/HPO 4 2-) buffer at 5.36 ppm.

| Processing and analysis of 2D relaxation data

For both T 1 and T 2 analysis, the free induction decays were first multiplied with a cosine squared weighting function before Fourier transformation. Data were then carefully phased and baseline corrected. The resulting data were analyzed using the Dynamics Center 2.2.4 package, which is a module of Topspin 3.5 and which is able to perform advanced mono-and biexponential curve fitting. This software fits relaxation data using a combination of Simplex and Levenberg-Marquardt algorithms.

and T 2 effects for a more accurate quantification in CPMG experiments

The signal measured after a CPMG sequence when taking into account T 1 and T 2 relaxation is given by (Equation 1):

where S o is the true intensity of the signal. The signal loss is therefore given by (Equation 2):

where t a is the acquisition time plus the recycle time and t b is the length of the CPMG element. Under our experimental conditions: t a = 4.33 seconds and t b = 93.71 ms.

The correction factor that should be used to correct for the relaxation effects is therefore (Equation 3):

| Human brain tumor specimens

Six human brain specimens (normal cerebral cortex [n = 2], oligodendroglioma [n = 2], and glioblastoma [n = 2]) were prepared following our standard protocol in order to evaluate the performances of our database for metabolite assignment. Additionally, the T 1 and T 2 relaxation times of 17 brain metabolites (20 chemical groups) were determined in these samples. These relaxation times were measured to confirm the applicability of our database to biological samples.

| RESULTS

For each of the 76 metabolites contained in the database, the full set of experiments described in Methods was acquired. The 1D and 2D spectra were manually assigned using the information contained in the different spectra as well as the Human Metabolome Database 9 and simulated NMR spectra. [START_REF] Binev | Prediction of 1H NMR coupling constants with associative neural networks trained for chemical shifts[END_REF] The complete list of 1 H/ [START_REF] Aue | Homonuclear broad-band decoupling and 2-dimensional J-resolved NMR spectroscopy[END_REF] C chemical shifts and 1 H-1 H J-couplings for the 76 metabolites is presented in Table 1. It also lists the T 1 and T 2 values of the different protons of (1)

S Measured = S o 1 -exp - t a T 1 × exp - t b T 2 (2) S Loss = S o 1 -1 -exp - t a T 1 exp - t b T 2 (3) C Factor = 1∕ 1 -exp - t a T 1 × exp - t b T 2
the molecules in the metabolic database. An example of the T 1 /T 2 measurement for the different signals of ascorbate is shown in Figure 1. As expected for small molecules in solution, the values of the T 1 and T 2 are very similar. However, some protons like the Lactate_CH exhibit a distinct behavior, with T 2 being much shorter than T 1 . These relaxation parameters are particularly important because they allow to estimate the amount of magnetization lost because of relaxation and to compute a correction factor. The comparison of the T 1 /T 2 values obtained for the metabolites in solution and in biopsy samples is presented in Table 2. For all metabolites, both T 1 and T 2 values in biopsy specimens are slightly shorter than in solution. On average, the T 1 values obtained in biopsy samples are 17.3% shorter whereas the T2 values are 27.2% shorter. The average loss of magnetization attributed to T 1 /T 2 relaxation effects during a CPMG experiment is computed to be 16% in solution and 17% in biopsy specimens. This result is quite remarkable and provides, for the first time, an estimate of the correction factor that should be used to obtain a more accurate quantification. The fact that the correction factors are so similar, despite slight variations in T 1 T 2 , can be explained using Equation . A shorter T 1 will decrease signal loss, whereas a shorter T 2 will increase it. The simultaneous decrease of T 1 and T 2 in biopsy specimens explains that similar signal losses are observed. These values are applicable, provided our own experimental conditions are followed (i.e., CPMG block of 93 ms and recovery delay of 4.33 seconds). However, the T 1 and T 2 values presented in Table 1 allow to estimate the loss in signal for any of these 2 experimental parameters. As a validation sample, we have used 4 mixtures at different concentrations of 4 metabolites. The results presented in Figure 2 show that a high precision in the concentrations can be obtained, provided the corrections attributed to relaxation effects are taken into account.

F I G U R E 1 A 500-MHz

| DISCUSSION

HRMAS NMR metabolomics is a promising technique to evaluate the metabolic status of the tissue of a patient. Depending on the tissue and the disease studied, different metabolic pathways and different metabolites may be affected. [START_REF] Battini | Metabolomics approaches in pancreatic adenocarcinoma: tumor metabolism profiling predicts clinical outcome of patients[END_REF][START_REF] Battini | High-resolution magic angle spinning (1)H nuclear magnetic resonance spectroscopy metabolomics of hyperfunctioning parathyroid glands[END_REF][START_REF] Imperiale | In vivo detection of catecholamines by magnetic resonance spectroscopy: a potential specific biomarker for the diagnosis of pheochromocytoma[END_REF] For instance, alterations in the following pathways have been studied by HRMAS: glycolysis, tricarboxylic acid glutamate/glutamine cycle, 1-carbon metabolism, membrane metabolism, and amino acid synthesis. HRMAS NMR metabolomics has been shown to be potentially beneficial in personalized treatment in the following fields of medicine [START_REF] Armitage | Monitoring cancer prognosis, diagnosis and treatment efficacy using metabolomics and lipidomics[END_REF][START_REF] Vander Heiden | Targeting cancer metabolism: a therapeutic window opens[END_REF][START_REF] Everett | NMR-based pharmacometabonomics: a new paradigm for personalised or precision medicine[END_REF] :

(a) Evaluation of the quality of grafts in liver or lung transplants. 5,21,22 (b) Detection of oncometabolites for identification of DNA mutations: 2-hydroxyglutarate in glioma for isocitrate dehydrogenase mutation, [START_REF] Choi | 2-hydroxyglutarate detection by magnetic resonance spectroscopy in IDH-mutated patients with gliomas[END_REF] succinate in pheochromocytomas to distinguish succinate dehydrogenase (SDHx)-and non-SDHx-related tumors. 24 (c) Detection of biomarkers to evaluate tumor aggressiveness. Past studies have shown which metabolic pathways are involved in tumor growth. [START_REF] Bathen | Magnetic resonance metabolomics of intact tissue: a biotechnological tool in cancer diagnostics and treatment evaluation[END_REF][START_REF] Loening | Quantification of phosphocholine and glycerophosphocholine with 31P edited 1H NMR spectroscopy[END_REF] In oligodendrogliomas, the degree of tumor hypoxia obtained at an early state of disease was better correlated with the patient prognosis than the morphological criteria. [START_REF] Erb | Toward improved grading of malignancy in oligodendrogliomas using metabolomics[END_REF] Distinct metabolic signatures were found for lung adenocarcinoma (related to phospholipid metabolism and protein catabolism) and squamous-cell carcinoma (glycolytic and glutaminolytic profiles). [START_REF] Rocha | NMR metabolomics of human lung tumours reveals distinct metabolic signatures for adenocarcinoma and squamous cell carcinoma[END_REF] In pancreatic adenocarcinoma, the assessment of ethanolamine concentration can be clinically relevant as a single metabolic biomarker to predict long-term survival. 6 (d) Evaluation of tumor exeresis during surgery: The intraoperative metabolic analysis of biopsy specimens can help better characterize tumor margins. This novel application has been presented as a proof of concept in colon cancer. The value of T 1 and T 2 for the pure metabolite is also reported. The last 2 columns show the signal loss expected during a CPMG experiment for the 17 metabolites in a biopsy and in a pure solution, taking into account both T 1 and T 2 effects.

The simplified metabolite pathways used in our laboratory is presented in Supporting Information Figure S1. The analysis of the variations in concentration of the different metabolites is usually performed by multivariate or network analysis. [START_REF] Cicek | ADEMA: an algorithm to determine expected metabolite level alterations using mutual information[END_REF] This approach allows to study the up-and downregulation of different pathways in order to find potential biomarkers. In this context, rapid and robust identification as well as reliable quantification are critical.

When performing a study by HRMAS NMR on a large cohort of samples, the detailed metabolite assignment is usually performed only on a few samples. For these samples, the combination of 1D and 2D NMR spectra is particularly powerful, and the database presented in this article makes this assignment relatively straightforward. A unique feature of our metabolic database is that it contains all the spectra used classically to study biopsy specimens: 1 H 1D, Several softwares, like Metabominer [START_REF] Xia | MetaboMinersemi-automated identification of metabolites from 2D NMR spectra of complex biofluids[END_REF] and AssureNMR (Bruker BioSpin) can provide a semiautomated identification using J-Resolved, TOCSY, and 1 H- [START_REF] Aue | Homonuclear broad-band decoupling and 2-dimensional J-resolved NMR spectroscopy[END_REF] C HSQC spectra. [START_REF] Xia | MetaboMinersemi-automated identification of metabolites from 2D NMR spectra of complex biofluids[END_REF] Once the different peaks have been assigned in these samples, the assignment can simply be reported in the 1D spectrum of the remaining samples. Because the conditions for sample preparation and data acquisition follow a strict protocol, there is little variation in peaks position attributed to pH or temperature. Our database was applied to analyze multiple biopsy specimens obtained from tissue banks or from routine surgery. It is clearly a powerful tool that helps provide an accurate and fast metabolic profile. To illustrate this point, Figure 3A shows a typical example of the use of the database when assigning the 1D 1 H CPMG spectrum of a biopsy sample. The CPMG is the experiment that is used systematically to characterize biopsy specimens. In this example, the assignment of Allocystationine is highlighted. In case of severe overlap, the assignment can also be confirmed by superimposing the different 2D spectra of the biopsy sample with the corresponding spectra of the database. Figure 3B illustrates this approach by showing the assignment of Allocystationine in the de 2D J-Resolved spectrum of a human biopsy. This approach allows to assign the different signals contained in a biopsy sample with a very high degree of confidence.

For biomedical applications, it is essential to obtain quantitative information on metabolites present in tissues, that is, in the cells, extracellular space, and microenvironment. This information is usually presented as the number of nmoles of metabolites per mg of biopsy material. The reliability of the analysis will depend on the quality of the tissue sample controlled in histopathology (with a preference for the HRMAS sample itself or the mirror sample). The T 1 and T 2 data presented in this article allow to improve the quantification of metabolites in biopsies by computing, for each chemical group of a given metabolite, a correction factor that takes into account relaxation losses during the CPMG experiment. Under our experimental CPMG conditions, this loss of signal ranges from 11% to 23%. Similar relaxation times measurements have been performed in vivo at 17.2T in the rat brain to improve quantification of metabolites. [START_REF] Lopez-Kolkovsky | Metabolite and macromolecule T1 and T2 relaxation times in the rat brain in vivo at 17.2T[END_REF] A very important feature of our HRMAS database is that it contains spectra recorded on a known quantity of metabolites. This implies that the CPMG spectrum of the database can be used to directly estimate the amounts of metabolites contained in the CPMG spectrum of a biopsy sample. The scaling factor used to match the peaks in the biopsy spectrum can be used to calculate the concentration of that specific metabolite. Techniques of quantification of metabolites in biopsy specimens are becoming very performant; however, they still need to be improved, notably in terms of speed and automation. In our laboratory, these quantification routines are currently being implemented in the Chenomx and AssureNMR software. This database has the potential to become a reference to identify metabolite peaks in tissues outside our laboratory, provided data are collected under the exact experimental conditions. It should be useful at both 500 and 600 MHz, which are the 2 magnetic fields most commonly used in this type of study.

| CONCLUSION

HRMAS NMR metabolomics is becoming an important technique in clinical research and could potentially improve personalized medicine. The database presented in this study considerably simplifies the assignment of signals in the spectra of biopsy specimens. The spectra of the database are all provided along with the number of nmoles used for the standard sample, which makes them useful for quantification applications. Additionally, the loss of magnetization attributed to T 1 and T 2 effects observed during CPMG experiments run under our experimental conditions is reported for metabolites in solution (average loss 16%) and in biopsy specimens (average loss 17%).

T A B L E 1

 1 Abbreviations: s = singlet; d = doublet; t = triplet; m = multiplet; dd = doublet of doublet ; q = quadruplet. * Referenced to TSP (3-(trimethylsilyl) propionic-2.2, 3.3-d4 acid, sodium salt) at 0.00 ppm in 1 H and -0.25 ppm in 13 C.

2. 4 |aspects 2 . 4 . 1 |

 4241 Quantitative Correction factors to take into account T 1

  1D 1 H CPMG HRMAS spectrum of ascorbate along with graphs showing the measurement of T 1 and T 2 for the 3 different chemical groups (peak 3 at 4.52 ppm: b and c; peak 2 at 4.02 ppm: d and e; peak 1 at 3.74 ppm: f and g)

F

  I G U R E 2 A 500-MHz 1D 1 H CPMG HRMAS spectrum of a mixture containing 4 metabolites (acetate, ascorbate, choline, and hypotaurine) (A) described in Methods. The measurement of T 1 (B) and T 2 (C) of hypotaurine at 3.35 ppm is shown. The graph (D) shows the correlation between the number of nanomoles measured and expected in the 4 different solutions for hypotaurine. The nanomoles measured values have been corrected for T 1 and T 2 effects | 81 RUHLAND et AL.

1 H

 1 1D CMPG, 2D J-Resolved, 2D TOCSY, and 2D HSQC.

F

  I G U R E 3 (A) A 500-MHz 1D 1 H CPMG HRMAS spectrum of a human brain tumor sample (right temporal glioblastoma) and several spectra of the database superimposed on top. The resonances of allocystathionine are highlighted with blue squares. The green squares show the spectral region at 3.3 ppm, where there is a strong signal overlap that will require additional 2D experiments. (B) A 500-MHz 2D J-resolved HRMAS spectrum of a human brain tumor sample (right temporal glioblastoma) superimposed with the allocystathionine spectrum (the spectrum has been slightly shifted to the right in order to visualize the good agreement between the peaks)

2.1.2 | Standard metabolites mixture for the validation of concentration measurements

	concentration measurements, taking into account T 1 and T 2
	effects.
	The content of the 4 solutions is listed below:
	• Mixture 1: acetate, ascorbate, choline, and hypotaurine
	(19.68, 19.44, 19.92, and 20.4 nmoles, respectively);
	• Mixture 2: acetate, ascorbate, choline, and hypotaurine
	(32.16, 31.92, 32.64, and 33.36 nmoles, respectively);
	• Mixture 3: acetate ascorbate, choline and hypotaurine
	(41.04, 40.8, 41.52, and 42.48 nmoles, respectively); and
	• Mixture 4: acetate, ascorbate, choline, and hypotaurine
	(53.52, 53.28, 54.24, and 53.04 nmoles, respectively).
	Four solutions containing 4 metabolites were prepared in the
	same D 2 O/PBS buffer at different concentrations to validate

  15222594, 2019, 1, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/mrm.27696 by Cochrane France, Wiley Online Library on [16/10/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

	RUHLAND et AL.	| 77

  15222594, 2019, 1, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/mrm.27696 by Cochrane France, Wiley Online Library on [16/10/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

	RUHLAND et AL.	| 79

  [START_REF] Piotto | Towards real-time metabolic profiling of a biopsy specimen during a surgical operation by 1H high resolution magic angle spinning nuclear magnetic resonance: a case report[END_REF]

  T 1 and T 2 relaxation times measured for 17 metabolites and 20 chemical groups on 6 different biopsies: normal cerebral cortex

	(n = 2), oligodendroglioma (n = 2), glioblastoma (n = 2)
							Solution
					Solution T 2	Biopsy CPMG	CPMG
		Biopsy T 1 (s)	Solution T 1 (s)	Biopsy T 2 (s)	(s)	signal loss	signal loss
	Acetate_1.92ppm	2.49 ± 0.15	3.20	1.42 ± 0.05	2.61	0.23 ± 0.02	0.28
	Alanine_1.48ppm	0.92 ± 0.02	1.05	0.70 ± 0.04	0.83	0.13 ± 0.01	0.12
	Ascorbate_4.52ppm	1.17 ± 0.03	1.78	0.63 ± 0.07	1.47	0.16 ± 0.01	0.14
	Choline_3.22ppm	1.27 ± 0.08	1.45	1.09 ± 0.06	1.15	0.11 ± 0.01	0.11
	Creatine_3.94ppm	0.71 ± 0.03	0.80	0.57 ± 0.06	0.69	0.15 ± 0.01	0.13
	Creatine_3.03ppm	1.34 ± 0.05	1.49	0.89 ± 0.06	1.22	0.14 ± 0.01	0.12
	GABA_1.90ppm	0.83 ± 0.06	1.17	0.40 ± 0.06	0.92	0.21 ± 0.02	0.12
	Glutamate_2.09ppm	0.68 ± 0.02	0.68	0.39 ± 0.02	0.58	0.21 ± 0.01	0.15
	Glutamine_2.14ppm	0.62 ± 0.01	0.72	0.45 ± 0.02	0.58	0.19 ± 0.01	0.15
	Glutamine_2.45ppm	0.72 ± 0.02	0.79	0.55 ± 0.03	0.64	0.16 ± 0.01	0.14
	Glycerophosphocholine_3.24ppm	0.57 ± 0.03	0.59	0.50 ± 0.02	0.46	0.18 ± 0.01	0.19
	Glycine_3.56ppm	1.49 ± 0.05	2.05	0.89 ± 0.12	1.85	0.15 ± 0.01	0.16
	Lactate_1.33ppm	0.99 ± 0.05	1.15	0.65 ± 0.06	0.77	0.15 ± 0.01	0.13
	Lactate_4.12ppm	1.68 ± 0.14	3.43	0.53 ± 0.05	0.57	0.23 ± 0.01	0.39
	Myo-Inositol_4.05ppm	0.83 ± 0.02	1.19	0.45 ± 0.02	0.92	0.19 ± 0.01	0.12
	NAA_2.02ppm	1.02 ± 0.06	1.09	0.60 ± 0.05	0.76	0.16 ± 0.01	0.13
	Phosphocholine_3.23ppm	0.62 ± 0.02	0.65	0.50 ± 0.02	0.53	0.17 ± 0.01	0.16
	Taurine_3.42ppm	1.33 ± 0.09	1.64	0.81 ± 0.07	1.34	0.14 ± 0.01	0.13
	Valine_1.04ppm	0.67 ± 0.06	0.68	0.52 ± 0.01	0.61	0.17 ± 0.00	0.14
	Scyllo-Inositol_3.34ppm	0.87 ± 0.02	1.40	0.53 ± 0.03	0.93	0.17 ± 0.01	0.14

T A B L E 2

15222594, 2019, 1, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/mrm.27696 by Cochrane France, Wiley Online Library on [16/10/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

ACKNOWLEDGMENT

We thank Dr Clément Heude for his help in analyzing the T 1 and T 2 data of biopsy specimens.

(Continues)

T A B L E 1 (Continued) 15222594, 2019, 1, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/mrm.27696 by Cochrane France, Wiley Online Library on [16/10/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

| 75 RUHLAND et AL.

(Continues)

T A B L E 1 (Continued)

SUPPORTING INFORMATION

Additional supporting information may be found online in the Supporting Information section at the end of the article.

FIGURE S1

Simplified metabolite pathways used in our laboratory


Anicet Dansou
email: anicet.dansou@insa-strasbourg.fr

Saïda Mouhoubi
email: saida.mouhoubi@insa-strasbourg.fr

Cyrille Chazallon
email: cyrille.chazallon@insa-strasbourg.fr

Anicet Suffix

Saïda Suffix

Cyrille Suffix

Optimizations of a fast multipole symmetric Galerkin boundary element method code

Keywords: SGBEM -FMM -Parallelization -Sparse matrix -Fatigue crack propagation SGBEM, FMM, Parallelization, Sparse matrix, Fatigue crack propagation 1 Introduction

Optimizations of a fast multipole symmetric Galerkin boundary element method code

U N C O R R E C T E D P R O O F Numerical Algorithms
the problem, the accuracy of the stress field results ahead of singularities, and the simplicity of the re-meshing during the propagation. An interesting approach of the BEM is the symmetric Galerkin BEM (SGBEM). The SGBEM yields symmetric coefficient matrix and presents an important gain from a numerical point of view as it reduces the computational cost and also permits BEM/FEM coupling (see references [START_REF] Costabel | Symmetric methods for the coupling of finite elements and boundary elements (invited contribution)[END_REF][START_REF] Ganguly | Symmetric coupling of multi-zone curved galerkin boundary elements with finite elements in elasticity[END_REF][START_REF] Springhetti | Weak coupling of the symmetric galerkin bem with fem for potential and elastostatic problems[END_REF]). Another advantage of the SGBEM over the traditional collocation BEM is the ability to treat hypersingular and singular integrals solely by means of standard continuous elements. For this reason, the SGBEM can easily capture the crack tip behavior and provide smoother solution in the neighborhood of geometric discontinuities. The SGBEM is used in many works for crack problems (see for example

Frangi [START_REF] Frangi | Fracture propagation in 3d by the symmetric galerkin boundary element method[END_REF] for a simple fatigue crack growth simulation; Roberts et al. [START_REF] Roberts | Sgbem modeling of fatigue crack growth in particulate composites[END_REF] and Kitey et al. [START_REF] Kitey | Modeling of crack growth through particulate clusters in brittle matrix by symmetric-galerkin boundary element method[END_REF] for crack growth in particulate composites; Xu et al. [26] for 2D crack propagation; Tavara [START_REF] Távara | Sgbem for cohesive cracks in homogeneous media[END_REF] for cohesive crack growth in homogeneous media; Nguyen et al. [START_REF] Nguyen | An isogeometric symmetric galerkin boundary element method for two-dimensional crack problems[END_REF] and Nguyen et al. [START_REF] Nguyen | Isogeometric symmetric galerkin boundary element method for three-dimensional elasticity problems[END_REF] for exploiting the advantages of isogeometric analysis).

The usual slow evaluation of double integrals in the SGBEM can be accelerated by the fast multipole method. Initially introduced by Rohklin [START_REF] Rokhlin | Rapid solution of integral equations of classical potential theory[END_REF], this algorithm considers one group of particles and represents it by an intermediate pole. As all the interactions with this group are transferred via this pole, the overall number of operations is greatly reduced. By coupling the SGBEM with the fast multipole method (FMM) and an iterative solver, the complexity of the method is significantly reduced:

O(N) for the storage requirements and O(NlogN) for the operation count [27].

Therefore, the range of boundary analysis can be extended to large-scale practical issues with good performance (see, for example, application of FMM in elastodynamics [START_REF] Chaillat | Fast Multipole Method for 3-D elastodynamic boundary integral equations[END_REF]). Trinh [25] presented a fast multipole accelerated symmetric Galerkin BEM (FM-SGBEM) code to simulate crack problems.

Over recent decades, computers have evolved a lot. Parallel architecture machines have become standard. Parallel computing techniques can significantly increase the performance of existing serial codes. Many researchers have used parallelization to accelerate the BEM (see for example [START_REF] Adelman | Fmm/gpu-accelerated boundary element method for computational magnetics and electrostatics[END_REF][START_REF] Greengard | A parallel version of the fast multipole method[END_REF][START_REF] Gu | Accelerated parallel computation of field quantities for the boundary element method applied to stress analysis using multi-core cpus, gpus and fpgas[END_REF][START_REF] Ptaszny | Parallel fast multipole boundary element method applied to computational homogenization[END_REF]). In this work, a parallel implementation of the code presented by Trinh [25] is achieved. The matrices manipulated in FM-SGBEM are almost always sparse. Several techniques can reduce the storage space of sparse matrices. The simplest method is the coordinate format (COO).

COO stores each non-zero value of the matrix with its coordinates in three vectors.

The most popular method is the compressed sparse row (CSR) format, suggested in [START_REF] Brameller | Sparsity: its practical application to systems analysis London[END_REF] and [START_REF]Sparse matrices and their applications[END_REF]. In this work, an efficient method is presented: the upper bounded incremental coordinate (UBI-COO) method. This new method takes advantage of the simplicity of the COO format and the performance of the CSR format to limit the memory required during the matrix construction phase.

This paper is organized in the following way: Section 2 introduces the SGBEM, the FMM and the initial fast multipole symmetric Galerkin BEM code for crack propagation; Section 3 presents the different optimizations: data re-using for fast matrix computation, parallel implementation, and optimized use of sparse matrix formats. Numerical examples and performances of the optimized code are presented in Section 4.

U N C O R R E C T E D P R O O F Numerical Algorithms

FM-SGBEM and initial crack propagation code

Symmetric Galerkin BEM

Let us consider a fractured solid Ω subjected to prescribed tractions t D on the boundary St and displacement constraints u D on Su. The boundary of Ω (including the crack S c ) is thus defined as S = S t S u S c . S c is conceived as a locus of displacement discontinuity (Fig. 1). The jump of the displacements can be computed Q3 as:

u(x) = u(x + ) -u(x -) (1)
where u(x + ) and u(x -) are respectively the displacement of the upper and lower faces of the crack S c = S - c S + c . The direction of the normal of the crack is by convention, pointing from S -to S + . The boundary integral formulation for this problem is written as follows:

⎧ ⎨ ⎩ B uu (u, ũ) + B tu (t, ũ) + B uu ( u, ũ) = F u ( ũ) B ut (u, t) + B tt (t, t) + B ut ( u, t) = F t ( t) B u u (u, ũ) + B t u (t, ũ) + B u u ( u, ũ) = F u ( ũ) (2)
We introduce here B tt , B tu and B u u for example the bilinear forms:

B tt (t, t) = S u S u t k (x)U k i (x, x)t i (x)dS x dS x B tu (t, ũ) = - S u S T t k (x)T k i (x, x) ũi (x)dS x dS x B u u ( u, ˜ u) = S c S c [R u] iq (x)B ikqs (r)[R ũ] ks (x)dS x dS x (3)
u, t and u are respectively the unknown on S t , S u and S c ; U k i (x, x) and T k i (x, x)

are respectively the ith displacement and traction of x due to a point load at x in Fig. 1 An elastic cracked domain

U N C O R R E C T E D P R O O F Numerical Algorithms
the direction of the kth coordinate axis. For x and x ∈ R 3 , they are called Kelvin fundamental solutions and are written as:

U k i (x, x) = 1 16πμ(1 -ν)r [r i rk (3 -4ν) + δ ik ] (4)
T k i (x, x) = - 1 8π(1 -ν)r 2 n j (x)[3r i rk rj + (1 -2ν)(δ ik rj + δ jk ri -δ ij rk )] (5) having set r = x -x, r = r , r = r/r
The formulations [START_REF] Frangi | Fracture propagation in 3d by the symmetric galerkin boundary element method[END_REF] are written in regularized form which contains only weakly singular double integrals: O r -1 . The regularization procedure involves the Stokes theorem together with the indirect regularization. The surface curl operator R arises as a result of this manipulation and is defined as:

[Ru] ks (x) = e jf s n j u k,f (x) [START_REF] Costabel | Symmetric methods for the coupling of finite elements and boundary elements (invited contribution)[END_REF] while the weakly singular fourth-order tensor B ikqs is given by:

B ikqs (r) = 1 8π(1 -ν)r [2δ qs ri rk + 2(δ ik δ qs -2νδ is δ kq -(1 -ν)δ iq δ ks )] (7)
Weak continuity requirements (C 0,α ) are enforced on u, u, ũ, ũ which is less restrictive than the collocation approach (C 1,α ). The SGBEM can therefore deal with hypersingular and other singular integral boundary equations only by means of standard continuous elements. Besides, the Galerkin discretized matrix is symmetric since the role of x or x can be exchanged and the fundamental solutions are symmetric.

Fast multipole method

The fast multipole method is an alternative technique to enhance the performance of a boundary integral analysis. The bottlenecks in BEM are due to the presence of Kernel function: the same calculation is repeated from one observation point to another thus entailing a high amount of operations. Furthermore, since the interactions between points are normally non-zero, the resulted matrix is fully populated. In FMM, intermediate points (called poles) are used to represent distant particle groups and then a local expansion is introduced to evaluate the distant contributions in the form of a series. The principle of the FMM can be illustrated in Fig. 2: We need to compute the interaction between 2 groups of points x and y (respectively, on S x and S y ). Supposing that we have n points on S x and m points on S y , we should therefore need m.n operations by conventional approach. FMM, on the other hand, uses the point O to represent S y ; the contributions from S y are thus carried out and transferred to every point x via O; the total number of operations is now reduced to only m + n which is much smaller than m.n. In conjunction with an iterative solver, FMM can generally reduce the computational complexity of a BEM problem from O(N 2 ) to O(NlogN).

In SGBEM, the FMM reformulates the kernels U k i , T k i , B ikqs into multipole series, which achieves a complete separation of the variables x and x. For this purpose, the relative solution vector r = x -x (see Fig. 3) is decomposed into x 0 , x0 . With these notations, the multipole expansion of 1/r (see [27]) is given by:

U N C O R R E C T E D P R O O F Numerical Algorithms
1 r = ∞ n=0 n m=-n (-1) n R nm (x ) n n =0 n m =-n S n+n ,m+m (r 0 )R n m (x ) R nm (y) = 1 (n + m)! P m n (cosα)e imβ ρ n S nm (y) = (n -m)!P m n (cosα)e imβ 1 ρ n+1 (8)
(ρ, α, β) being the spherical coordinates of the argument y and P m n denoting the Legendre polynomials, and with the overbar denoting complex conjugation. For given poles x 0 , x0 , the above expansion ( 18) is convergent for any x 0 , x0 such that:

x < xr 0 and x < x + r 0 (9)

Let (x 0 ) and ˜ (x 0 ) ⊂ ∂Ω denote two subsets of ∂Ω such that (18) holds for any

x ∈ (x 0 ) and x ∈ ˜ (x 0 ). Then, the contribution of surfaces (x 0 ), ˜ (x 0 ) to the bilinear form B tt (t, t), denoted by B tt (x 0 , x0 ), is given by:

B tt (x 0 , x0 ) = (x 0 ) ˜ (x 0 ) t k (x)U k i (x, x)t i (x)dS x dS x ( 10
)
Fig. 3 Decomposition of the position vector

U N C O R R E C T E D P R O O F Numerical Algorithms
and can be evaluated by replacing the kernel U k i by its multipole expansion, and likewise for the other bilinear forms. For simplicity, only the contribution B tt (x 0 , x0 ) is detailed here. The treatment of the other bilinear forms follows the same approach.

By substituting ( 18) into ( 20), the contribution B tt (x 0 , x0 ) can be written as:

B tt (x 0 , x0 ) = ∞ n=0 n m=-n (-1) n ∞ n =0 n m =-n M1 knm (x 0 )S n+n ,m+m (r 0 )M 1 kn m (x 0 ) + M1 knm (x 0 )r 0k S n+n ,m+m (r 0 )M 2 n m (x 0 ) + M2 nm (x 0 )S n+n ,m+m (r 0 )M 2 n m (x 0 ) ( 11
)
In terms of the multipole moments:

M 1 knm (x 0 ) = S u R nm (x )t k (x )dS x M 2 nm (x 0 ) = S u R nm (x )x k t k (x )dS x ( 12
)
associated to the pole x 0 and:

M1 knm (x 0 ) = S u δ ik -(3 -4ν) xk ∂ ∂ xi R nm (x )t i (x )dS x M2 nm (x 0 ) = (3 -4ν) S u ∂ ∂ xi R nm (x )t i (x )dS x ( 13
)
associated to the pole x0 . Equation ( 11) can be recast into the following equivalent form:

B tt (t, t) = ∞ n=0 n m=-n (-1) n M1 knm (x 0 )L 1 knm (x 0 ) + M2 nm (x 0 )L 2 nm (x 0 ) ( 14
)
in terms of the local expansion coefficients, related to the multipole moments by the following multipole-to-local (M2L) relation:

L 1 knm (x 0 ) = n n =0 n m =-n S n+n ,m+m (r 0 )[M 1 kn m (x 0 ) + r 0k M 2 n m (x 0 )] L 2 nm (x 0 ) = n n =0 n m =-n S n+n ,m+m (r 0 )M 2 n m (x 0 ) (15)

Octree structure

To optimize the acceleration permitted by ( 14), a hierarchical octree structure of elements is introduced. For that purpose, a cube containing the whole boundary, called level-0 cell, is divided into eight cubes (level-1 cells), each of which is divided in the same way. The cell subdivision is continued until the number of elements in a cell is

U N C O R R E C T E D P R O O F Numerical Algorithms
Fig. 4 An octree structure smaller than a given value (which is called max elem). Any given boundary element is deemed to belong to one cell of a given level only, even if is geometrically shared by two or more same-level cells (see Fig. 4).

Initial crack propagation code

The initial FM-SGBEM code and its performance are well detailed in [START_REF] Trinh | Modelling multizone and multicrack in three-dimensional elastostatic media: a fast multipole galerkin boundary element method[END_REF]25]. This Fortran code inherits a number of innovative algorithms from the BE community such as (i) the singular integration schemes by Andrä and Schnack [START_REF] Andrä | Integration of singular galerkin-type boundary element integrals for 3d elasticity problems[END_REF][START_REF] Frangi | 3d fracture analysis by the symmetric galerkin bem[END_REF], (ii) the index of severity [START_REF] Rezayat | On time-harmonic elastic-wave analysis by the boundary element method for moderate to high frequencies[END_REF], (iii) the nested Flexible GMRES which makes use of the nearinteraction matrix [START_REF] Chaillat | Fast Multipole Method for 3-D elastodynamic boundary integral equations[END_REF]; and (iv) the extension of the BIEs to multizone configurations [START_REF] Gray | Symmetric galerkin boundary integral formulation for interface and multizone problems[END_REF]. Subroutines of matrix-vector operations are taken from BLAS library. Flexible GMRES and GMRES scripts are downloaded from www.cerfacs.fr.

Figure 5 resumes the main phases of the code. In a nutshell, the aim is to solve a linear system (K.X = b) with an iterative solver: the Flexible GMRES. The matrix K which corresponds to (2) is composed of double surface integrals B tt , B tu and B u u presented in (3). B tt for example, is double integral over two surfaces Su.

The generic double surface integral takes the following aspect: The matrix K can be separated in two parts: K near and K f ar . K f ar consists of the integrals over far enough surfaces (according to ( 9)), calculated with the

I (S e , S f ) = S e S f f (x)G(x, y)g(y)dS y dS x (16) U N C O R R E C T E D P R O O F Numerical Algorithms
U N C O R R E C T E D P R O O F Numerical Algorithms
FMM (see [START_REF] Nguyen | An isogeometric symmetric galerkin boundary element method for two-dimensional crack problems[END_REF] for B tt ). K near consists of the integrals over near surfaces. Singularities will occur when these two surfaces S e and S f are coincident, adjacent by edge, or adjacent by vertex. The singular integrals are evaluated using special schemes (see [START_REF] Andrä | Integration of singular galerkin-type boundary element integrals for 3d elasticity problems[END_REF]25]). The regular integrals are evaluated with normal quadrature rule:

I (S e , S f ) = e f f (x(η))G(x(η), y(ξ ))g(y(ξ ))J y (ξ )J x (η)dξ dη ( 17
)
where

e ∈ [-1, 1] × [-1, 1] and f ∈ [-1, 1] × [-1, 1]
. This integral can be approximated by:

I (S e , S f ) N pge i=1 N pgf j =1 f (η i )G(η i , ξ j )g(ξ j )J y (ξ j )J x (η i )A j ξ j A i η i ( 18
)
where η i and A i η i denote the abscissas and weights of the Gaussian points for exterior elements; ξ j and A j ξ j denote the corresponding parameters for the interior elements.

N pge and N pgf are the number of Gaussian points for exterior and interior elements respectively.

With the separation, the linear system to solve can be written as in [START_REF] Rokhlin | Rapid solution of integral equations of classical potential theory[END_REF]. It is important to note that due to the FMM, K f ar is not stored, and the matrix-vector product K f ar * X is computed directly. When the convergence is achieved, the crack propagates, and the elastostatic code is repeated.

K near + K f ar .X = b ( 19
)
3 Improvements and optimizations

Fast computation of b and K near

In the initial algorithm, during each cycle, a layer of new elements is added to the geometry and the system (especially the matrix K near ) needs to be recomputed. If one rebuilds the coefficient matrix, the interactions between pairs of old elements will be repeated and will require wasteful operations because nothing changes in the calculation of those interactions. Therefore, starting from the second cycle, the interactions between pairs of old elements are re-used. Only the parts of the matrix that are related to the newly added elements are computed (see Algorithms 1 and 2).

Re-using the interactions between pairs of old elements for the computation of K near is a simple idea, but one must be sure to keep in the same conditions as the initial configuration. For example, the octree structure must be fixed from one cycle to another so that near elements stay near elements and the same thing goes for far to the element list of cell c.

In the initial co de, quater-point elements are used at the crack front, and they are Q5 returned back to normal QUA 8 elements before adding new quater-point elements at the new crack-front. So, those elements are modified and their contributions in the matrix are not the same after the propagation. That means their contributions can not be simply re-used. To solve this problem, the contributions of modified elements are saved format in another matrix. After each increment, the contribution of modified elements is removed from the old matrix (algorithm 2, line 11). Then, those elements are set as new elements (algorithm 2, line 9) so that their new contributions are computed (algorithm 2, line 15). The treatment of crack front elements is therefore based on [START_REF]Sparse matrices and their applications[END_REF]. It is important to note that the position of the unknowns can change from one cycle to another. Permutation operations must be performed in these cases on the matrices of the previous cycle. All of these matrix manipulations are performed while maintaining a compressed matrix format: the Compressed Sparse Row (CSR).

K i = K i-1 -K f ront i-1 + K oldf ront i + K f ront i (20) U N C O R R E C T E D P R O O F Numerical Algorithms
With the fast computation, the cost of re-constructing the coefficients matrix K near can be greatly reduced especially in the case of few number of cracks. The effect of the fast computation is shown in Table 1 for different models. In Table 1, T pre is the cumulated preparation time from cycle 2 to cycle 10. Figure 11 presents for model C8, the duration of the preparation phase for each cycle.

U N C O R R E C T E D P R O O F Numerical Algorithms

Parallel implementation with OpenMP

Parallelization is a technique for dividing a large problem into small problems that can be solved simultaneously. The aim is to solve the initial problem in the smallest possible time. A multiprocessing parallelization is achieved in this work by Q6 using OpenMP [START_REF] Chandra | Parallel Programming in OpenMP[END_REF]. OpenMP is an Application Program Interface (API) for parallel computing on shared memory architecture. It simplifies writing multi-threaded applications by using compiler directives and library routines.

The goal here is to speed up the existing code by avoiding big changes. A simple observation of Trinh's [START_REF] Trinh | Modelling multizone and multicrack in three-dimensional elastostatic media: a fast multipole galerkin boundary element method[END_REF] results (see Table 2) shows that the solving phase is time-consuming. To reduce this duration, it is necessary to reduce the number of iterations or the duration of one iteration. This work focuses on the duration of one As illustrated in Fig. 5, Sub MVP consists of the computation of matrice vector products. Code profiling shows that K f ar .w j evaluated with [START_REF] Nguyen | An isogeometric symmetric galerkin boundary element method for two-dimensional crack problems[END_REF] for B tt , is the most time-consuming part of Sub MVP due to multiple imperfectly nested loops. In multizone configurations, the number of these deeply nested loops can reach 12: This part of the code is completely reorganized and parallelized. Loop invariants are identified and moved out of loops. Vectorization work is performed for the innermost loops. For this part, the parallelization of the first loop is not interesting due to the large amount of data that would be in private. Also, the problems considered here have a maximum of three zones. Thus, the second and the third outer loops are parallelized. In the second loop, the number of iterations is the number of octree cells ncells, while in the third it is the number of elements in a cell nel. These numbers vary from one problem to another and depend on the octree construction, especially on the parameter max elem. A parameterized nested parallelism is performed with OpenMP (see Listings 1 and 2). The user can enable or disable one or both of the parallel loops and if enabled, the number of threads can be given. For small problems, the parameter max elem can have a high value (100, 200, 1000, ...). nel is thus high and ncells is small. For these cases, the third parallel loop should be activated. For large-scale problems, max elem is generally limited (50, 30, 15,...) by the computer RAM memory. nel is thus small, ncells is high, and the second parallel loop should be activated.

The speedup and the efficiency of the parallelization are shown in Table 3 for the model C216 and for one iteration. Simulations are done on a 20-core Intel Xeon E5-2630v4 processor running at 2.2 GHz. After the acceleration of the solving phase, the preparation phase is also accelerated by the parallelization of the time-consuming part of subroutine Sub b knear. The speedup and the efficiency of the parallelization are similar to those of the solving phase (Table 3). Table 4 shows the global speedup and efficiency due to the parallel implementation for the model C216. The code can be more parallelized but it will become more complex and extension work will be difficult. Although the code is not entirely parallelized, the parallelization results are very good.

Upper bounded incremental coordinate method

In large-scale simulations, the memory usage requires special attention, especially in a context of parallel computing. In the initial code, the CSR format is used to store the matrices after computation, but dense format (DNS) is used before and during the computation (see subroutine Sub b knear). Using DNS for construction causes large allocated but not used memory. Since the construction is in parallel, dense format causes pic of allocated memory. To avoid this, an upper bounded incremental coordinate method (UBI-COO) is designed for the construction of the matrices. Based on Sparsekit subroutines written by Youcef Saad [START_REF] Saad | Sparskit, a basic tool kit for sparse matrix computations[END_REF], necessary subroutines for the manipulation of the matrices in COO or CSR format are written. The coordinate format (COO) is well known for constructing sparse matrices.

A comparison between DNS and UBI-COO is presented in Fig. 12. Using DNS is simple: a dense matrix is allocated and is converted to CSR format at the end of the construction. With the UBI-COO, the coordinate format is used in an incremental way. A parameter fixes the maximal number of data in the COO. When the limit is reached (Fig. 12: step 1.6 and 2.6), the COO matrix is converted to CSR (step 1.7 and 2.7) and the CSR is cumulated with an existing CSR (step 1.8). At this step, multiple entries are also cumulated. The COO is then re-initialized for the rest of the construction. In fact, the number of non-zero is not known before the computation, so the dimension of the arrays which contain the coordinates and the non-zero values is not known. So, two parameters (p1 and p2) are used to achieve incremental COO. the COO size reaches p3, the COO matrix is converted to CSR and cumulated with a temporal CSR matrix and the COO arrays are re-set and then the multiple entries are cumulated in the temporal CSR matrix.

This upper bounded incremental coordinate method erases memory peaks.

Figure 13 presents the virtual memory needed using DNS and UBI-COO during the matrix computation for the model C8. There is no memory variation during the solution phase, so only one iteration is performed in order to focus on the preparation phase (matrix computation). It can be noticed that the maximum memory needed is greatly reduced. It can also be noticed that the duration of the construction phase is reduced (for cycles 2 and 3) because less data are manipulated.

Numerical examples

This section presents the results of all the optimizations presented in this paper. The calculation times are measured on an Intel Xeon (20 cores, 2.2 GHz) computer with 128 Go of RAM. T tot is the total time including pre-processing (input reading, octree construction, etc.) and post-processing (results writing in files). This duration is compared with Trinh's code [START_REF] Trinh | Modelling multizone and multicrack in three-dimensional elastostatic media: a fast multipole galerkin boundary element method[END_REF] duration noted T old tot . Table 5 shows computational data for static analyses. Table 6 shows computational data for propagation analyses. For the cube with 8 cracks (model C8), the evolution of the total time according to the cycle number is presented in Fig. 14. Other optimizations will be investigated in future, for example the use of hierachical matrix representation and the coupling of the FM-SGBEM code with the FEM.

The coupled code will make use of the advantages of the FEM and thus will extend

Fig. 2

 2 Fig. 2 Illustration of the FMM

Fig. 5

 5 Fig. 5 Initial FM-SGBEM code for crack propagation

3. 1 Fig. 6

 16 Fig. 6 Circular crack mesh 48 elements

Fig. 7 Fig. 8 Fig. 9 Fig. 10

 78910 Fig. 7 Circular crack mesh 768 elements

Fig. 11 Fast

 11 Fig. 11 Fast K near computation: C16 result

Listing 1 A parameterized nested parallelism: part 1 Listing 2 A

 112 parameterized nested parallelism: part 2

Fig. 12

 12 Fig. 12 DNS and UBI-COO principles

Fig. 14 Evolutions of total time 5 Conclusion

 145 Fig. 14 Evolutions of total time

  the code to inhomogeneous and anisotropic materials and non-linear constitutive behavior.

Table 1

 1 Fast K near computation: results

	No.	Model	N init dof s	N cycles	N end dof s	T pre (s)	T old pre (s)	Speedup
	1	C1	16,593	10	17,889	1051	478	2.2
	2	C4	17,754	10	22,938	2892	440	6.6
	3	C8	19,302	10	29,670	4925	857	5.7
	4	C16	22,398	10	40,830	14,387	3059	4.7

Table 2

 2 Bi-material cube with crack array: Trinh's results[START_REF] Trinh | Modelling multizone and multicrack in three-dimensional elastostatic media: a fast multipole galerkin boundary element method[END_REF] So the time-consuming parts of the routine Sub MVP are parallelized by using OpenMP directives.

	No.	N dof s	T pre (s)	N iter	T sol (s)	T tot (s)	T sol /T tot (%)
	1	401,412	5457	79	44,319	50,986	87
	2	683,148	12,197	66	79,464	95,584	83
	3	1,061,928	11,903	102	190,944	206,114	93
	iteration. Almost all the time of one iteration (more than 90%) is spent in the routine
	Sub MVP.

Table 3

 3

	Parallelization:
	efficiency (Sub MVP)	N th	T Sub MV P (s)	Speedup	Efficiency(%)
		1	221	-	-
		2	113	2.0	98
		4	6 4	3 . 5	86
		8	3 4	6 . 6	82
		12	24	9.4	78
		16	18	12.0	75
		20	17	13.3	67

Table 4

 4

	Parallelization: global
	efficiency	N th	T tot (s)	Speedup	Efficiency(%)
		1	3 771	-	-
		2	2 217	1.7	85
		4	1 241	3.0	76
		8	760	5.0	62
		12	603	6.3	52
		16	534	7.0	44
		20	528	7.1

Table 5

 5 Static tests

	No.	Model	N dof s (s)	T tot (s)	T old tot (s)	Speedup
	1	C8	19,302	169	1469	8.7
	2	C64	40,974	611	6417	10.5
	3	C1000	403,206	4478	72,107	16.1
	4	C1728	684,942	8721	119,249	13.7
	5	C1000	1,075,206	15,446	166,808	10.8
	6	C8000	3,112,206	53,288	848,162	15.9

Table 6

 6 Propagation tests: cube with crack array

	No.	Model	N init dof s	N cycles	N end dof s	T tot (s)	T old tot (s)	Speedup
	1	C8	19,302	10	29,670	889	31,396	35.3
	2	C64	40,974	8	123,918	3445	220,685	64.1
	3	C512	214,350	8	656,718	48,556	2,500,000	51.5
	4	C2744	1,078,134	3	1,868,406	399,600	8,000,000	20.1

JrnlID 11075 ArtID 781 Proof#1 -13/07/2019

Funding information This work is supported in part by the French National Research Agency (SolDuGri project ANR-14-CE22-0019) and in part by the region "Grand-Est, France."


T Landes
email: tania.landes@insa-strasbourg.fr

M Heissler
email: mathias.heissler@bas-rhin.fr

M Koehl

T Benazzi
email: theo.benazzi@insa-strasbourg.fr

T Nivola

UNCERTAINTY VISUALIZATION APPROACHES FOR 3D MODELS OF CASTLES RESTITUTED FROM ARCHEOLOGICAL KNOWLEDGE

Keywords: Restitution, Uncertainty, visualization, Archeological knowledge, Point Cloud, Assessment, reliability

In the cultural heritage field, several specialists like archaeologists, architects, geomaticians, historians, etc. are used to work together. With the upcoming technologies allowing to capture efficiently data in the field, to digitize historical documents, to collect worldwide information related to the monuments under study, the wish to summarize all the sources of data (including the knowledge of the specialists) into one 3D model is a big challenge. In order to guarantee the reliability of the proposed reconstructed 3D model, it is of crucial importance to integrate the level of uncertainty assigned to it. From a geometric point of view, uncertainty is often defined, quantified and expressed with the help of statistical measures. However, for objects reconstructed based on archaeological assumptions, statistical measures are not appropriate. This paper focuses on the decomposition of 3D models into levels of uncertainties (LoUs) and on the best way to visualize them through two case studies: the castle of Kagenfels and the Horbourg-Wihr Castellum, both located in Alsace, France. The first one is well documented through still ongoing excavations around its remains, whereas the second one disappeared under the urbanization of the city. An approach enabling, on the 3D models, not only to quantify but also to visualize uncertainties coming from archaeological assumptions is addressed. Finally, the efficiency of the approach for qualifying the proposed 3D model of the reconstructed castle regarding its reliability is demonstrated.

INTRODUCTION

The process of 3D modeling of cultural heritage buildings requires the combination of data coming from several sources. Depending on the conservation state of the building, the geometric acquisition based on surveying measurements is primordial first step. However, with the purpose of restituting the building in its original stage, it is often necessary to complete the model with supplementary material or knowledge, coming usually from archaeologists or architects working on those sites.

Restitution projects use increasingly the process of digital modeling as a means of interactive presentation and interpretation. Models are no longer limited to 2D renderings, but are now based on advanced 3D or 4D modeling techniques. This brings up new questions about the handing of the model quality. In order to evaluate the quality of a 3D model, an error budget might be elaborated. It is rather easy to assess the quality of the model from a geometric point of view, when restitution relies on measurements. However, when there is a lack of tangible data and the restitution is mainly based on archaeological knowledge or hypotheses, the way to quantify the uncertainties related to the reconstructed object remains a challenge. For some authors, a lack of certainty might even be detrimental as the models imply a kind of truth, which may hinder further interpretation and analysis (Sylaiou et al., 2009).

Concerning evaluation and representation of uncertainties in 3D models base on archaeological knowledge, papers proposed in the literature do not solve this problem universally. The methods used are various and usually adapted to the specific site under study. Depending on the project, the available data and the objectives followed, the most effective means of representation of the geometry and of the uncertainty of the 3D model must be chosen. This paper presents firstly a state of the art dealing with the quantification and the representation of uncertainties for 3D models based on archaeological knowledge. Then, applied to two case studies, an approach enabling not only to quantify but also to visualize uncertainties will be addressed.

RELATED WORK

The accuracy of the restitution process has given rise to much debate within the archaeological community (Barrat, 2016). The first difficulty comes from the fact that there is often a lack of consistent methodology for restitution, mainly because of the lack of a standard in data production. The reconstruction of partially or totally destroyed structures might rely on data from heterogeneous sources [START_REF] Rodríguez-Gonzálvez | 4D reconstruction and visualization of cultural heritage: Analyzing our legacy through time[END_REF]. Indeed, metric data acquired on site might be completed with multitude of historical documents that often have non-metric properties (historical maps, old photographs, drawings, sketches, paintings) as well as archaeological knowledge based on deductions (existence of elements based on castles constructed in the same period by probably the same architects).

Therefore the quality, accuracy and completeness of restitution depend on the way these heterogeneous data are combined. Apart from the fact that these data are non-metric, they also often lack objectivity and are usually the result of an interpretation made at a given moment. Moreover, they might not be up-to-date due to later discoveries. There is therefore of crucial importance to work in close collaboration with archaeologists and heritage architects during the whole process of 3D reconstruction.

While it is fairly easy to quantify the uncertainty of a measurement or a 3D survey, it is more difficult to qualify nonmetric data. In the engineering field, uncertainties of metric elements are expressed through statistical values (standard deviations, root mean square errors, histograms, etc.). When representing elements (in 2D or 3D), like archaeological sites which are no longer measurable (disappeared or destroyed) or which existence result from hypothesis, a solution consists in using a visualization process. Not only the chronological dating might in some cases be uncertain, but also the position of the discovered object. Obviously missing data is a major component of uncertainty [START_REF] Zuk | Visualizing Temporal Uncertainty in 3D Virtual Reconstructions, The 6th International Symposium on Virtual Reality[END_REF]. In addition, the degree of uncertainty varies considerably from one project to the other (Barrat, 2016). Any proposal relating to archaeological data therefore contains uncertainties. This uncertainty must be integrated into any visualization not only for transparency purposes but also to improve the cognitive task of spatio-temporal understanding.

For visualizing uncertainties, we might use several colors or textures, degrees of transparency [START_REF] Bertin | Semiology of Graphics: Diagrams, Networks[END_REF]. For instance, several degrees of transparency allow to emphasize the parts of the model which are assumptions (transparent) from those which existence is proven (opaque). However, experiments will show that the use of transparency makes the model sometimes difficult to understand, especially when it is of complex geometry with a lot of superimposed elements. Another solution consists in mixing several levels of details (LoD), i.e. a simplified geometry for objects which knowledge is uncertain against very high LoD and high resolution texturing for objects which level of knowledge is very high. Whatever the adopted solution, when we move from quantified uncertainty to visualized uncertainty, we obviously simplify the uncertainty to make it compatible with the visual channels, which are limited.

Whereas [START_REF] Gutierrez | Dust and light: predictive virtual archaeology[END_REF] and Freiman and Gillings (2007) focus on realism at the expense of representing uncertainties, Dell'Unto et al (2013) partially abandon visual realism in favour of the development of a colour coded uncertainty scale. The authors suggest representing the virtual hypothesis with markers indicating different levels of consistency, in order to map how the different reconstructed structures relate to the archaeological evidences. [START_REF] Kensek | Survey of methods for showing missing data, multiple alternatives, and uncertainty in reconstructions[END_REF] presents a review of the methods used in many disciplines to represent uncertainty in reconstructions, which include colouring schemes and shadings, hatching, and line types, rendering techniques or levels of transparency to distinguish easily missing and reconstructed data. Some of these methods have been successfully applied in the field of architectural and archaeological virtual reconstruction.

Although the issue of assessing and representing the uncertainty of more or less hypothetical elements composing a reconstructed 3D model is discussed in many bibliographical references, it remains impossible to solve this problem universally. However, a classification in several levels of uncertainties (LoU) seems to be realistic. This solution requires a decomposition of the uncertainty related to the element to assess, ranging from pure assumptions (high LoU) to a very high level of archaeological knowledge, i.e. the existence of the element (low LoU).

In order to assess the suggested methodology of uncertainty decomposition and visualization, our approach has been applied on two study cases: the first one is a castle which has been restituted in 3D based on existing remains and objects found during several decades of excavations; the second one is a military construction from which nothing tangible remains. In order to understand the complexity of assigning a level of uncertainty, it is important first to describe the data on which both modelling processes have been based.

STUDIED CASTELS

Kagenfels' castle

The first case of study concerns a castle located in the forest of Obernai-Bernardsviller, Alsace, France. First stones of the Kagenfels' castle have been laid in the middle of the XIII th century, around 1262. Several phases of construction gave birth to the castle in its final version, until the end of the XVI th century. In 1999, when archaeological and deforestation work began, the site of Kagenfels was like an impressive cone of scree, where only few masonry remains were visible. The castle is almost entirely ruined, but its recovered history revealed that the site has probably never been looted. Therefore, the cone of scree was of crucial importance, because the architectural elements were supposed to be present on site.

Survey of the site and the castle's remains

With the aim of reconstructing the castle in 3D, even if only a part of the structure remains, first of all, it was necessary to record the site in its current state. In this context, photogrammetric and terrestrial laserscanning campaigns have been carried out. Terrestrial laserscanning with Faro X330 scanner allowed obtaining a point cloud of 40 million points on the area. Terrestrial photogrammetric survey with a Canon EOS 700D (CMOS image sensor of 18 MP) has provided 850 photographs. Since the surrounding area gives also valuable information on the defensive system, not only the castle's remains but also about 1 square kilometer of the surrounding area has been captured through about 2000 photos by UAV (DJI Phantom 3, with CMOS image sensor of 12 MP). These data helped us to capture the upper parts of the remains, which could not be scanned from the ground. After processing of the whole datasets, a 3D model of the ruins and their surroundings has been obtained (Figure 1). The geometric accuracy of the final mesh model is about 3 cm.

Geo-localized inventory of architectural remains

Through the study of the collapses of this castle of modest size, it was conceivable to restore graphically the elevations of the disappeared facades. An exhaustive and geo-localized inventory of the architectural elements has therefore be undertaken. More than 1500 blocks have been found in the field around the remains of the Kagenfels' castle during the last 23 years of excavations. Every block in this 2500 square meters area has been recorded and located on a map (Figure 2). This work has been done in 2D on that map, but also in depth in the screes and in height compared to the level of abandonment of the castle. This three-dimensional inventory makes it possible to develop a rigorous analysis allowing to restore the initial location of each architectural element. This original method, implemented during 19 years, is one of the specificities of the Kagenfels site. It is based on the current spatial distribution of "series" of scattered blocks coming from the same architectural ruined site. [START_REF] Heissler | Rapport de synthèse, ruines du château du Kagenfels[END_REF]; each color represents a specific "series" of blocs All blocks and lapidary fragments belonging to the same architectural element, like doors, windows or arrow loops for instance constitute a "series". Each discovered block or fragment thus receives a double identification, composed of the discovery number and a serial number referring to the part of the castle it belongs to. Each series of blocks is described in detail in a specific summary sheet describing the nature of the architectural element, its characteristics, its dimensions, the three-dimensional localization of the blocks and fragments assigned to it, and obviously photographs of the elements.

The compilation of several summary sheets of series gathering the elements relating, for example to the north facade of the house, enables to conclude to the presence of a facade of four floors. This method of recording has already made it possible to offer graphic representations of the north, east and west facades of the house, as well as several towers of the castle including multiple arrow loops. This inventory gives, indirectly, information on the partitions and the internal organization of the house. As a consequence, 16 years ago the heritage architect in charge of the excavations was already able to propose a graphical representation of the ruined elevations (Figure 3). The graphical and documentary databases available today for the restitution of the castle are summarized in Figure 3. Regarding spatial organization and architectural details, the ruins of other neighbouring castles (Birkenfels, Dreistein, Ortenbourg) present many similarities with Kagenfels. They constitute local references for typologies of windows, especially the large windows with multiple rectangular fenestrons, for example. In the past, same typologies of windows were used quasi-systematically in a same geographical area.

Production of the 3D models

By cross-referencing and analyzing all data collected in the field for almost 20 years, i.e. geometrical, historical, constructive and aesthetic rules, and archaeological data, as well as the metric data acquired during the surveying campaign, a reliable 3D reconstructed model of the castle as it was in the XVIth century has been constructed in the Sketchup software (Figure 4). The whole processing chain dedicated to the modeling of the castle and its surroundings has required approximately 350 production hours. This model, firstly without texture, results from a very efficient but laborious and time consuming collaboration between geomaticians and a heritage architect who works for 25 years with volunteers at the preservation of the castle. For communication purpose, also a rendered model has been realized in the Lumion software (Figure 5) and integrated in a virtual reality (VR) system to allow immersive visits. Apart from the impressive impact, the immersive virtual visit of the restituted castle allowed to correct some parts of the model, like positions of arrow loops, which could be refined after simulating the sniper position.

Figure 5: photo-realistic rendering of the Kagenfels' 3D model [START_REF] Benazzi | Restitution 4D du Château du Kagenfels par combinaison de l'existant et d'hypothèses archéologiques pour une visite virtuelle du site[END_REF] The way the levels of uncertainty have been integrated to the process of restitution is explained in section 4.

Horbourg-Wihr's castellum

A What is known today about the geometry of this castellum is listed in [START_REF] Fuchs | Le castellum : une étude lacunaire. Actes d'ARCHIHW 2-1996 : Horbourg-Wihr à la lumière de l'archéologie[END_REF]:

-The building forms a quasi-regular quadrangle and its sizes are of 170 meter in the East, 172 meter in the South, 166 meter on the West and 172 meter in the North, by measuring from the external face of the rampart; and its surface is 2,87 hectares.

-Four round towers are located in the angles of the castellum.

-Two simple entrances are present in the middle of the north and south sides of the rampart and two double entrances are located in the middle of east and west sides.

-Two rectangular towers are placed on both sides of every entrance of the castellum for a total of eight towers on the whole perimeter of the rampart.

-Eight intermediate half-round towers complete the defense system and are positioned between entrances and angle towers.

-Nothing is known about the internal arrangement of the camp except the existence of a temple having preceded the construction of the castellum.

In summary, the available data are plans and sketches realized by the architect C. Winkler, the reports and photos stemming from archaeological excavations, a drawing from a restitution work of P.Y. Videlier (2016), point clouds from the acquisition of a few retrieved blocks, a publication describing the castellum, some information about its size by [START_REF] Fuchs | Le castellum : une étude lacunaire. Actes d'ARCHIHW 2-1996 : Horbourg-Wihr à la lumière de l'archéologie[END_REF] and finally the theoretical plan of E.A. Herrenschneider (1894).

The planimetric sizes of the angle-towers, intermediate towers and entrance-towers could be determined from archaeological excavations and from theoretical plan drawn up from these.

Concerning the temple, the archaeological excavations of 1884/85, then of 2004 allowed to bring up-to-date its foundations. A capital of the Tuscan order was found by E.A. Herrenschneider in 1884 and elements of ledge and molding were also found in re-use. -The thickness of walls and a part of the rises of certain towers were preserved on the camp of Boppard, which is of the same time-period as that of Horbourg-Wihr. That's why the information to restore the thickness of the walls of towers has been taken from this site and that of Yverdon-les-Bains.

Other restitutions of castella like the restitution of the castrum of Yverdon-les-Bains, or of the museum of Bad Kreuznach or Kaiseraugst allowed to propose a model of the internal organization of the camps, for which no information was highlighted during archaeological excavations with the exception of that of the temple.

Production of the 3D models

The virtual model of the castellum of Horbourg-Wihr has been realized based on the various types of data presented previously and from architectural hypotheses. Indeed, the available data being incomplete, it is not possible to restore buildings in detail and in a homogeneous way. To model these buildings in 3D, it was thus essential to lean on the knowledge and the hypotheses of the archaeologists. It is important to mention that the first aim of the restitution project was to highlight this historic disappeared site of the municipality of Horbourg-Wihr. The main objective was to make discover the site to the visitors and to the inhabitants. The 3D modelling allowed ending in a 3D model proposal, developed gradually in dialogue with archaeologists who studied and worked more than 20 years on the concerned building (Figures 7 and8). The objective was not only to propose a depiction in photorealistic images of the building, but also to set up interactive solutions of accessible valuation for the public. For that purpose, the solution chosen to discover the site was to realize a virtual visit, allowing the user to see the building since nine points of view in the city. This virtual visit allows to discover the site in an interactive way and to learn more on its history thanks to the integration of texts, images and videos (Panotour software). The 3D model has been put online, in order to allow the public to explore the castellum interactively. Figure 9 shows the results on the Sketchfab plateform, which gives also the possibility to visualize the model in VR in the scale 1 to 1.

Figure 9: 3D model of the Calstellum on Sketchfab platform [START_REF] Nivola | Modélisation 3D du Castellum de Horbourg-Wihr et exploitation de la réalité augmentée pour une mise en valeur dans la trame urbaine contemporaine[END_REF] In parallel to the 3D model productions of both castles described in this section, another stake of both projects was to succeed in modelling the building although archaeological knowledge and data were sometimes incomplete. Thus an approach allowing reporting uncertainties of restitution, bound to the reliability of the data has been developed and is presented in next section. It allows to be transparent regarding the modelling method and to keep a "traceability" of the choices made during the restitution.

UNCERTAINTY VISUALISATION APPROACH

Since both castles do not benefit from the same amount of available data, the way the uncertainties have been managed in each model is described for each one, separately. However the same visualization approach has been adopted, i.e. a representation in color gradients from red to green.

Kagenfels' castle

In order to evaluate the degree of uncertainty of the reconstructed castle, firstly the 3D model had to be structured in characteristic geometrical and semantical parts. To facilitate the future update of the model, each object has been modeled independently and in a structured way in a common repository.

Then we confronted the nature of the elements composing the castle (wall, roof, frame, etc.) with the level of knowledge associated to this element. This correspondence is a very challenging task.

In a first stage, for every recorded element (door, window frames, etc.), we planned to represent two uncertainties criteria simultaneously: a) uncertainty related to the geometry of the element (dimensions, shape); b) uncertainty related to their positioning on the structure.

Obviously, only one visual channel is available to represent these two informations. That's why a first experiment consisted in introducing a radiometric gradient (from red to green) to represent the geometric uncertainty and to use a degree of opacity (from transparency to opacity) to represent the positioning uncertainty. As expected, the resulting classes were not trivial to distinguish on the display. After investigations, we concluded that both variables are correlated. A well-known object results from the assembly of several fragments found on a specific location. Thanks to the "series", the location of the fragments related to the collapse sequence and given by the stratigraphic analysis allows to deduce the original location of the object on the castle. Therefore, we decided to merge both variables and to define five LoUs, ranging from simple supposition (highest LoU) to a very high level of knowledge (lowest LoU). Finally, a representation highlighting these five levels decomposed in a color gradient from red to green is proposed.

In this context of classifying the elements of the castle according to the level of knowledge assigned to them, the collaboration with the archaeologist in charge of the project was essential.

The work of assessing the five LoU of more than 300 elements composing the 3D model has been reported in an assessment table, from which Table 1 presents one line. For the specific case of a wall, Table 1 shows the way the knowledge associated to a wall might be considered as a "pure assumption" (red), as "weak" (orange), as "moderate" (yellow), as "high" (green) and as "very high" (dark green).

In order to apply the result of that assessment work on the 3D model, it has been decided to duplicate the final model depicted in Figure 5, to remove all the textures assigned and to apply one of the five colors considered to represent the uncertainty of the element (Figure 10).

Also the elements constructed based on cross-referencing to other castles of the same period had to be affected to one of these classes. For example, figure 10 shows that all timber and gables are assigned to LoU5, whereas the castle's walls remain and belong therefore to LoU1. Table 3. Links between LoDs and LoUs for sub-components of the double entrance.

DISCUSSION

The challenges of these projects, from an archaeological point of view, were to achieve the best possible restitutions of buildings, although knowledge and data were sometimes incomplete or missing. Both final rendered deliverables look impressive, but the degree of knowledge related to them is not comparable. Whereas the whole Castellum modeling was based mainly on assumptions, only some minor parts of the Kagenfels model are the product of subjective assumptions. Whereas for the castellum 4 levels of large classes were sufficient, for the Kagenfels the interval of uncertainty has been splitted into 5 levels of already well known objects.

In fact, although both projects used the decomposition in levels of uncertainties, the scales of knowledge are not the same for both castles. For the Kagenfels, the functional logic was in the center of the uncertainty assignment. Only components of the castle that must have existed have been represented. The highest level of uncertainty in the Kagenfels' model (LoU 5, qualified as "pure assumption") concerns the shape of the object but not its existence (wall, roof, etc.). The highest level of uncertainty in the Castellum's model (LoU4, qualified as "pure assumption") concerns the existence of the elements, like the buildings of the military camp. Indeed, we do not know anything about the internal spatial organization of the military camp.

If the decomposition in levels of uncertainties remains site dependent, the uncertainty visualization in a gradient of colored LoUs is an efficient solution and has been validated. With these visualizations, the 3D models can be considered as an archaeological synthesis tool allowing to identify directly the reliability of the reconstructed parts. This approach can be applied to every castle or historical monument under the condition that the table describing the level of knowledge (or LoUs) is joined to the model.

CONCLUSION AND FUTURE WORKS

The purpose of this paper was to propose a way to show clearly the reliability of the elements composing a 3D model obtained based on excavations and archaeological assumptions. Our approach accounts for the uncertainties of restitution related to data reliability. A hierarchical classification allows data to be organized according to their LoU. It considers levels of uncertainties ranging from from pure assumptions (highest LoU) to a very high level of knowledge (lowest LoU).

The longing to be transparent on the created model and to maintain a traceability of reasoning that guided restitution are two essential factors in an archaeological restitution process. Moreover, the efforts put to produce these representations for two castles paid off, because the resulting models provide a solid foundation for upcoming updates, in case of future on-site discoveries Indeed, it is important to remember that digital models produced are not definitive and that new excavations might lead archaeologists to review some hypotheses concerning the architecture of buildings already modelled.

Figure 1 :

 1 Figure 1: 3D model as textured mesh model of the Kagenfels' remains in its current stage (end 2017), on about 80 m x 60 m.

Figure 2 :

 2 Figure 2: Location plan of inventoried blocks (Mathias Heissler, 2007); each color represents a specific "series" of blocs

Figure 2 :

 2 Figure 2: Proposition of graphical restitution of the castle in 1550 (by Mathias Heissler, 2002)

Figure 3 :

 3 Figure 3: Collected data as a support the modeling (a) iconographic documents; (b) and (c) photographs; (d) digitalization on orthophotos; (e) photographs of similar sites or roof structures; (f) and graphical restitution propositions

Figure 4 :

 4 Figure 4: 3D reconstructed model suggested for the Kagenfels castle, as it was in 1550.

  second case of study concerns a castellum of the IV th century situated in Horbourg-Wihr, Alsace, France. It disappeared under the urbanization of the city of Horbourg-Wihr but is known thanks to numerous archaeological works led for the 19th century, with important excavations realized by E.A Herrenschneider and architect C. Winkler during years 1884-85. Their investigations allowed to draw up a plan of the castellum. Many more excavations were since realized and others are still programmed to try to learn more about the history of the city. Figure 5-(a) represents the theoretical 2D plan of the castellum of Horbourg-Wihr drawn up by E.A Herrenschneider in 1894. The same excavations also concerned the 16th century castle of Horbourg-Wihr, which appears in the northeast corner of the castellum Figure 6-(b) is a graphical representation of castellum proposed by architect C. Winkler in 1905.

Figure 6 .

 6 Figure 6. Graphical data of the castellum; (a) Theoretical plan of the castellum and the castle of Horbourg-Wihr, prepared by E.A. Herrenschneider (1894); (b) Horbourg-Wihr's castellum elevation, according to architect C. Winkler (1905).

Figure 6 :

 6 Figure 6: Schematic of the castellum and 16 th Century castle (topographic map dating from 1900 repositioned on the cadastral background of Horbourg-Wihr)

Figure 7 :

 7 Figure 7: 3D model of the castellum without texture.

Figure 8 :

 8 Figure 8: Castellum main component library

Figure 10 :

 10 Figure 10: Uncertainty visualization based on five levels of uncertainty (LoUs), for the castle of Kagenfels, Alsace (southeast point of view)

Figure 11 :

 11 Figure 11: Uncertainty visualization based on four levels of uncertainty (LoUs). Case of the castellum of Horbourg-Wihr, Alsace.To deepen the comprehension of the 3D model, several sheets have been completed for each of the main components of the model: internal buildings, double entrances, single entrances, rampart, temple, corner towers, intermediate towers and for the entrance towers. These illustrated sheets include a description of the component, an exploded view of the component, the presentation of the data used to model the component and its origin, a bursting view of the component and the various sub-components, the main dimensions of the building, the interpretation that was made by the archaeologists, and finally a coloured view of the component according to the four LoUs in relation to the LoD of each component. Figure12shows a proposal for the modelling of the double entrance door. Table3illustrates the links between LoU and LoD for each sub-component.

Table 1 :

 1 Description of the level of knowledge (and level of uncertainty) which might be associated to a wall belonging to a floor.

			Level of knowledge
		Pure	Weak	Moderate	High	Very high
		assumption	knowledge	knowledge	knowledge	knowledge
		(LoU5)	(LoU4)	(LoU3)	(LoU2)	(LoU1)
		No remain,	Remain of	Remain of	Proven	Proven
		nor	low height	more than 1	existence	existence
		iconography	(< 1m). But	m. But no	with	with known
		or	no data	data about	uncertaintie	height (or
	Wall	document.	about initial	initial	s : uncertain	assessable)
		But	height nor	height, nor	height but	and
		supposed to	shape or	shape or	shape /	shape/functi
		exist for	function of	function of	function	on known
		defensive,	the wall	the wall	known
		structural or
		formal
		necessities

Table 2 .

 2 Description of the level of uncertainty which might be associated to the 3D model of the Castellum.

At this stage, the resulting 3D digital model can be used as an archaeological synthesis tool to identify directly what is known and not (Figure

11

).

The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLII-2/W9, 2019 8th Intl. Workshop 3D-ARCH "3D Virtual Reconstruction and Visualization of ComplexArchitectures", 6-8 February 2019, Bergamo, Italy This contribution has been peer-reviewed. https://doi.org/10.5194/isprs-archives-XLII-2-W9-409-2019 | © Authors 2019. CC BY 4.0 License.

This contribution has been peer-reviewed. https://doi.org/10.5194/isprs-archives-XLII-2-W9-409-2019 | © Authors 2019. CC BY 4.0 License.