Clément Jailin
email: clement.jailin@ge.com

Stéphane Roux

David Sarrut

Simon Rit

Projection-based dynamic tomography

Keywords: Computed tomography, Motion compensation, Tomographic reconstruction, Breathing motion

come

Introduction

Being able to image faithfully and precisely a body or a specimen is of utmost importance for medical diagnoses or for non-destructive evaluation. Among many different 3D imaging techniques (Magnetic Resonance Imaging (MRI), Ultrasound imaging, ...), Xray computed tomography (CT) is frequently used, because of, e.g., its spatial resolution, its rapidity, and its contrast between soft and hard biological tissues. The actual limitation is essentially set by the biological damage due to x-ray radiation.

Initially developed for medical imaging [START_REF] Kalender | X-ray computed tomography[END_REF][START_REF] Wang | An outlook on X-ray CT research and development[END_REF], tomography is now widely used in a broad range of applications ranging from biology [START_REF] Fogelqvist | Laboratory cryo X-ray microscopy for 3D cell imaging[END_REF] to material science [START_REF] Baruchel | X-Ray Tomography in Material Sciences[END_REF][START_REF] Maire | Quantitative X-ray tomography[END_REF][START_REF] Salvo | 3D imaging in material science: Application of X-ray tomography[END_REF]. The reconstruction of the 3D volume from the sinogram (collection of all acquired projections) is an inverse problem [START_REF] Kak | Principles of computerized tomographic imaging[END_REF]. The technique relies on the strict satisfaction of conditions, in particular concerning the geometry of the set-up and the motion of the sample as a perfectly controlled rigid rotation with prescribed axis and angles. Displacement in tomography can be either fought to avoid motion blur in the image or identified to feed models. If the motion happens during the imaging process, it leads to a poor quality reconstruction with blurry edges [START_REF] Milanfar | A model of the effect of image motion in the radon transform domain[END_REF]. This displacement may come from various reasons: motion of the sample itself (e.g., due to the experiment itself, heartbeat and breathing [START_REF] Vandemeulebroucke | Spatiotemporal motion estimation for respiratory-correlated imaging of the lungs[END_REF], creep, etc.), unexpected motion of the setup, uncalibrated geometry or related to quick phase and/or texture changes [START_REF] Babin | Fast X-ray tomography analysis of bubble growth and foam setting during breadmaking[END_REF][START_REF] Limodin | In situ and real-time 3-D microtomography investigation of dendritic solidification in an Al-10 wt.% Cu alloy[END_REF]. Different solutions have been found to reduce the scanning motion: using smart sampling in time for periodic motions (e.g., based on phase measurements [START_REF] Vandemeulebroucke | The POPI-model, a point-validated pixelbased breathing thorax model[END_REF][START_REF] Mokso | Four-dimensional in vivo X-ray microscopy with projection-guided gating[END_REF]), accelerating the scanning procedure (e.g., faster setup, brighter sources [START_REF] Michiel | Fast microtomography using high energy synchrotron radiation[END_REF], multiple sources [START_REF] Wang | A scheme for multisource interior tomography[END_REF][START_REF] Flohr | First performance evaluation of a dual-source CT (DSCT) system[END_REF], developing ultrafast scanning procedure in the synchrotron [START_REF] Maire | 20 Hz X-ray tomography during an in situ tensile test[END_REF][START_REF] Dobson | 4-D imaging of sub-second dynamics in pore-scale processes using real-time synchrotron X-ray tomography[END_REF]), using, in medical imaging, sedative, beta-blocker to reduce the heartbeat rate (especially for pediatric use [START_REF] Rashed | Effects of dose reduction on diagnostic image quality of coronary computed tomography angiography in children using a third-generation dual-source computed tomography scanner[END_REF]), requesting the patient not to move, or to stop breathing, etc. On the contrary, if this motion is measured and controlled, it gives access to extremely rich quantitative data that are precious pieces of information for diagnosis [START_REF] Yue Pan | Assessment of lung biomechanics in copd using image registration[END_REF], model identification and data assimilation (e.g., in material science [START_REF] Maire | Quantitative X-ray tomography[END_REF][START_REF]Full-field measurements and identification in solid mechanics[END_REF][START_REF] Buljac | Digital volume correlation: Review of progress and challenges[END_REF], bio-mechanics [START_REF] Moireau | Sequential identification of boundary support parameters in a fluid-structure vascular model using patient image data[END_REF][START_REF] Han | A hybrid patientspecific biomechanical model based image registration method for the motion estimation of lungs[END_REF][START_REF] Argani | Deformation and dynamic response of abdominal aortic aneurysm sealing[END_REF]). From volumes imaged without motion (from static or 4D acquisitions), spatio-temporal deformable registration has received considerable attention in the literature, in material science [START_REF] Hild | Three-dimensional analysis of a tensile test on a propellant with digital volume correlation[END_REF][START_REF] Buljac | Digital volume correlation: Review of progress and challenges[END_REF] and in medical imaging, in cardiac image analysis [START_REF] Schaerer | A dynamic elastic model for segmentation and tracking of the heart in MR image sequences[END_REF][START_REF] Peyrat | Registration of 4D cardiac CT sequences under trajectory constraints with multichannel diffeomorphic demons[END_REF], but also for respiratory-correlated imaging of the thorax [START_REF] Vandemeulebroucke | Automated segmentation of a motion mask to preserve sliding motion in deformable registration of thoracic CT[END_REF][START_REF] Vandemeulebroucke | Spatiotemporal motion estimation for respiratory-correlated imaging of the lungs[END_REF][START_REF] Zhang | A biomechanical modeling guided CBCT estimation technique[END_REF]. When the sample is imaged with motion, or with an uncalibrated tomograph, different methods have been developed in the literature to get a satisfactory reconstruction quality. They are often composed of two successive procedures: (1) evaluation of the motion (or calibration parameters), based on either additional sensors, or full field measurement methods on the acquired images and (2) correction of the reconstruction procedure using the estimated motion. After having measured the displacement field, the reconstruction algorithm can be enriched with motion. Techniques using a motion-compensated Filtered Backprojection (FBP) algorithms have been developed [START_REF] Prummer | Cardiac C-arm CT: Efficient motion correction for 4d-fbp[END_REF][START_REF] Prummer | Cardiac C-arm CT: a unified framework for motion estimation and dynamic CT[END_REF][START_REF] Rit | Comparison of analytic and algebraic methods for motioncompensated cone-beam CT reconstruction of the thorax[END_REF] as well as iterative algorithms [START_REF] Wang | Simultaneous motion estimation and image reconstruction (SMEIR) for 4D cone-beam CT[END_REF][START_REF] Dang | A pilot evaluation of a 4-dimensional cone-beam computed tomographic scheme based on simultaneous motion estimation and image reconstruction[END_REF][START_REF] Guo | Reconstruction of a high-quality volumetric image and a respiratory motion model from patient CBCT projections[END_REF]. Each (filtered) back-projection is warped with the measured spatial displacement field (estimated at the corresponding acquisition time) so that the sum of all those warped back-projections reconstructs the volume in a static configuration. Instead of developing techniques based on measurement between multiple volumes (e.g., in [START_REF] Ruhlandt | Four dimensional material movies: High speed phase-contrast tomography by backprojection along dynamically curved paths[END_REF][START_REF] Schryver | Motion compensated micro-CT reconstruction for in-situ analysis of dynamic processes[END_REF][START_REF] Zang | Space-time tomography for continuously deforming objects[END_REF][START_REF] Zang | Warp-and-project tomography for rapidly deforming objects[END_REF][START_REF] Mira | 3d hermite transform optical flow estimation in left ventricle ct sequences[END_REF]), the projection-based strategy is very well fitted to the displacement measurement as it provides a spatial displacement field for each frame [START_REF] Guillard | Dynamic X-ray radiography reveals particle size and shape orientation fields during granular flow[END_REF]. Developed in the medical imaging field and called 2D/3D image registration [START_REF] Wang | Simultaneous motion estimation and image reconstruction (SMEIR) for 4D cone-beam CT[END_REF][START_REF] Zeng | Respiratory motion estimation from slowly rotating X-ray projections: Theory and simulation[END_REF][START_REF] Zeng | Estimating 3-D respiratory motion from orbiting views by tomographic image registration[END_REF][START_REF] Hinkle | 4D CT image reconstruction with diffeomorphic motion model[END_REF][START_REF] Markelj | A review of 3D/2D registration methods for image-guided interventions[END_REF][START_REF] Delmon | In-room breathing motion estimation from limited projection views using a sliding deformation model[END_REF][START_REF] Suzuki | Projection-based motion estimation for cardiac functional analysis with high temporal resolution: a proof-of-concept study with digital phantom experiment[END_REF], a similar approach has also been developed in material science and called Projection-based Digital Volume Correlation, P-DVC [START_REF] Leclerc | Projection savings in CT-based digital volume correlation[END_REF][START_REF] Taillandier-Thomas | A soft route toward 4D tomography[END_REF][START_REF] Jailin | In situ µCT-scan Mechanical Tests: Fast 4D Mechanical Identification[END_REF][START_REF] Jailin | Fast tracking of fluid invasion using time-resolved neutron tomography[END_REF]. In those approaches, an initially well reconstructed volume is warped in space and time such that its projections match the recorded deformed projections. However, knowing a reference static volume is not accessible in all medical acquisitions.

It is proposed to combine both the projection-based measurement and the motioncompensated algebraic reconstruction to achieve the desired dynamic tomography. An initially motion-blurred reconstructed volume is used as the reference volume and leads to an evaluation of the displacement field. The measured motion is then accounted for in a motion-compensated ART algorithm and these two steps are performed iteratively. In [START_REF] Jailin | Self-calibration for lab-µCT using space-time regularized projection-based dvc and model reduction[END_REF][START_REF] Van Leeuwen | Automatic alignment for three-dimensional tomographic reconstruction[END_REF][START_REF] Zhao | Motion compensated dynamic mri reconstruction with local affine optical flow estimation[END_REF], an affine transformation model was used to calibrate the acquisition system or to reconstruct moving samples. In [START_REF] Burger | A variational reconstruction method for undersampled dynamic x-ray tomography based on physical motion models[END_REF][START_REF] Jailin | Dynamic tomographic reconstruction of deforming volumes[END_REF], the coupled measurement/compensated-reconstruction framework was developed with 2D synthetic examples of moving phantom. Those synthetic test cases resulted in high quality reconstructions as well as the full space-time kinematics to be identified. In [START_REF] Wang | Simultaneous motion estimation and image reconstruction (SMEIR) for 4D cone-beam CT[END_REF][START_REF] Dang | Simultaneous 4D-CBCT reconstruction with sliding motion constraint[END_REF], the simultaneous motion estimation and image reconstruction (SMEIR) is developed. After grouping the projections by phase, multiple 'approximately still' reconstructions of the volume (ART algorithm regularized with total variation) are used to measure and initialize the major part of the 4D motion. This motion is then precisely updated from projection-based approach (regularized with a 3D mesh and a bio-mechanical model in [START_REF] Zhang | A biomechanical modeling guided CBCT estimation technique[END_REF][START_REF] Zhong | 4d cone-beam ct reconstruction using multi-organ meshes for sliding motion modeling[END_REF][START_REF] Zhou | General simultaneous motion estimation and image reconstruction (g-smeir)[END_REF]). In the study we propose, the (quasi-)periodic aspect of the motion was not exploited (to potentially take into account all possible form of displacements) and no initial 'static' volumes could be reconstructed. We propose in the present paper a dynamic tomography method based on the consistency of the projection. A fully coupled projection-based measurement method with a motion-compensated ART algorithm embedded in a multi-scale framework is developed. The technique is applied on the numerical simulation of breathing from static chest CT and the cone beam (CB-) CT of a breathing patient acquired for image guidance of radiotherapy at the Léon Bérard Cancer center (Lyon, France). The proposed method simultaneously provides, with nothing more than a standard acquisition, the accurate reconstruction of the lungs and their 4D regularized motion along the approximately 35 irregular respiratory cycles taking place during the acquisition. A quality criterion based on a projection residual error allows the volume improvement and the measured displacement field to be quantified.

Coupled reconstruction-measurement problem

Acquisition with motion

The collected intensity for each detector position r = [r, z] (where z is parallel to the specimen rotation axis, and r is perpendicular to it) and rotation angle, θ, is to be related to the relative beam intensity attenuation. The Beer-Lambert law relates the line integral of the material absorption along the X-ray path L(r, t) from source to detector at position r, to the recorded intensity I(r, t) at time t (and rotation angle θ(t)). In the present work, the projection operator Π is computed using the ASTRA toolbox as defined in [START_REF] Van Aarle | The ASTRA Toolbox: A platform for advanced algorithm development in electron tomography[END_REF] and, more precisely, its GPU implementation,

p(r, t) ≡ -log(I(r, t)/I 0 (r)) = L(r,t) g(x, t) dx ≡ Π θ [g(x, t)] (1)
with g(x, t) the linear attenuation coefficient, referred to as the volume in the following and I 0 (r, t), the intensity recorded without sample along the beam. Possibly, typical tomography artifacts [START_REF] Boas | CT artifacts: causes and reduction techniques[END_REF][START_REF] Vidal | Investigation of artefact sources in synchrotron microtomography via virtual X-ray imaging[END_REF] (e.g., beam hardening, intensity variations, etc) can be corrected at this stage. The collection of N θ projections p(•, t) for each angle θ(t) is called the sinogram.

With the acquisition noise η, the projection can now be written p(r, t) = Π θ [g(x, t)] + η(x, t). When motion v(x, t) is considered during the experiment, the projection of the deformed volume g(x, t), defined for every voxel x = [x, y, z] is written

p(r, t) = Π θ [g(x + v(x, t))] + η(x, t).
(

) 2
For given estimated reconstruction f and displacement field u, consistency of the model [START_REF] Van Leeuwen | Automatic alignment for three-dimensional tomographic reconstruction[END_REF] can be evaluated from the so-called projection residuals, ρ(r, t; f, u)

ρ(r, t; f, u) = Π θ [f (x + u(x, t))] -p(r, t). (3)
Reconstructed volume and motion estimate are both determined from the minimization of a functional Γ(f, u) which is chosen as the quadratic norm of these residuals (that can be shown to be optimal for a white Gaussian noise affecting p). The mean square error minimization is hence written

Γ(f, u) = r,t (ρ(r, t; f, u)) 2 . (4)
This functional also provides a validation metric of the procedure (e.g., quality of the reconstruction and estimated displacement field). The general coupled problem is solved as a fixed point problem, considering two staggered steps iteratively, as sketched in Figure 1, such that

f = Argmin f Γ( f , u), (5)
u = Argmin ū Γ(f, ū). (6)

Motion-compensated reconstruction

The first part of the two step procedure consists in reconstructing the 3D structure of the sample from the acquired sinogram and the estimated displacement field (initialized at 0 unless additional sensors information or initial guess can be exploited). The use of an ART reconstruction algorithm is natural as it reduces to equation [START_REF] Maire | Quantitative X-ray tomography[END_REF] in the absence of motion u = 0. The SART algorithm [START_REF] Anders | Simultaneous algebraic reconstruction technique (sart): a superior implementation of the art algorithm[END_REF][START_REF] Kak | Principles of computerized tomographic imaging[END_REF] was chosen here. At each iteration, indexed by (n), the volume is warped with the estimated displacement field and its projection is compared with the recorded one, and their difference defines the projection residual field (see equation [START_REF] Fogelqvist | Laboratory cryo X-ray microscopy for 3D cell imaging[END_REF]). The latter is normalized and back-projected ∆f (n)

(x, t) = B θ [ρ (n) (r, t)]
with the backprojection operator B θ . Finally, the correction term is unwarped to the frame of the undeformed state, so that it matches the reference configuration and it is added to the volume, written in a reference configuration, at t = 0 with u(x, 0) = 0:

f (n) (x, 0) = f (n-1) (x, 0) + ∆f (n) (x -u(x, t), t).
The volume warping with the u(x, t) field is performed with tri-linear interpolation. A convergence criterion has to be chosen as in the ART procedure. However, the value of the cost functional Γ cannot be used directly in this case as the reconstruction is imperfect. A stopping criterion based on the incremental variation of Γ (or an escape condition limiting the maximum number of iterations) can be set. Because of the difficulty to choose a stopping criterion on the functional evolution and because minimizing the computation time was not a target, it was decided in this paper to limit the number of iterations. This limit N DynART was set to be large enough to allow convergence of the reconstruction. The algorithm of the motion-compensated reconstruction is shown in the Supplementary materials.

While several other algorithms are frequently used to improve the quality of the reconstructed volumes (e.g., , total variation), we did not use them to observe the benefit of our motion correction alone.

As additional metrics to evaluate the sharpness improvement of the reconstructed volumes, we computed (1) the L2 norm of the gradients:

1/N x • Ω |∇f (x)| 2 and (2)
the 10%-90% rise response of the image edges [START_REF] Higaki | Visualization of simulated small vessels on computed tomography using a model-based iterative reconstruction technique[END_REF] (distance between 10% and 90% of the edge intensity).

Space-time measurements

Constrained kinematics: The second part of the two-step procedure is the projectionbased displacement field measurement following a P-DVC strategy. This procedure uses as the reference volume the previously reconstructed volume and minimizes, now with respect to the displacement field, the same cost function Γ, i.e., the RMS of the projection residuals. DVC algorithms are known to be ill-posed [START_REF]Full-field measurements and identification in solid mechanics[END_REF], consequently the displacement field has to be regularized. Among many possible choices to regularize the displacement field, our choice is to write the displacement field on a reduced kinematic basis composed of a limited number of degrees of freedom

u(x, t) = ij a ij φ i (x)σ j (t) (7)
with a ij unknown amplitudes associated with φ i (x) and σ j (t) respectively space and time functions that characterize the space-time evolution. The displacement space can be expressed as the collection of independent subsets (called local DVC [START_REF] Smith | Digital volume correlation including rotational degrees of freedom during minimization[END_REF]), or from a global approach [START_REF] Besnard | finite-element" displacement fields analysis from digital images: application to portevin-le châtelier bands[END_REF], such as using a finite-element mesh [START_REF] Zhong | 4d cone-beam ct reconstruction using multi-organ meshes for sliding motion modeling[END_REF] or b-spline interpolations, etc.

The time parameter t used in the σ j function can be, for example, substituted with the signal of additional sensors measuring a force (as used in most mechanical experiments) or a phase signal (accessible in four-dimensional cardiac and thoracic imaging). Finally, both space and time variations can be coupled using a model-based kinematics (called "Integrated Digital Image Correlation" [START_REF] Mathieu | Estimation of elastoplastic parameters via weighted FEMU and integrated-DIC[END_REF]) or using knowledge from prior similar experiments then reduced with singular value decomposition [START_REF] Charbal | Spatiotemporal regularization for digital image correlation: Application to infrared camera frames[END_REF]. In [START_REF] J R Mcclelland | A generalized framework unifying image registration and respiratory motion models and incorporating image reconstruction, for partial image data or full images[END_REF], the authors show that the identification of a regularized breathing motion from different 3D dynamic volumes leads to better results allowing to work with partial data.

The choice of both space and time functions φ i (x) and σ j (t) and their numbers are very important. Too many degrees of freedom leads to ill-conditioning and too few introduces model errors.

Relaxed kinematics: Because the imposed time evolutions (5) may not be perfect, it can be updated without constraints. The discrete time is written t k where k labels the different instants where a radiograph is acquired. The updated discrete temporal evolution n is hence written as

σ n (t k ) = b n k (8)
with b n k time unknowns identified with the space unknowns from a fixed point procedure and initialized with the breathing signal. It is noteworthy that this displacement mode should converge to the one that maximizes the residual field correction (although it is highly influenced by the signal initialization). After having identified a first space-time couple, the breathing time signal is re-identified from the estimated breathing space motion, time-wise, at the next iteration. This allows slight amplitude corrections in the time evolution.

Resolution: The non-linear problem is solved using a Gauss-Newton's algorithm, i.e., successive linearizations around the current solution (detailed equations are presented in the Supplementary material), until a stopping criterion is reached. In the following application, a single linearization was performed before updating the volume. All components of the gradient of the cost function with respect to all degrees of freedom are computed using first order centered finite differences. Because of a small number of degrees of freedom, each linearized problem can be easily inverted (this is to be contrasted with other convergence procedures based on multiple search directions on separated spaces [START_REF] Passieux | A space-time PGD-DIC algorithm[END_REF]).

The residual field of this procedure ρ, as defined previously, informs on the quality of the procedure and chosen model. From this residual field, both model and procedure can be (un)validated and, if required, the model can be enriched with additional degrees of freedom.

Details of the P-DVC algorithm are given in the Supplementary information.

Multi-scale procedure

The key to be able to solve the previous two step framework, even with large displacements, is to apply a multi-scale pyramidal coarse-to-fine procedure. To improve numerical stability and local minima, a pyramidal approach starting with the coarse scale is introduced. After the displacement field has been estimated at a coarse level (tolerant to displacements of large magnitude but inaccurate), the finer details of the next pyramidal level are considered and the displacement field is corrected (with better precision, and lower tolerance to large magnitude). The P-DVC problem is convex only when the displacement is small as compared to the characteristic length of the projected texture (e.g., the auto-correlation length). This multi-scale analysis provides an "adaptive convexity" in the neighborhood of the solution. To generate the coarse description, projections are convoluted with a Gaussian kernel whose characteristic length can be adapted to the displacement amplitude. Thereby, high frequency patterns of the projection and reconstruction are erased and the P-DVC functional is driven by the global shape motion. The projections can also be down-sampled, in space and time to accelerate the procedure as the resulting images carry about the same information as the filtered one. For each scale, the procedure is computed until stagnation of the residual fields and displacement variation. The space kinematics is modeled as that of a finite element mesh on the pulmonary pleurae, meshed from the reconstructed volume intensity (obtained by thresholding and a closing morphological operation). Because the initial volume is corrupted by motion blur, the position of the envelop nodes is updated while the volume edges are more accurately reconstructed with the procedure iterations.

Experimental studies

Breathing simulations from multiple static chest reconstructions

The above proposed method is applied to a test case that has been constructed to be realistic and to provide a known ground truth. It is based on 10 real reconstructions of a patient's chest with phase-triggered 4D CT-scans [START_REF] Vandemeulebroucke | Spatiotemporal motion estimation for respiratory-correlated imaging of the lungs[END_REF][START_REF] Vandemeulebroucke | The POPI-model, a point-validated pixelbased breathing thorax model[END_REF] of the respiratory motion acquired at the CREATIS laboratory, Léon Bérard cancer center and available online ‡. ‡ https://github.com/open-vv/popi-model, with data from [START_REF] Vandemeulebroucke | The POPI-model, a point-validated pixelbased breathing thorax model[END_REF] These volumes were obtained from a time sampling based on the acquisition of a respiratory surrogate signal (Pneumo Chest pressure belt). Each one consists of 256 × 256 × 141 voxels. From this set of tomographic images, a synthetic moving CT-scan is assembled as follows:

• A respiratory phase (rounded to the closest of the 10 available phases) as a function of acquisition time (or equivalently projection angle) is first constructed as shown in Fig. 2(a). The same breathing state is used for 8 successive projections in order to have a realistic number of breaths in the entire experiment. This choice has almost no effect on the following developments.

• The synthetic sinogram is constructed by projecting at each angle (chosen to be uniformly sampled at 1 degree interval, over one full rotation, resulting in 360 radiographs) the volume corresponding to the respiratory phase.

While those volumes where acquired using an helical CT scan, the simulated moving patient is generated with a beam covering most of the patient's chest (i.e., CBCT conditions). For simplicity, faster computations and to avoid removing the top and bottom part of the volume due to the cone beam geometry, a parallel beam geometry is chosen here. The difference between parallel and cone projections only lies in the projection operator Π θ but has no further consequence on the current procedure.

To illustrate the motion, Figure 2(b) shows a composite image, where one such projection (90°) at respiratory phase # 1 is built with green color, and, for the same angle, the projection of the respiratory phase # 7 is attributed the complementary magenta color. Hence, all regions not affected by the motion appear as gray shades, but moving parts appear either in magenta or green. This shows that the displacement is essentially vertical and localized at the diaphragm. The obtained sinogram is used as an input of our procedure.

Because the bottom part of the scan moves with the breathing, the mean intensity of the image, κ(t) computed over a region that contains mostly the bottom part of the image, varies and is assumed, at first order, to be related to the phase of motion. This measured signal κ(t) is hence used to regularize the time evolution. In addition, the absolute |κ(t)| and squared κ(t) 2 values of this signal and two periodic functions of the rotation angle, sine and cosine, (that could be related in real cases to scan rotation issues) are added to the time function library. For the present case study, the space kinematic model is embedded into a multi-scale framework: at first, a simple C8 element (i.e., 8 noded-cube element with tri-linear interpolations) composed of 24 degrees of freedom (i.e., a total of 96 dofs for a ij for the 4 time steps) is considered. Then, when a coarse reconstruction has been obtained, a finer mesh composed of 579 nodes with T4 elements (i.e., 4 noded-tetrahedra elements with tri-linear interpolation) is used resulting in 6,948 degrees of freedom. The root mean square (RMS) of the projected residuals reaches a plateau after 4 iterations (an iteration being both the reconstruction with compensated motion and the projection-based measurement). The volume was reconstructed with N DynART = 10 at each scale.

Clinical CBCT with motion

The method is applied to a patient cone-beam CT (CBCT) scan acquired for 3D image-guided radiotherapy at the Léon Bérard cancer center. The CBCT scanner is integrated with the linear accelerator (Elekta Synergy; Elekta Oncology Systems Ltd., Crawley, West Sussex, UK). The flat panel recorded radiographs with 512 × 512 pixels of 0.8 × 0.8 mm 2 at a frame rate of 5.5 fps over a full revolution. The distance "source to isocenter" is 1000 mm, the distance "source to detector" is 1536 mm. Image acquisition was performed at 120 kV, 25 mA and 1 mAs/projection. The sinogram is composed of 676 projections equally spaced over a 360°revolution. The flat panel detector of size 41 cm×41 cm is offset with respect to the line defined by the source and the center of rotation and images only a part of the entire chest (a projection is displayed in figure 9 showing the missing area). A partial reconstruction procedure is hence required. The different geometrical parameters required for the cone beam reconstruction are given by the tomograph (position of the source and detector for each radiograph).

During the acquisition, the phase of approximately 35 irregular breaths of the patient are extracted from the projections. Because of the corresponding large magnitude (10-20 mm at the diaphragm) of breathing motion, standard reconstruction procedures lead to a blurry reconstruction around the diaphragm area.

The sought space-time kinematics is expressed in a reduced basis composed of selected space and time functions. A scheme of the regularization is shown in figure 3. The unstructured mesh is composed of 192 nodes with 532 T4 elements. In addition, the lungs may slip on the rib cage [START_REF] Vandemeulebroucke | Automated segmentation of a motion mask to preserve sliding motion in deformable registration of thoracic CT[END_REF][START_REF] Delmon | Recalage déformable de projections de scanner X à faisceau conique[END_REF][START_REF] Zhang | A biomechanical modeling guided CBCT estimation technique[END_REF] and hence a discontinuity of the tangential velocity may appear. Reciprocally, non-zero normal displacements of the mesh envelop are removed by projecting the full displacement on axes normal to the mesh. The time kinematics could be modeled with 5 simple time functions, including (1-2) a sine and cosine of the rotation angle, (3-4) constant and linear evolution and (5) a breathing signal obtained from an independent intensity tracking measurement on projections. This latter time evolution is expected to be the one containing most of the breathing motion. As such, it was chosen to update this breathing signal (5) with the relaxation procedure. With a normalization for the amplitude (by convention), 675 degrees of freedom (dofs) are hence added in the identification procedure. The 4D motion is finally composed of 192 nodes thus 576 spatial dofs, 4 time functions and 675 temporal dofs from the time relaxation generating a total of 2979 space-time dofs. The space mesh and time evolutions are shown in figure 3.

Results

The measurement part was performed on CPU: 20 cores @ 2.20 GHz while all projection operations are done on GPU: Nvidia Tesla K80 (with ASTRA toolbox). What is computationally demanding is volume warping at each projection/back-projection when performing the motion-compensated SART procedure and the generation of the projected sensitivity fields (e.g., image variation with respect to a parameter change).

Results on the numerical simulation

Slices of the reconstructed volume before and after are shown in figure 4. After the displacement field correction, the diaphragm is visible while the reconstruction without motion correction was more fuzzy. The 10%-90% rise response of the diaphragm edge decreased from 11.2 voxels (22.4 mm) to 2.9 voxels (5.8 mm). (i.e., the part of the initial data that is not explained by the reconstruction approach), and this is precisely captured in the projection residual (see figure 6 before and after correction). The projection residual field consists of noise and artifacts affecting the acquisition but includes also all model errors (approximate kinematic model, inaccurate interpolations, etc). It can be seen that the residual values are very high in the beginning. On the bottom part, large positive and negative values are visible resulting from the oscillating respiratory motion. At the end of the procedure, most of the oscillating motion has been erased meaning that the estimated displacement field was trustworthy. The RMS was 2.9 times smaller after correction.

Displacement field

The reconstruction quality can be evaluated by comparison with respect to the reference 3D volume at the same respiratory phase. The estimated displacement field can be compared with "standard" DVC procedures (because we have access to the reference 3D volumes). Let us underline that both of these quality evaluations are generally not accessible, but they motivated the construction of the present test case. The estimated 4D displacement field can be compared with full 3D DVC procedures using the 3D reference volumes. Using the Correli-3.0 software [START_REF] Leclerc | [END_REF] that follows a similar intensity based measurement in 2D, the displacement field was obtained with the same 579-node mesh but independent time states. The comparison of the vertical displacement field using the dynamic procedure and DVC is presented in figure 7 and shows very similar results. The mean and standard deviation of the difference of the vertical fields is respectively -0.21 and 0.92 voxels. This comparison has however to be considered as informative and not as a validation because (i) the cost functions of the two methods are different and (ii) some nodes have a low sensitivity because of the poor contrast and hence they are susceptible to noise, although this poor determination of the displacement has no perceptible consequence on the reconstruction quality precisely because of the low contrast of these zones. Reconstructed volume The quality of the reconstructed volumes, for the 10 respiratory phases, can be evaluated from their comparison with the initial reference volumes used to generate this test-case. The norm of the difference is mostly controlled by a region containing the diaphragm (see figure 8). In order to quantify the error, these norms are normalized by the noise variance on the same domain. This noise is estimated from the residual obtained after a fine 3D registration between two arbitrary 3D references (two different respiratory phases). This value corresponds to the best residual that can be obtained considering an extremely fine full 3D kinematic correction, and hence the normalized residual norm for a perfect correction should approach unity. It can be seen that the corrected dynamic procedure produces a residual error much lower than the standard static reconstruction. The oscillating signature of the red curve due to the periodic motion has been, for the most part, cancelled. The error level flattens to about 1.5 and is no more than 50% higher than a perfect correction (with a much finer mesh). This result is to be compared to the error obtained with no motion correction which can reach values higher than 3.

Application to clinical CBCT data

Figure 9 shows the re-projection of the volume at angle θ(t k ) = 15°and the projected mesh envelop that defines the region of interest at this angle. Although the missing part of the detector now appears (as the full volume was reconstructed from all angles), the residual field of this area cannot be computed as the projections are not defined. However, thanks to the time regularization, the sensitivity of the motion on those masked areas is obtained from the other angles where the information becomes visible. The projected residual field and displacement field increment become stationary in 10 iterations at scale 3 and 4 iterations at scale 1. N DynART = 10 and N DynART = 30 were chosen as stopping criterion for the reconstruction respectively at scale 3 and 1. At each iteration, the nodes of the mesh are re-positioned onto the lung external boundary. The procedure takes approximately 6 hours to reach stationarity at scale 3 (projections reduced by a factor 1/3 in each dimension) and 24 hours at scale 1.

The estimated 4D displacement field is shown figure 10. As expected, the vertical displacement is very large along the diaphragm (18 voxels at scale 1 corresponding to approximately 15 mm) and follows essentially the breathing motion. The normal displacement has been erased allowing sliding of the pleural cavity [START_REF] Zhong | 4d cone-beam ct reconstruction using multi-organ meshes for sliding motion modeling[END_REF]. The refined breathing signal is in phase with the original one, estimated from projections and used for initialization, however, amplitudes are observed to be markedly different. The reason may be that the initialization was performed with a tracking of the diaphragm in the projections, thus submitted to perspective effect due to the cone beam while rotating.

Reconstruction improvement

The correction naturally focuses on the large displacement areas, i.e., the bottom of the chest around the diaphragm. It can be noted that the diaphragm is reconstructed much more accurately after the displacement correction. Initially smooth and blurry, the edges get much sharper, as expected. Some details of the bronchus become visible. The comparison of the two reconstructions is shown figure 11 and highlights the quality improvement. The plot in figure 11 shows the intensity profiles at the bottom of the lungs (visible in the slices above), before (dashed curve) and after (plain curve) which highlights the correction of the motion blur. Initially spread over 12 voxels (from z = 152 to z = 164 voxels), also measured by the ± 6 voxel displacements of the diaphragm figure 11, the diaphragm transition is reconstructed sharper. The 10%-90% rise response decreased from 13.2 voxels (10.6 mm) to 1.4 voxels (1.2 mm). Segmenting the lungs for example would be much easier and reliable after correction.

The metric based on the gradient norms gave 2.8•10 -7 and 3.4•10 -7 for respectively the standard and corrected image. As an order of magnitude, when filtered by an averaged squared kernel of size (3x3), the the L2 norm of the gradients from the corrected filtered image is 2.7 • 10 -7 which is the value for the standard unfiltered image. This improvement highlights the sharpness evolution of the reconstruction.

The true metric to validate our procedure is the projected residual field that corresponds to the reconstruction model error. This residual sinogram can be shown in 3D in figure 12 before and after the dynamic tomography procedure. The reconstruction inconsistency due to motion is clearly visible in the left image with the positive and negative patterns. Those error parts are mostly corrected when taking the kinematics into account and the RMS decreases by a factor of 2.31 in the region of interest described by the mesh volume (expressed in % of the gray level amplitude of the projections).

What remains in the left corner part of the final residual field is an error located at the intersection between the diaphragm and the bottom rib (the acute lung corner in the bottom left of the coronal slice in figure 11). The kinematics at this position is complex (slip and recovery) and was not perfectly captured. Moreover, this part of the diaphragm is not always visible depending on the breathing and rotation angle. This residual field study gives indications on where and how to improve the kinematic model. The model could be further enriched as suggested by the interpretation of the final residual, adding space or time dofs.

Discussion and Conclusions

The present work follows the trend that has been outlined in the literature review of the introduction on motion corrected tomography. An iterative algorithm is proposed to determine both morphology and kinematics with no additional acquisition cost as compared to a classical scan. The procedure consists of two staggered parts within a pyramidal coarse-to-fine multi-scale framework [START_REF] Edward H Adelson | Pyramid methods in image processing[END_REF]: (i) a motion compensated procedure extended from a SART algorithm and (ii) a displacement field measurement from the reconstructed volume and the initial sinogram.

The procedure is tested on (1) a numerical simulation of breathing from static chest CT and (2) a CBCT of a breathing patient during radiotherapy at different phases of the respiratory motion whose reconstruction leads to large motion blur artifacts. The multi-scale procedure, performed for the volume reconstruction, allows large motions to be captured. Finally, a large part of the motion blur artifacts was removed from the reconstruction, leading to small root mean squared errors computed on the projection residual fields. The reconstructed volume and measured displacement field of the numerical application were compared to ground truth, thereby validating the procedure.

The originality of the present work lies in a very accurate description of the kinematics, using a finite element mesh that offers many degrees of freedom. This progress leads to two main difficulties:

• First, the mesh is to be suited to the anatomy, however, initially the chest morphology is unknown, so that not only the 3D image and its motion are to be determined simultaneously, but also the high discretization of the motion and suited boundary conditions. In our specific application, the thoracic cage is mostly static, and hence the velocity field at the boundary of the meshed lungs is constrained to be parallel to lung boundary, allowing for a possible slip. This calls for several iteration steps along which the morphology is progressively described with a better fidelity, the motion discretization is rendered more and more accurately, ending with a fine mesh. The simultaneous determination of these three items -morphology, mesh and motion -is an original problem.

• Second, with the above spatial description of the displacement field, involving many degrees of freedom, each of which being susceptible to vary in time all along the scan, makes the problem severely ill-posed. This difficulty is circumvented by a model reduction technique whereby "modes" are considered. Each mode is a separated form being the product of one spatial motion (supported by the mesh description) and a temporal part that describes the amplitude evolution of this mode. Several types of modes are introduced, some where the temporal shape function is prescribed (such as sine or cosine of the rotation angle, or a low order polynomial in time), and some more complicated ones, where the time evolution is that of breathing. The latter can be inferred initially by the apparent motion of the diaphragm, and -through iterations -space and time parts are progressively finetuned to match the observation. By the suited selection of these modes, the number of degrees of freedom is severely reduced, and makes the problem tractable. In some way, one may thus add a fourth category of unknowns, in addition to morphology, mesh and motion, that are the considered modes.

It is interesting to control the procedure by studying the projected residual fields that contains all features not explained by the projection model (e.g., detector noise, non-linear absorption, unsuited kinematic model e.g., due to the chosen discretization, etc.). This residual field allows for (un)validating the kinematic model. This feedback is essential when encountering an unexpected kinematics. Regularization should be considered as evolutionary and any feature that can be understood from the residuals should help to enrich the model up to the stage where the remaining residual is deemed acceptable. Additional degrees of freedom allowing to model and correct the heartbeat motion is an interesting perspective.

The proposed method is totally non-intrusive as it does not require more than a standard acquisition (i.e., a single revolution, and a standard acquisition protocol). Therefore, this method can be applied to any type of scanner. Of course, all additional information (from sensors, from other modalities, from experts, etc.) would help the analysis and accelerate the optimization of the process, making the problem better posed. From an hardware side, using faster scans would allow reducing the motion blur that may happen in the projections and would improve the identify temporal motion model.

In the proposed approach, it is not required to initiate the computation with a subset of radiographs that could be considered as (approximately) still, (because they would originate from a similar respiratory cycle phase). This is not a small difference with respect to existing approaches such as SMEIR, because for instance, the mesh that is used for describing the motion, has to align perfectly the organs in order to allow for free slip at boundaries in [START_REF] Zhong | 4d cone-beam ct reconstruction using multi-organ meshes for sliding motion modeling[END_REF][START_REF] Zhang | A biomechanical modeling guided CBCT estimation technique[END_REF]. But when no decent 3D volume has been reconstructed yet, it is not possible to proceed with a fine adjustment of the mesh. The volume, mesh, and displacements are all three very intricated and need to be progressively improved better and better defined. As a result, the proposed method does not require a (quasi-)periodic motion nor extraneous information such as the measurement of the phase. Without this initialized motion, its measurement from truncated partial projections became impossible without the introduction of time regularization.

One current limitation for clinical practice is the substantial important computation time that includes the computation of all the sensitivity fields. In addition to hardware improvement (faster computational resources, GPU, parallel computation, etc.), the lightening of the computation operations (development of reduced breathing models, minimization of the projected gradients updates, etc.) will highly decrease the computation time. During the iterative process, improving the reconstruction quality and segmenting the moving organs with image processing (spatial smoothness, total variation) or Deep-Learning techniques would also allow to better mesh and model the breathing motion leading to faster convergences.

This procedure opens up new avenues in dynamic tomography enabling to image moving patients, even if the displacement is large and with no additional data acquisition compared to a standard CT-scan protocol. The clinical development proposed in [START_REF] Rit | On-the-fly motioncompensated cone-beam ct using an a priori model of the respiratory motion[END_REF] with an a-priori motion model could be extended without the need of the prior 4D-CT scan. Moreover, any additional sensors and detectors (e.g., force or displacement measurements) could be used transparently in the above method to constrain further the analysis. Using dual source scanners is also a promising perspective, especially when the sources are orthogonal [START_REF] Wang | A scheme for multisource interior tomography[END_REF][START_REF] Flohr | First performance evaluation of a dual-source CT (DSCT) system[END_REF] so that displacement fields could be completely estimated for each acquisition time (i.e., providing a better conditioning of the problem). Finally, any additional regularization, either for the volume reconstruction (e.g., assumptions on discrete phases [START_REF] Batenburg | DART: a practical reconstruction algorithm for discrete tomography[END_REF][START_REF] Roux | Efficient binary tomographic reconstruction[END_REF], total variation (TV) regularization [START_REF] Emmanuel | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF], dictionary of commonly encountered tissue morphology, etc.) or for the kinematic model can be easily included within the proposed algorithm. Such developments are very relevant for practical applications. Nevertheless, they were not included in the presented test case, in order to focus on the key ingredients for the success of the reconstruction. However, additional acquired or prior information is expected to be beneficial for faster convergence, more accurate results, or for further reducing the radiation dose (e.g., acquiring fewer projections for an equivalent result).

II. Dynamic reconstruction algorithm

The global dynamic reconstruction algorithm consists of two main steps, detailed hereafter, which are run alternately, up to stationarity.

Algorithm 1 Proposed motion-compensated reconstruction

Initialization n ← 1 Initialization f (1) ← 0 Choose a permutation, π, over N t acquisition times for i ← 1 to N DynART do for t ← 1 to N t do τ ← π(t)

f (n) (x, τ ) = f (n) (x + u(x, τ )) warp ρ(r, τ ) ← p(r, θ(τ )) -Π θ(τ ) [ f (n) (x, τ )] ∆ f (n+1) (x) = B θ(τ ) [ρ(r, τ )] ∆f (n+1) (x) = ∆ f (n+1) (x -u(x, τ )) unwarp f (n+1) (x) = f (n) (x) + ∆f (n+1) (x

Figure 1 .

 1 Figure 1. Scheme of the proposed procedure that requires only a standard scan acquisition of the moving scene. It consists of two staggered blocks: 1) motion measurement based on the current approximation of the reconstructed volume and 2) motion-corrected reconstruction based on the current determination of the motion.Initially, the motion may either be set to 0 or benefit from prior knowledge.

Figure 2 .

 2 Figure 2. Synthetic test case generated from 10 real CT-scans. (a) Chosen sequence of the respiratory phase as a function of acquisition angle. The entire test is composed of 4.5 respiratory cycles. (b-c) Composite image formed by two projections for the same angle 90°, and two different respiratory phases, where each projection is given a complementary color (green and magenta). (The two states correspond to the red dots in figure (a).) (c) Composite image from a 3D slice generated from the same respiratory phases.

Figure 3 .

 3 Figure 3. Space-time regularization of the motion composed of (a) a 192-nodded mesh with tetrahedral elements and (b) 4 simple time evolutions and one fully identified temporal function initialized with the breathing signal.

Figure 4 .

 4 Figure 4. Sagittal and coronal slices of the reconstructed chest without (left column) and with (right column) motion correction (the vertical red lines indicates the slice positions). It can be observed (especially in the zoomed area) that the diaphragm is better reconstructed and sharper (as pointed out by the two white arrows). On the sagittal view, a horizontal plane that polluted the reconstruction disappeared (white arrow). The bottom plot shows 2 intensity curves (in arbitrary units) extracted from the left and right volume on the zoomed area (yellow lines). The dashed curve, corresponding to the non-corrected volume, is smooth. The solid curve is sharper thanks to the motion correction.

Figure 5 .

 5 Figure 5. (a) Scale 1 mesh with 579 nodes is plotted highlighting the lungs (in red); (b) Color-coded vertical displacement field is shown for the acquisition at angle 48°. The displacement is expressed in voxels (Note that voxels are slightly anisotropic, of size [1.95, 1.95, 2.0] mm along [x, y, z]).

Figure 6 .

 6 Figure 6. Projection residual fields (difference between the acquired sinogram and the set of projections computed from the reconstructed volume at the corresponding deformation stage) before (left) and after (right) correction by the estimated displacement field. The oscillating positive and negative patterns are the signature of the 4 respiratory motion periods. After correction by the estimated displacement field, a large part of the error has been corrected.

Figure 7 .

 7 Figure 7. Estimated vertical displacement field (in voxels) for angle 48°using (a) a standard DVC procedure and (b) the dynamic approach using the same color scale. The red parts correspond to the motion of the diaphragm .

Figure 8 .

 8 Figure 8. Norm of the differences between the reference volumes and (a) the corrected volume using the dynamic reconstruction procedure and (b) using a standard static ART reconstruction.

Figure 9 .

 9 Figure 9. (a) Acquired projection p(r, t k ), with θ(t k ) =15°and (b) simulated projection from the volume reconstruction. The rotation axis is highlighted in red. The blue line is the envelop of the projected mesh. The projection is incomplete and involves an insensitive area. The top and bottom part of (b) are subjected to the projected volume due to missing data in areas outside the field-of-view.

Figure 10 .

 10 Figure 10. (a) Vertical displacement field of the mesh envelop at the first maximal exhale position (projection 7), expressed in voxels (1 voxel ↔ 800 µm) and (b) measured (black) vs estimated (dashed blue) shape of breathing signal. The amplitude of this motion is on the corresponding spatial field.

Figure 11 .

 11 Figure 11. (a) Comparison of two slices of the reconstruction, using (left) a standard ART procedure and (b) the proposed dynamic reconstruction method with the identified displacement field. It is to be noted that the diaphragm area is much sharper in the right corrected slices and bronchus (highly subjected to motion artefacts) are more visible. The red lines indicate the slice locations and the white segment corresponds to the place where intensity profiles are extracted. On the bottom plot, the intensity profiles (in arbitrary units) on the diaphragm area, (dashed) before and (plain) after the dynamic reconstruction procedure. The diaphragm interface is much sharper for the corrected procedure.

Figure 12 .

 12 Figure 12. Sinograms of the projected residuals before (left) and after (right) the kinematics identification and correction. The plot is focused on the diaphragm part. The color amplitudes and transparencies are identical for the two plots. The high positive and negative patterns (signature of the 35 breaths) visible in the left image have been mostly erased.

  ) Implement additional constraints on f (n+1) (e.g., positivity)n ← n + 1 end for end for Algorithm 2 Projection-based kinematic measurement u(x, t) ← u 0 (x, t) Displacement initialization for Measurement iteration ← 1 to N iteration do p ← pλ Gaussian filtering f λ (x) ← Motion-compensated reconstruction(p, u) Algo. 1 Compute ρ l (r, t) Projection residual for i ← 1 to N s do for j ← 1 to N t do f (x, t) ← f λ (x + u(x, t)) Volume advection χ i (x, t) ← Π θ(t) [φ i (x)∇ f (x, t)] S ij (r, t) ← χ i (x, t)σ j (t)

				Projected sensitivities
	end for
	end for	-1
	ã ←	S ⊗ S	ρS	Motion identification
	r,t	r,t
	a ← ã			Erasing normal displacements
	for i ← 1 to N Sensitivity of the first mode
	end for	-1
	b ←	S ⊗ S	ρ S	Time identification
	r	r
	end for

s do S(x) ← Π θ(t) [a i1 φ i (x)∇ f (x, t)]

Acknowledgements

The authors would like to thank Dr. Matthieu Vitse and the SciViz group, for having performed most of the 3D rendering in this article.

Ethical statement

This research study was conducted retrospectively using non identifiable human subject data. The POPI study has been approved by the local research ethics committee and data have been anonymously open to the community since 2007. Applicable law and standards of ethic have been respected in accordance with the principles embodied in the Declaration of Helsinki and in accordance with local statutory requirements.

Supplementary material I. Projection-based measurements

The P-DVC measurement consists of solving, with a Newton-Raphson algorithm, the following equation

with µ ∈ U = P × S, P and S being respectively the space generated by the φ(x) space functions and time evolution σ(t)

The projected error field,

with p(r, t) the acquired projections, is linearized, with f (x, t) = f (x + u(x, t))

Then, with u(x, t) = ij a ij φ i (x)σ j (t), the minimization with respect to a leads to the solution of the linear system

with S the projected sensitivity fields (signature in the detector space of the dof sensitivity fields)


Marc Neveu

Jose C Aponte

Julie Castillo-Rogez

Bethany L Ehlmann

Heather B Franz

Christopher H House

Erwan Mazarico

Alex Pavlov

Vassilissa Vinogradoff

Kris Zacny

José C Aponte

Exploring Solar System Organic Chemistry Evolution through the Surfaces of Ceres and Large Asteroids

établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

• Accessing and characterizing solar system organic material is crucial for understanding the processes by which our Solar System became habitable.

• Recently exposed surfaces in craters of large asteroids (e.g., 10 Hygiea, 24 Themis, Ceres), with a large diversity of ages and abundant organic compounds, provide a natural laboratory to study Solar System organic chemistry and the potential for habitability throughout time. For the dwarf planet Ceres, surface exploration is possible at Discovery and New Frontiers costs. Ceres has likely hosted habitable environments, possibly prebiotic ingredients, and perhaps even incipient life at some time in its history.

• Ceres' organic material is abundant enough to be seen from orbit. Ceres' surface spans (1) geologically old areas comprising ~10 wt. %t carbon, (2) a few locations with higher concentrations of more diverse organic compounds, and (3) bright areas (mostly sodium carbonate), some at most a few 10 5 -10 6 years old and perhaps much younger. This young surface age is an opportunity to examine organic matter evolution, including both products recently exposed from Ceres' interior and those modified by radiolysis damage.

• Radiative processing has been a major factor shaping organic evolution throughout Solar System history, both degrading compounds (e.g., molecular biomarkers) and increasing organic complexity. Laboratory investigations of both aspects, including on products of aqueous reactions, are needed to advance our understanding of habitable worlds.

Follow the organics

Organic material enables life. Accessing and characterizing solar system organic material is key to understanding the processes by which Solar System worlds became habitable and, for some, inhabited. Understanding organic evolution from Solar System formation to the present requires identifying precursor molecules and elucidating pathways of synthesis and modification.

Chief among these ar e photo-and radiolytic processes and aqueous alteration or metamorphism (Chang & Bunch 1986;Shock & Schulte 1990;Aponte et al. 2017), which produce organic compounds relevant to prebiotic chemistry . Solar radiation was the largest energy source in the young Solar System, and affected organic matter before accretion in planetesimals and then on the surfaces of planetary bodies. Organic material synthesized and modified by radiolysis, photolysis, or photoexcitation, further reacted in planetary bodies through the interaction of water and rock fueled by radiogenic and impact heat.

Ceres and like asteroids, record keepers of organic evolution

Understanding processes that drive the synthesis and evolution of organic material that contributes to making worlds habitable requires leaving Earth. On our planet, the processes that led to the organic makeup that preceded the emergence of life but paved the way for it have been hopelessly difficult to disentangle from the imprint of biology. Arguably, the worlds best suited for investigating abiotic organic evolution are those on which such organic synthesis is ongoing (to minimize biases introduced by subsequent modification or destruction processes) and where organic compounds are abundant (to ensure detection and characterization of a broad range of compounds, which helps elucidate synthesis pathways and the full suite of products). Although organic synthesis can proceed in atmospheres, including likely that of early Earth, the level of chemical complexity that led to life likely required additional endogenic processing of compounds predating accretion or photochemical products in liquid water.

Very few accessible (e.g., surface or plume) solar system locales meet all three criteria (Fig. 1), in part because they act against one another. For example, Earth's ongoing endogenic abiotic organic synthesis feeds life; its consumption results in a small proportion of abiotic organic compounds. A record of solar system organic matter is preserved, albeit likely degraded, in carbonaceous chondrites, but less-modified products of recent syntheses are in all likelihood absent. All three conditions are likely met in the interiors of ocean worlds but masked by atmosphere-surface physical and chemical interactions or a thick, diluting hydrosphere. Large protoplanets, the survivors of those that contributed to delivering water and organic material to the terrestrial planets, have organic material that is both endogenous and abundant. Among carbonaceous protoplanets, Ceres is the best-known representative and perhaps the only one both with a record of ancient organic material and where endogenic synthesis might continue today , as detailed below.

Figure 1. Large asteroids such as Ceres are the only solar system bodies preserving a record of organic synthesis spanning both early solar system material and the present day. They may be our only window into solar system organic chemistry evolution.

Surface regions on Ceres span ages from the early solar system to the present day:

• Most of the dark surface (Fig. 2a) was likely emplaced early in solar system history and contains up to 20% carbon of organic origin. Its mineralogy "appears to be dominated by products of rock-fluid interactions, such as phyllosilicates -some of which are NH 4 -bearingand carbonates.

[...] Ceres's mineral and elemental data can be explained by the presence of carbonaceous chondritic-like materials (~50-60 vol%), possibly due to infalling asteroids, admixed with aqueously altered endogenic materials that contain higher-than-chondritic concentrations of carbon. [...] Ceres's surface may contain up to 20 wt% of carbon, which is more than five times higher than in carbonaceous chondrites. The coexistence of phyllosilicates, magnetite, carbonates and a high carbon content implies rock-water alteration played an important role in promoting

Figure 2. Ceres records billion-year-old organic synthesis in the presence of liquid water on most of its surface (a). Products of ongoing synthesis may be available at select areas such as Cerealia facula in Occator (b). Even measurements made from orbit distinguish different types of organic compounds at the surface, such as aliphatic species at Ernutet crater (c).

• Other areas of undetermined age are organic hot spots. Ernutet crater, despite not showing any outstanding morphology, is the site of a patch rich in aliphatic (chained) and possibly aromatic (ringed) organic compounds that straddles part of its rim (De Sanctis et al. 2017, 2019). This is surprising because organic compounds tend to degrade into 6-carbon ring structures (ultimately, graphite), and/or become disordered (amorphous) carbon due to space weathering (see next section). The presence of not fully aromatic or disordered organic carbon seen from orbit thus suggests fresh sourcing from recent excavation from the subsurface where these Other protoplanets, too (Fig. 3) . Although Ceres is the only large carbonaceous airless world studied in detail, it may not be the only host of abundant endogenic organic material with ages spanning solar system history. The third and fourth largest main-belt objects, Pallas and Hygiea , also have surface compositions (inferred from remote measurements) reminiscent of carbonaceous chondrites. These 500 km-class objects are highly unlikely to preserve any liquid by radiogenic heating alone (since even Ceres is at best on the brink of being fully frozen; Similarly, sampling carbonaceous protoplanet compositions that differ from Ceres' may require visiting other objects. This brings about a trade between the diversity of organic syntheses that can be accessed and the scope of the mission(s) required to access it. In rough order of increasing science return and increasing cost, possible options include:

• Multiple remote-sensing missions to protoplanets other than Ceres (Discovery-class), provided evidence for geologically seemingly ancient and recent organic material can be uncovered with upcoming telescopes (e.g., extremely large ground-based telescopes and JWST ). The quality of the data obtained on such objects even with existing assets (Rivkin et al. 2019;Vernazza et al. 2020; Fig. 3 ) suggests the potential of next-generation telescopes to obtain such evidence.

• A single-site in situ investigation at Ceres (Discovery-class) with the capability to determine the organic matter composition of one or more compositional unit(s) of Ceres' surface.

• A mobile in situ investigation at Ceres (New Frontiers-class, see PMCS report) with the ability to carry out detailed organic matter compositional analysis.

• Sample return from Ceres (New Frontiers-class, see PMCS report). The PMCS report features sampling from a single site in fresh evaporites but additional sites may be feasible.

Where the "knee in the curve" lies in science return per dollar remains to be demonstrated by mission concept studies. The Ceres PMCS report provides essential input in this regard, as well as in the types of measurements required and the spaceflight readiness of measurement methods.

Radiation processing of surface organic material

Today, water-rock interaction has all but ceased but its organic products exposed at asteroidal surfaces continue to undergo radiative processing. Disentangling radiation processing post-surface emplacement from prior drivers of organic evolution presents a challenge to understanding the emergence and persistence of habitable environments in the Solar System.

Radiation giveth organic material. Organic compounds can be formed from carbon-bearing ices (e.g., CO, CO 2 , CH 4 , HCN) through chemical reactions driven by radiolysis or UV photolysis. This has been studied most extensively for extremely cold environments, such as comet-forming regions of the Solar System, molecular clouds, and the interstellar medium (e.g., Takano et al. 2007;Gontareva et al. 2009;Hudson and Loeffler 2013;Bergantini et al. 2017;Bibang et al. 2019). Radiolysis by galactic cosmic rays may have formed organic material at the dark poles of the Moon and Mercury (Crites et al. 2013). In addition, complex biochemical compounds have been synthesized experimentally from UV irradiation of precursor molecules in the presence of different mineral samples relevant to the surfaces of small bodies in the Solar System (Simakov & Kuzicheva 2005).

Radiation taketh organic material away. Ion irradiation of carbonaceous material analogous to interstellar materials clearly causes destruction of aliphatic C-H features over time (Godard et al. 2011). An application to Ceres is shown in Fig. 4 , in which experimental results on amino acid degradation by γ-ray irradiation are combined with estimates of radiation fluxes at Ceres to calculate the fraction of molecule survival over time. Preservation is increased at a few cm depth because by 5 cm, most solar irradiation is absorbed, while galactic cosmic radiation (GCR) decreases slowly between 5 and 10 cm. These results suggest that abundances of key organic compounds such as amino acids, if emplaced within a few cm of Ceres' surface in ice-poor regolith, should remain within the same order of magnitude for tens of millions of years; i.e., not significantly changed (relative to the limits of detection of measurement techniques) in recent surface locations such as Occator crater where time scales preclude significant impact gardening. They also show that the mineralogy in which organic compounds can be embedded can significantly affect their survival; experiments to date have only scratched the surface in exploring this essential parameter space. Pavlov et al. (2012Pavlov et al. ( , 2014)); the GCR flux was increased by 10% compared to the flux at 1 AU due to less modulation by the heliosphere. No water was assumed; for organic molecules embedded in ice, degradation would be faster. The experimental radiolysis data used were obtained from γ-ray irradiation; protons tend to increase degradation for the same dosage. These data are for local noon temperatures; degradation is expected to be substantially slower at Occator's ≈ 160 K average diurnal temperatures (Hayne & Aharonson 2015).

Reconstructing the history of organic material: experiments needed

Compositional and contextual measurements of surface organic material are underpinned by observational techniques of astronomy and geology, but reconstructing the formation and evolution of organic material also requires experimental chemistry. Even when in situ or sample return investigations establish the concentration and molecular distribution of organic matter, experimental irradiation and hydrothermal processing of material analogous to that inferred to be present on Ceres and like asteroids is needed to firmly assess the origin of such compounds. The need for more analog experiments . A valuable attempt to reconstruct the early and current organic content observed on worlds such as Ceres requires carrying out systematic studies of the formation of organic compounds through a combination of irradiation and hydrothermal processing experiments. Since organic matter can be intimately mixed with minerals (typically phyllosilicates but also salts), as observed both in carbonaceous chondrites and on Ceres, we recommend scrutinizing the effect of minerals (e.g., ammonium phyllosilicates and salts in the case of Ceres) on the survivability of organic matter on irradiated planetary surfaces.

An example avenue would involve these steps: (1) simulate icy grains under conditions relevant to dense molecular clouds and protoplanetary disks, followed by (2) continued synthesis and modification of organic mixtures via hydrothermal processes post-accretion, and (3) re-irradiate to mimic the expected time a given surface has been exposed. An alternative avenue could involve continued irradiation at surface conditions in the presence of realistic analogous mineralogy. Generalizing experimental results to predict organic modification for different ice and mineral compositions, radiation doses, and hydrothermal conditions (e.g., temperatures) requires experimental modeling ( Fig. 4 ).

Isotopic clues.

The stable isotopic composition of organic material represents another key component in understanding its origin and history. Analyses of organic compounds in chondrites have revealed a wide range of values in isotope ratios, particularly those of H, C, N, and O (e.g., Kerridge 1985;Elsila et al. 2012), which often differ from those typical of terrestrial biomarkers. This observation is helpful in identifying products of abiotic organic chemistry on the meteorite parent bodies by screening out terrestrial contamination. Stable isotopes also offer a means to test hypotheses concerning the origins and modification histories of organic compounds and water on Ceres, e.g., whether the volatile inventory represents a single source incorporated during primary accretion or whether some compounds are derived from materials delivered to the surface later, such as by comet impacts. Stable isotope ratios should be key measurements for in situ missions to and sample return from Ceres (see white paper by Miller et al.). Incorporating isotopic measurements in the above supporting experimental avenues is thus essential.

Understanding the processes by which Solar System worlds became habitable requires exploration of Ceres (and like asteroids) and organic chemistry experiments essential to the interpretation of mission measurements.

  widespread carbon chemistry"(Marchi et al. 2019). The overall homogeneity of Ceres' surface suggests a global emplacement at an early time when liquid water was ubiquitous in Ceres' interior, which likely required heating by short-lived radionuclides in the first few million years (e.g.,Neveu and Desch 2015; Ammannito et al. 2016; Travis et al. 2018; Castillo-Rogez et al. 2019).• Select areas of the surface appear geologically very recent and likely bear organic material. Examples are the dome of Cerealia facula in the center of Occator crater ( Fig.2b) as well as other bright, recently exposed areas in Haulani and Oxo craters, and others. Some of these areas are at most a few 10 5 -10 6 years old (Nathues et al. 2017; Krohn et al. 2018; Scully et al. 2018) and perhaps much younger. These areas are rich in sodium-and ammonium-carbonate salts and chlorides (De Sanctis et al. 2020), and some bear water ice which is unstable below ~50º latitude on Ceres (Combe et al. 2016). Predicted sub-percent abundances of organic compounds (Neveu et al. 2017; Castillo-Rogez et al. 2018) are too low for remote detection but copious for in situ or sample return investigations, as demonstrated at Enceladus (e.g., Postberg et al. 2018). The present-day subsurface aqueous environments in which endogenic synthesis may be taking place (Hesse & Castillo-Rogez 2019; Ruesch et al. 2019) may be kept liquid by radiogenic or impact heating.

  compounds were synthesized (Bowling et al. 2020). Depending on the structure and composition of these aliphatic compounds, maximum abundances at Ernutet range between 5 area% and 65 area% (Kaplan et al. 2018); recent results suggest 15 wt.% (Vinogradoff et al. 2019).

  Neveu and Desch 2015; Travis et al. 2018; Castillo-Rogez et al. 2019), so any geologically recent or current organic synthesis on these bodies hinges on their recent bombardment by impactors able to generate local melting. However, the absence of bright spots in best-resolution images (at which Occator faculae are obvious; Vernazza et al. 2020) of Pallas and Hygiea does not bode as well for the presence of recently synthesized organic compounds. On the other hand, the existence of collisional families for both Pallas (Lemaître & Morbidelli 1994) and Hygiea(Zappalà et al. 1995) offers the prospect of sampling (old) material from their interior. These objects and others (e.g., the 200-km Themis ) show compositional variations relative to Ceres(Larson et al. 1983; Campins et al. 2010; Rivkin & Emery 2010; Vernazza et al. 2017; Rivkin et al. 2019; see white paper by Landis et al.). Such diversity in carbonaceous planetesimal compositions would help inform the diversity of conditions in early organic synthesis beyond the information provided by carbonaceous meteorites (e.g., within geologic context).

Figure 3 .

 3 Figure 3. Spatially resolved images and spectra of the largest protoplanets acquired in the past decade have enabled the identification of prime locations for investigating the origin and fate of organic compounds that make worlds habitable. Ceres image obtained by the Dawn mission; Pallas and Hygiea images obtained by Vernazza et al. Sampling a single versus multiple sites. At the sensitivity and spatial resolution of Dawn measurements, areas of Ceres' surface with newly synthesized compounds (e.g., Occator's faculae) and enriched in organics (e.g., Ernutet) are far apart, though the ubiquitous regolith at all locations contains ancient organic material at abundances greater than carbonaceous chondrites. Mobility on Ceres is expensive in terms of propellant but hopping between sites separated by tens of km is feasible within the New Frontiers cost cap (Castillo-Rogez et al. PMCS report).

  Example notable experimental results . While not fully representative of asteroid or Ceres surfaces, irradiation of carbonaceous chondrites produces spectral variations in the visible and near-infrared with inverse reddening/darkening trends for objects rich or poor in carbon(Lantz et al. 2015, Brunetto et al., 2018). This provides a means to constrain the exposure age of planetary surfaces by spectroscopy. In experiments simulating the inner structure of icy dust grains exposed to high-energy radiation, polypeptides were formed from single amino acids(Gontareva et al. 2009). Past aqueous alteration experiments have produced hundreds of molecules of a

	given compound starting from one precursor (Cody et al. 2011; Vinogradoff et al. 2018; Isono et
	al. 2019), leading to sequestration of organic molecules inside phyllosilicates (at least 7%)
	(Vinogradoff et al. 2020), to the production of insoluble polymeric carbon structures (Kebuwaka
	et al. 2013), and amino acids (Kebuwaka et al. 2017).

References: Ammannito et al. (2017) Science 353 , aaf4279 . Aponte et al. (2017) ACS Earth & Space Chem 1 , 3 . • Bibang et al. (2019) Low Temp. Phys. 45 , 590 . • Bowling et al. (2020) EPSL 534 , 116069 . • Brunetto et al. (2018) PSS 158 , 38 . • Campins et al. (2010) Nature 464 , 1320 . • Castillo-Rogez et al. (2018) MaPS 53 , 1820 . • Castillo-Rogez et al. (2019) GRL 46 , 1963 . • Castillo-Rogez et al. (2020) Astrobiology 20 , 269 . • Chang & Bunch (1986) In Clay Minerals and the Origin of Life , p. 116-129. • Cody et al. (2011) PNAS 108 , 19171 . • Combe et al. (2016) Science 353 , aaf3010 . • Crites et al. (2013) Icarus 226 , 1192 . • De Sanctis et al. (2017) Science 355 , 719 . • De Sanctis et al. (2019) MNRAS 482 :2407 . • De Sanctis et al. (2020) Nat Astron , avail. online Aug. 10 • Elsila et al. (2012) MaPS 47 , 1517 . • Engel (2013) Proc SPIE 8865 , 886505 . • Ermakov et al. (2017) JGR Planets 122 , 2267 . • Gontareva et al. (2009) PSS 57 , 441 . • Hayne & Aharonson (2015) JGR Planets 120 , 1567 . • Hesse & Castillo-Rogez (2019) GRL 46 , 1213 . • Kaplan et al. (2018) GRL 45 , 5274 . • Kebukawa et al. (2013) ApJ 771 , 19 . • Kebukawa et al. (2017) Sci Adv 3 , e1602093 . • Kerridge (1985) GCA 49 , 1707 . • Krohn et al. (2018) Icarus 316 , 84 . • Lantz et al. (2015) A&A 577 , A41 . • Larson et al. (1983) Icarus 56 , 398 . • Lemaitre & Morbidelli (1994) Cel Mech Dyn Astron 60 , 29 . • Marchi et al. (2019) Nat Astron 3 , 140 . • Nathues et al. (2017) AJ 154 , 84 . • Neesemann et al. (2019) Icarus 320 , 60 . • Neveu & Desch (2015) GRL 42 , 10197 . • Neveu et al. (2017) GCA 212 , 324 . • Pavlov et al. (2012) GRL 39 , L13202 . • Pavlov et al. (2014) JGR Planets 119 , 1390 . • Postberg et al. (2018) In Enceladus and the Icy Moons of Saturn , p.129-162 . • Rivkin & Emery (2010) Nature 464 , 1322 . • Rivkin et al. (2019) JGR Planets 124 :1393 . • Ruesch et al. (2019) Nat Geo 12 , 505 . • Scully et al. (2019) Icarus 320 , 213 . • Shock & Schulte (1990) Nature 343 , 728 . • Simakov & Kuzicheva (2005) Adv Sp Res 36 , 190 . • Takano et al. (2007) EPSL 254 , 106 . • Travis et al. (2018) MaPS 53 , 2008 . • Vernazza et al. (2017) AJ 153 , 72 . • Vernazza et al. (2020) Nat Astron 4 , 136 . • Vinogradoff et al. (2018) Icarus 305 , 358 . • Vinogradoff et al. (2019) EPSC-DPS abstract 690-1 . • Vinogradoff et al. (2020) GCA 269 , 150 . • Zappalà et al. (1995) Icarus 116 , 291 .