Tianling Lyu

Wei Zhao

Yinsu Zhu

Zhan Wu

Yikun Zhang

Yang Chen

Limin Luo

Shuo Li

Lei Xing

Estimating dual-energy CT imaging from single-energy CT data with material decomposition convolutional neural network

Keywords: Dual-energy CT, deep learning, material decomposition, convolutional neural network

Dual-energy computed tomography (DECT) is of great significance for clinical practice due to its huge potential to provide material-specific information. However, DECT scanners are usually more expensive than standard singleenergy CT (SECT) scanners and thus are less accessible to undeveloped regions. In this paper, we show that the energy-domain correlation and anatomical consistency between standard DECT images can be harnessed by a deep learning model to provide high-performance DECT imaging from fully-sampled low-energy data together with singleview high-energy data. We demonstrate the feasibility of the approach with two independent cohorts (the first cohort including contrast-enhanced DECT scans of 5,753 image slices from 22 patients and the second cohort including spectral CT scans without contrast injection of 2463 image slices from other 22 patients) and show its superior performance on DECT applications. The deep-learning-based approach could be useful to further significantly reduce the radiation dose of current premium DECT scanners and has the potential to simplify the hardware of DECT imaging systems and to enable DECT imaging using standard SECT scanners.

Introduction

Material differentiation and quantification using a standard single-energy computed tomography (SECT) are extremely challenging because different materials may have the same CT value [START_REF] Mccollough | Dual-and multi-energy ct: principles, technical approaches, and clinical applications[END_REF]. To tackle this challenge, dual-energy CT (DECT) takes full advantage of the energy dependence of the linear attenuation coefficient by scanning the patients using two different energy spectra (Alvarez and Macovski, 1976a[START_REF] Kalender | Evaluation of a prototype dual-energy computed tomographic apparatus. i. phantom studies[END_REF][START_REF] Flohr | First performance evaluation of a dual-source ct (dsct) system[END_REF][START_REF] Johnson | Material differentiation by dual energy ct: initial experience[END_REF][START_REF] Boll | Calcified vascular plaque specimens: assessment with cardiac dual-energy multidetector ct in anthropomorphically moving heart phantom[END_REF][START_REF] Lee | A feasibility study of low-dose single-scan dual-energy cone-beam ct in many-view under-sampling framework[END_REF][START_REF] Petrongolo | Single-scan dual-energy ct using primary modulation[END_REF][START_REF] Xue | Accurate multi-material decomposition in dual-energy ct: A phantom study[END_REF]. This enables DECT imaging providing ACCEPTED MANUSCRIPT energy-and material-selective images, and having been very widely used in clinical practice for many applications, such as virtual monochromatic imaging [START_REF] Yu | Dualenergy ct-based monochromatic imaging[END_REF][START_REF] Pomerantz | Virtual monochromatic reconstruction of dual-energy unenhanced head ct at 65-75 kev maximizes image quality compared with conventional polychromatic ct[END_REF], differentiating intracerebral hemorrhage from iodinated contrast [START_REF] Phan | Differentiation of hemorrhage from iodinated contrast in different intracranial compartments using dual-energy head ct[END_REF], automated bone removal in CT angiography [START_REF] Sommer | The value of dual-energy bone removal in maximum intensity projections of lower extremity computed tomography angiography[END_REF][START_REF] Buerke | Dual-energy cta with bone removal for transcranial arteries: intraindividual comparison with standard cta without bone removal and tofmra[END_REF][START_REF] Morhard | Cervical and cranial computed tomographic angiography with automated bone removal: dual energy computed tomography versus standard computed tomography[END_REF][START_REF] Schulz | Automatic bone removal technique in whole-body dual-energy ct angiography: performance and image quality[END_REF], virtual non-contrast-enhanced imaging [START_REF] Takahashi | Dual-energy ct iodinesubtraction virtual unenhanced technique to detect urinary stones in an iodine-filled collecting system: a phantom study[END_REF][START_REF] Ferda | The assessment of intracranial bleeding with virtual unenhanced imaging by means of dualenergy ct angiography[END_REF], Graser et al., 2009[START_REF] Ho | Characterization of adrenal nodules with dual-energy ct: can virtual unenhanced attenuation values replace true unenhanced attenuation values?[END_REF][START_REF] Mangold | Virtual nonenhanced dual-energy ct urography with tin-filter technology: determinants of detection of urinary calculi in the renal collecting system[END_REF][START_REF] Toepker | Virtual non-contrast in secondgeneration, dual-energy computed tomography: reliability of attenuation values[END_REF] and urinary stone characterization [START_REF] Primak | Noninvasive differentiation of uric acid versus non-uric acid kidney stones using dualenergy ct[END_REF][START_REF] Boll | Renal stone assessment with dualenergy multidetector ct and advanced postprocessing techniques: improved characterization of renal stone composition-pilot study[END_REF][START_REF] Ascenti | Stone-targeted dual-energy ct: a new diagnostic approach to urinary calculosis[END_REF], Leng et al., 2015b). However, it is still an open and challenging task for clinical DECT imaging due to complex practical implementations, proprietary patents for major CT vendors, and lower popularity for DECT scanners compared to the standard SECT scanners.

Since the low-and high-energy CT images acquired from the DECT scanners have the same anatomical structures, there is substantial redundant anatomical information between the DECT images. For the scanned patients using the same DECT imaging protocols, the lowand high-energy CT images are also correlated in the energy-domain, resulting in information redundancies in the energy-domain (Leng et al., 2015a[START_REF] Zhao | Using edge-preserving algorithm with non-local mean for significantly improved image-domain material decomposition in dual-energy ct[END_REF]. Meanwhile, both DECT images are reconstructed using fully-sampled projection data which have to meet the classical Shannon-Nyquist theorem in angular-data sampling to reconstruct artifacts-free images. By fully exploiting the anatomical consistency and energy-domain correlation between the DECT images, it is possible to provide high-quality artifacts-free DECT images using conventional SECT images together with sparse sampling projection data at different energy levels.

CT imaging with the full use of as low as reasonably achievable (ALARA) principle has been commonly accepted in routine practice and further reducing radiation dose from CT scanning is clinically favorable and has been extensively studied for almost two decades. Deep learning (DL) has recently been proved to be a powerful tool for mapping complex relationships and incorporating existing knowledge into an inference model through feature extraction and representation learning [START_REF] Gulshan | Development and validation of a deep learning algorithm for detection of diabetic retinopathy in retinal fundus photographs[END_REF][START_REF] Esteva | Dermatologist-level classification of skin cancer with deep neural networks[END_REF][START_REF] Ting | Development and validation of a deep learning system for diabetic retinopathy and related eye diseases using retinal images from multiethnic populations with diabetes[END_REF][START_REF] Liu | Deep learning mr imaging-based attenuation correction for pet/mr imaging[END_REF], Xing et al., 2018, Zhao et al., 2019a[START_REF] Maier | Learning with known operators reduces maximum error bounds[END_REF][START_REF] Shen | Patient-specific reconstruction of volumetric computed tomography images from a single projection view via deep learning[END_REF][START_REF] Lee | Development of a deep neural network for generating synthetic dual-energy chest x-ray images with single x-ray exposure[END_REF], Zhao et al., 2019b). It has also been applied in low-dose CT [START_REF] Chen | Low-dose ct with a residual encoderdecoder convolutional neural network[END_REF][START_REF] Yang | Low-dose ct image denoising using a generative adversarial network with wasserstein distance and perceptual loss[END_REF][START_REF] Wolterink | Generative adversarial networks for noise reduction in low-dose ct[END_REF], Kang et al., 2018) and DECT imaging [START_REF] Liao | Pseudo dual energy ct imaging using deep learning-based framework: basic material estimation[END_REF], Zhang et al., 2019[START_REF] Poirot | physics-informed deep learning for dual-energy computed tomography image processing[END_REF][START_REF] Feng | Fully connected neural network for virtual monochromatic imaging in spectral computed tomography[END_REF][START_REF] Zhao | Obtaining dual-energy computed tomography (ct) information from a single-energy ct image for quantitative imaging analysis of living subjects by using deep learning[END_REF].

To reduce the radiation dose of DECT imaging, in this study, we synergically exploit the energy-domain correlation and anatomical consistency between DECT images by leveraging the deep learning approach and the seamless integration of the correlation and consistency in a data-driven DECT imaging process, and eventually push the sparse sampling to the limit of a single projection view and demonstrate the feasibility of high-performance DECT imaging using a deep learning approach termed fully low-energy and single high-energy DECT imaging (FLESH-DECT).

Material and methods

The flowchart of the proposed FLESH-DECT strategy is shown in Fig. 1. The input to the model is a single-view high-energy projection together with the low-energy image I low which is reconstructed using fully-sampled lowenergy projection data. For the low-energy image, it is firstly denoised with a denoising network to mitigate the impact of image noise. Instead of training a network directly mapping high-energy images from the low-energy images, we use a convolutional neural network (CNN) to perform material-decomposition-type operations and it is termed material decomposition CNN (MD-CNN). The input of the MD-CNN is the denoised low-energy image I de low while the output is a "material component" matrix A. The matrix A has the same image size as I de low but with multiple channels each of which corresponds to a pseudo-material-specific image. The values are the percentages of corresponding "basis material" on these pixels. We ensure that the sum of the percentages equals one (mass conversation) for each unique pixel. Furthermore, another CNN is used to pre-process the differences between the given high-energy projections and their corresponding low-energy projections. This projection-domain network is used to fill the gap between the denoised images and the non-denoised projection. Since CT forward projection can be regarded as linear summations of pixel values, we can use the least-squares method to solve the corresponding CT values of each "material" b di f according to the matrix A and the pre-processed projection difference. The estimated high-energy image is calculated as the summation of low-energy images, and the inner product of matrix A and vector b di f . Detailed formula derivation is described in the following subsections.

DECT Imaging

In CT imaging, the attenuation coefficient at each position can be represented as a linear combination of basis materials' attenuation coefficient (Alvarez and Macovski, 1976b).

µ = α 1 µ 1 + α 2 µ 2 + • • • + α m µ m (1)
where m is the number of basis materials, α i is the percentage of the i-th basis material and µ i is the attenuation coefficient of the i-th basis material. Since CT values in Hounsfield Unit (HU) can be represented as the linear transformation of attenuation coefficient µ with the following equation:

HU = 1000 × µ -µ water µ water (2)
Eq.( 1) can also be written as:

HU = α 1 HU 1 + α 2 HU 2 + • • • + α m HU m (3)
where HU i stands for the Hounsfield Unit CT value for the i-th basis material. Considering there are n pix = W×H pixels in an image slice, we can write Eq.(3) into the following matrix multiplication form

I = A • b T (4)
where I is the image vector sized n pix ×1 containing CT values at each pixel, A = [α i j ] n pix ×m is the material component matrix, α i j stands for the percentage of material j at pixel i, b = [HU i ] 1×m consists of HU values for each material.

In DECT, there are two different images I low and I high . The material component matrix A remains the same for both images because pixel compositions do not change between low-and high-energy scans. Therefore, we have the following equations

       I low = A • b T low I high = A • b T high (5)
By subtracting the high-energy equation from the lowenergy equation, we get

I di f = A • b T di f (6)
where I di f is the difference image between I low and I high , b di f is the difference between b low and b high . Let P high and P low be the given high-and low-energy projection measurements, and R be the projection matrix sized n ray ×n pix corresponding to the high-energy view. We have the following equation

R•A • b T di f = R•I di f = R•I high -R•I low = P high -P low (7)
For the CT imaging task with fully-sampled low-energy and single-view high-energy projection, the unknowns in Eq.( 7) are the material component matrix A and the corresponding difference values b di f . Assuming the material component matrix A is known, let M = R•A ∈ n ray ×m , P di f = P high -P low ∈ n ray ×1 , the difference values b di f can be calculated by solving the equation M•b T = P di f . In regular CT imaging, we have n ray >> m, the best b di f can therefore be found using least-squares method which can be computed with Cholesky decomposition, i.e.,

b di f = argmin b M• bT -P di f 2 2 = [(M T M) -1 M T P di f ] T (8)
The only task now is to estimate the material component matrix A from the low energy image I low .

Material decomposition-based dual-energy CT mapping

Due to its ability to learn complex relationships and incorporate existing knowledge into a nonlinear mapping model, a dedicated CNN model (termed MD-CNN) is used to estimate the material component matrix A. When designing the MD-CNN model, a major challenge is the lack of training labels due to unknown materials in the images and their percentages. To tackle this challenge, we train the MD-CNN indirectly. We firstly denoise the dualenergy image pairs, and the denoised low-energy images I de low are inputted into the MD-CNN to acquire material component matrices A DL . Meanwhile, we put the projection differences into another 1-D projection domain CNN to preprocess the projection data. Then, we compute b di f using Eq.( 8). The estimated denoised high-energy images I dl high can therefore be calculated as

I dl high = I de low + A DL •b T di f (9)
An image similarity loss is calculated between the denoised high-energy image I de high and the DL-estimated image I de DL , and the mean squared error (MSE) loss is used for the task:

L high = 1 n I de high -I de DL 2 2 (10)
Instead of getting material component labels, we focus on the target high-energy image and it is not necessary to specify each channel in A DL to represent real material or a linear combination of different materials. Meanwhile, since matrix A DL is supposed to be the material component matrix, it should be able to recover the input denoised low-energy image as well, resulting in the following loss function:

L low = 1 n I de low -A DL •b T low 2 2 (11)
The same strategy in Eq.( 8) is used to calculate the HU values b low for each "material" under low energy,

b low = argmin b A DL • bT -I de low 2 2 (12)
The final loss function is computed as the summation of L low and L high , i.e.,

L = L low + L high (13)
During the inference phase, the low-energy images are also firstly denoised and then inputted into the trained MD-CNN for A DL . The projection domain CNN preprocesses the projection difference P di f . The estimated difference images are calculated according to Eq.( 8). The difference between inference and training is that the estimated high-energy image is calculated as the summation of the original low-energy image I low and the difference image I di f in the inference phase.

Network details

Denoising CNN

We employed the denoising network in our previous work [START_REF] Zhao | Obtaining dual-energy computed tomography (ct) information from a single-energy ct image for quantitative imaging analysis of living subjects by using deep learning[END_REF] to reduce the DECT image noise. The network uses a plain structure that encompasses 13 convolution layers to learn the residual between the input image and the denoised image. The first 12 layers are convolution layers with kernel size 3 × 3. Each convolution layer is followed by a batch normalization layer (BN) and a rectified linear unit (ReLU) activation. The last layer is a convolution layer with kernel size 1×1 fusing the result. Fig. 2 shows the detailed structure of the denoising CNN. The denoised image is computed as the summation of the input image and the output from the last layer.

MD-CNN

For the material decomposition network, we employ a U-Net-type structure [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF] which has a large receptive field and is quite suitable for many medical image processing tasks. There are 10 normal 3 × 3 convolution layers in the proposed network (Fig. 3). Each convolution layer is followed by a BN layer and a ReLU activation layer. There are 3 resolution levels in total. For down-sampling, we use convolution layers with kernel size of 2 × 2 and stride equal to 2. Each strided convolution layer is also followed by a BN layer and a ReLU activation layer. For up-sampling, we use bilinear interpolation to double both image width and height. At the end of the network, a convolution layer with kernel size of 1 × 1 is added to fuse the channels. Since the values on each output channel are supposed to be the percentages of corresponding "basis material", we apply softmax after the convolution layer with kernel size of 1 × 1 to make sure that the sum of materials' proportion equals to 1 at each pixel.

Projection-domain CNN

To enhance the robustness of the least-squares problem in Eq.( 8), the projection-domain network (Fig. 4) is applied to slightly refine the inputting projection difference and to make its noise level to be consistent with that of the denoised DECT difference image, and we employ a concise 4-layer network for the task. A residual learning structure is used here for an easier startup at the beginning iterations. All basic convolution layers have a kernel with a size of 1 × 5 except for the last one which has a kernel with a size of 1 × 1. The first three convolution layers are followed by a BN layer and a ReLU activation layer.

Network training

The denoising CNN was implemented in MATLAB with MatConvNet framework [START_REF] Vedaldi | Matconvnet: Convolutional neural networks for matlab[END_REF]. Denoising was performed as a data preprocessing step for all DECT images. The material decomposition network and projection-domain network are implemented using Python with Tensorflow framework [START_REF] Abadi | Tensorflow: Large-scale machine learning on heterogeneous distributed systems[END_REF]. Those two networks were trained together in an end-to- end fashion. The parameters in the networks were optimized using ADAM algorithm [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF] with β 1 = 0.9 and β 2 = 0.999. The learning rate was set to 10 -3 during the training. The training set was randomly split into small batches in each epoch with a batch size of 8. The proposed network was trained for 200 epochs in total. We perform validation after each training epoch and the model with the best validation loss was selected as the final model for testing. The number of materials m was set to 10 in our experiments according to the loss on the validation set. We trained models with different m settings from 2 to 20 and the best validation loss came with m = 10. It has to note that the model was able to achieve reasonable results even with m = 3. As shown in Fig. 5, the validation losses are close when m > 2. The networks were trained and tested on a workstation with configurations as follows: CPU is Intel(R) Xeon(R) Gold 6130 CPU @ 2.10GHz; GPU is NVIDIA RTX 2080 Ti with 12 GB memory.

Projectors for different CT geometries

To calculate matrix M in Eq.( 8), a projector corresponding to the CT geometry is indispensable in our algorithm. There are mainly two types of 2D CT geometries, fan-beam and parallel-beam. The fan-beam geometry can be further divided into two sub-types, equiangular and equispaced. For each above-mentioned geometry, we developed a projector and trained a new model for evaluation.

Equiangular fan-beam

Equiangular geometry is mainly implemented with arc detectors which keep the angles between two adjacent detector pixels and the source-detector-distance the same. We tested our model in the anterior-posterior (AP) direction, but it can be easily extended to any other projection view. To acquire projection data using the 2D equiangular fan-beam geometry, we first calculate the intersection of each projection ray with each image row. The values at each intersection point are computed using linear interpolation. Suppose the image size is W×H, and a matrix I with a size of n ray ×H can be obtained using the following rebinning equation:

I (φ, y) = I((D -y) tan(φ), y) ( 14
)
where φ is the angle between the projection ray and the central ray (as shown in Fig. 6), and y is the image pixel position along the y-axis in the Cartesian coordinate system centered at the image center O. D is the distance between the projection source and the rotation center which is overlapped with the image center O. After the linear interpolation, we calculate the summation of each column in matrix I to obtain ray sum S which is weighted by the distance to yield the final projection P:

P(φ) = dy cos φ S (φ) ( 15
)
where dy is the image spacing in y-axis. We set D = 600mm, dy = 0.5mm, the number of detector channels n ray = 800 and the angle between adjacent channels ds = 7 11000 rad for all images.

equispaced fan-beam

Flat detectors with equal spacing between the adjacent channels are commonly used to implement the equispaced geometry. We used a similar strategy to implement the equispaced fan-beam projection in the AP direction. In this case, the rebinning procedure is computed with the following equation:

I (s, y) = I( s(D -y) L , y) ( 16
)
where s is the distance between the current channel and the detector center (direction included), and L is the source to detector distance, as shown in Fig. 7. The rest part of the implementation is the same as the equiangular projection. For the parameters, we have D = 600mm, L = 1100mm, dy = 0.5mm, n ray = 800 and the spacing between nearby channels ds = 0.78mm.

Parallel-beam

For the parallel-beam scenario, We assumed that the detector channel has the same spacing as the image pixels, and each X-ray projects exactly through an image column in the AP direction. Therefore, the projection in the AP direction can be computed as the summation of each image column.

Dataset

Training data for the denoising network

The AAPM Low-Dose CT Grand Challenge data was used to train the denoising network. This dataset consists of routine dose CT and the corresponding simulated low-dose CT data from 10 patients. The routine dose scanning voltage is 100 kV or 120 kV and the Xray tube current varies from 200mA to 500mA. The detector has 736 × 64 elements, and each element has a size of 1.2856 × 1.0947mm 2 . The source-to-axial distance is 59.5cm and the source-to-detector distance was 108.56cm. All the images were reconstructed to slice thickness of 1.0mm and 512 × 512 pixel. The pixel size varies from 0.66 × 0.66mm 2 to 0.78 × 0.78mm 2 . To simulate the low-dose CT data, Poisson noise was introduced into the routine dose to mimic a noise level that corresponded to 25% of the routine dose, and noisy projection was reconstructed to yield the low-dose image.

Siemens DECT image dataset

Clinical DECT images of 22 patients who underwent iodine contrast-enhanced DECT exams were collected for the study. All the exams were performed in Nanjing General PLA Hospital, China, with the approval of the institutional review board and patient consent forms. The DECT images (5753 slices in total) were acquired using a SOMATOM Definition Flash DECT scanner (Siemens Healthineers, Forchheim, Germany) after administering iodine contrast agent. The low-and high-energy of the DECT scans were 100 kV and 140 kV, respectively. All CT images were reconstructed using the filtered backprojection (FBP) algorithm provided by the commercial CT vendor. The dataset was split into a training set, a validation set and a testing set randomly with 16, 3, and 3 patients included respectively.

ACCEPTED MANUSCRIPT

Philips spectral CT dataset

We collected clinical spectral CT images of 22 patients with Philips IQon Spectral CT devices. All exams were performed in Nanjing Gulou Hospital, China, with the approval of the institutional review board and patient consent forms. Philips IQon Spectral CT applies spectral detector technology to implement dual-energy scans. Virtual mono-energetic images at different X-ray energy levels can be reconstructed from the data collected. In our experiments, images at 70 keV, 110 keV, and 150 keV were used to evaluate the proposed method. Noncontrast-enhanced brain virtual mono-energetic images (2463 slices in total) were collected for each energy level and all images were reconstructed using a build-in algorithm provided by the CT vendor. Similar to the Siemens dataset, the Philips dataset was also split into a training set, a validation set and a testing set randomly with 16, 3 and 3 patients included respectively.

Results

We firstly focus on the results of equiangular geometry on the Siemens Dataset, and comparison results using different geometries are presented in the middle of this section. Results on the Philips dataset are described at last.

Fig. 8 shows the original DECT images and the 140 kV images predicted using the proposed method for a testing patient. The first, second and third columns show the original 100 kV images, the original 140 kV images, and the predicted results, respectively. The first, third, and fifth rows show CT images in transverse, sagittal, and coronal planes, respectively. The second, fourth, and sixth rows show difference images with respect to the corresponding real high-energy CT images in transverse, sagittal, and coronal planes, respectively. As can be seen, the proposed DL-derived high-energy images are highly consistent with the original high-energy images. There are some differences at sharp boundaries which also appear in the difference images between original high-and low-images. Those differences may be motion introduced differences between original DECT images because there are approximately 90 degrees out of phase for the low-and highenergy data acquired using a dual-source DECT scanner. When inputting the low-energy image into the model, Figure 8: Example results on a testing slice and the difference with respect to the corresponding real high-energy image. The first, third and fifth rows display the images on axial, sagittal and coronal view. From left to right are real 100 kV image, real 140 kV image and result from the proposed method, respectively. A kidney lesion in transverse image is labeled using red arrow. The second, fourth and sixth rows are the corresponding differences with real 140 kV images. The CT images are displayed with a window width=300 HU and center=50 HU while the difference images are displayed under window width=300 HU and center=0 HU.

the model performed prediction based on the anatomical structure of the low image and can not reflect the change with respect to the original high-energy image.

Quantitative metrics were calculated to evaluate the accuracy of the predicted high-energy CT images. We used the well-established metrics MSE, PSNR and SSIM to as- sess the image similarity between the real and the predicted high-energy images. Additionally, more than 100 region-of-interest (ROIs) were randomly selected for each testing volume on homogeneous areas (e.g. liver and stomach). We calculated the mean HU value differences on those ROIs and the results show that the averaged HU error between the predicted and the original high-energy image is smaller than 2.09 HU. For the computation time, it took around 2.5 seconds for the proposed method to process 300 slices on our server. All quantitative results are shown in Table 1. Since most of the differences come from noise, we also compared the denoised DL-predicted images with the denoised high-energy images. In this case, the DL-predicted images were calculated by adding the denoised low-energy image to the DL-estimated difference image. For those denoised images, the proposed method achieved an average MSE of 171.09, PSNR of 44.33 and SSIM of 0.9848.

MANUSCRIPT

In our previous work [START_REF] Zhao | Obtaining dual-energy computed tomography (ct) information from a single-energy ct image for quantitative imaging analysis of living subjects by using deep learning[END_REF], we proposed the DL-DECT method to estimate I di f directly from I low . In this work, with the introduction of the additional high-energy single-view projection, we can further enhance the accuracy of the predicted high-energy image and eventually the accuracy of the material-and energy-specific images. To show the benefit of the additional high-energy projection, we quantitatively compared the proposed FLESH-DECT method with the DL-DECT method. Virtual non-contrast (VNC) images and iodine maps were derived to demonstrate the clinical utility of the proposed method. Fig. 9 depicts the VNC images and iodine maps reconstructed using different methods on the transverse, sagittal and coronal planes. Both the DL-DECT and FLESH-DECT algorithms provided highquality VNC and iodine images that were consistent with the images generated by original DECT images. Quan- titative metrics on VNC and iodine images are shown in Table 2. The results demonstrate FLESH-DECT can provide high-quality material-specific images and it outperforms the DL-DECT method.

From the VNC images and Iodine maps, we found that the DL-derived VNC images and iodine maps show a remarkably reduced noise level compared with those generated from original DECT images. Here we also compared the noise level by calculating the standard deviation 3. The mean standard deviations of the images derived using FLESH-DECT are close to that of the DL-DECT method which is much lower than those from the original DECT images. This result shows that the FLESH-DECT method maintains the denoising feature.

ACCEPTED MANUSCRIPT

We tested the proposed method with several 2D CT geometries. Example results on a testing volume using different geometries are displayed in Fig. 10. All CT images are displayed with window width=300HU and center=50HU while difference images are displayed with window width=300HU and center=0HU. As can be seen, all models can provide competitive results that are highly consistent with the original high-energy image. For better comparison, we also calculated quantitative metrics on those results (Table 4). The differences between the results obtained using the proposed method with different geometries are marginal and all of them are superior to our previous DL-DECT method which does not utilize the additional single-view projection. Overall, the proposed method reduces mean-squared error (averaged for all testing cases) from 1858.32 to 898.35 while increasing PSNR from 33.83 to 36.89 and SSIM from 0.8641 to 0.8744. The computation time for different methods are displayed and compared in Table 5. Note that computation time for the denoising CNN is not included in the results. The denoising network takes about 0.01 seconds for each slice and the denoising time is the same for all methods. Compared to the DL-DECT method, the proposed method speeds-up the computation time by 2-fold which can be MANUSCRIPT attributed to the reduced number of weights and the simple network structure. There are some differences among the time using different geometries which means the computational cost of the proposed method depends on the projector.

To further examine the proposed method, crossvalidation on the Siemens dataset was also derived. Seven-fold cross-validation was used with 3 testing cases for each fold (except for the last one with 4 cases). The results are displayed in Table 6. In cross-validation, the average mean-squared-error value of all 22 cases achieves 823.63 while average PSNR achieves 37.36 and average SSIM achieves 0.8862. All those quantitative metrics are close to the results on our testing set which shows the robustness of our proposed method. MSE values in each cross-validation group are also analyzed and the results are depicted in Fig. 11. For each group, the average MSE value lies between 700 and 1000. Clinical doctors were invited to evaluate the artificially generated images. Three experienced radiologists from Jiangsu People's Hospital were involved in this process. The radiologists were supposed to provide a score from 1 to 5 (1 for useless, 5 for high clinical use) to real 100 kV, 140 kV images from practical CT machine as well as artificial 140 kV images generated using the proposed method. Images of all 22 patients (together with all generated results in cross-validation) were evaluated. Fig. 12 shows the clinical evaluation results. The average score of the artificially generated 140 kV images is the same as that of real 140 kV images because the generated images look almost the same as the real ones. We also tested the proposed method on the Philips spectral CT image dataset. One independent model was trained between each energy level pair (e.g. between 70 keV and 110keV). Therefore, we got 3 models for the Philips dataset in total. For each model, 1828 slices were used as training data while 299 slices were used for validation and 336 for testing. All the hyper-parameters were set the same as those selected on the Siemens dataset. One difference that should be noticed is that we did not per-MANUSCRIPT form image denoising on Philips dataset because of its low noise level. The results are concluded in Table 7. For the task of estimating 150 keV images from 70 keV images, the average MSE between our results and the target images is only 13.45. Average PSNR value and SSIM value achieve 55.84 and 0.9985, respectively. The models for the other two tasks (estimating 110 keV images from 70keV images and 150 keV from 110 keV) perform even better on those metrics since the differences between real low-and high-energy images are smaller. Low noise level plays an important roll in the improvement in quantitative measurements. Besides, anatomical structures in non-enhanced brain images are much simpler than those in abdominal images. Results on one example slice are displayed in Fig. 13. Brain tissue has similar CT values under different X-ray energy levels and the differences mainly come from the bone area. The proposed method performs well on the new Philips dataset acquired using a different DECT scanning protocol.

ACCEPTED MANUSCRIPT

Discussion

There are substantial redundant information and correlation in both anatomical structure and energy-domain between the low-and high-energy DECT images. By incorporating the redundancy and the correlation into a deep learning model, it is possible to provide materialand energy-specific images using standard SECT scanners, which has the potential to alleviate the need for premium DECT scanners. Besides, compared to the standard fully low-and high-energy sampling DECT mechanism, the use of sparse sampling at the second energy level can significantly reduce the radiation dose of DECT imaging.

Our results show both the DL-DECT and FLESH-DECT methods achieve high-performance DECT imaging by using the input low-energy CT data, and quantitative analysis shows FLESH-DECT outperforms DL-DECT in terms of HU accuracy and calculation speed. The superior performance of the FLESH-DECT can be attributed to the additional single-view high-energy projection. Different from the DL-DECT method which directly infers a high-energy image using the incorporated prior knowledge, the proposed FLESH-DECT method uses the learned knowledge to fit the measured high-energy projection. Namely, the high-energy projection introduces a penalty to constrain the projection generated by the predicted high-energy images to be consistent with the measurement, which in turn enhances the accuracy of the predicted images. The single-view high-energy projection can be obtained shortly before or after the standard SECT low-energy data acquisition, and existing SECT systems are able to implement these scanning protocols without modifying the hardware.

In the proposed method, the single-view high-energy projection and the corresponding low-energy projection view are supposed to be strictly aligned. We moved the X-ray source under high-energy to mimic patient movements and mismatches between projection source positions. The results are depicted in Fig. 14. The model is robust to shift in axial direction to some extent. It outperforms DL-DECT model even when there is a 10mmshift. The model is also robust concerning the movements in horizontal directions (projection source moving left or right). Changes in distance between patient and projection source are not accepted as they introduce severe mismatch on organ boundaries in projections. The proposed method is unlikely to introduce artifacts when there is mild mismatch between projections at different energylevels. The material maps (in which structural information is included) are generated according to the given lowenergy image. The anatomical structure will be well preserved when there is no artifact in the low-energy image. The projection difference only affects the CT values of each "material".

It should be noticed that the high-energy projection in the proposed method does not need to cover the whole area the corresponding low-energy projection scanned. Since in most cases, the number of equations in Eq.7 overwhelms the number of unknowns, reducing the number of equations (which means reducing the number of projections) has little influence on the results. From our experiments, MSE between the real image and the generated image will not change severely if we cut the projection channel number into half or even quarter, though full projections do perform slightly better (see Fig. 15). FLESH-DECT is suitable for different geometries. In this study, we have tested the method using 2D geometries (fan-beam and parallel beam). However, the method can be applied straightforwardly to 3D geometries by extending the networks and the projection operators into 3D scenarios. Since the reconstructed size is smaller than the ACCEPTED MANUSCRIPT projection view size in the rotation axis for a 3D case, projections generated using the reconstructed image volume cannot match the measured projection. To solve this issue, one can use the middle part of the measured projection which does not propagate through the region beyond the image volume.

The proposed method employs the MD-CNN to generate "material-decomposition" maps A from 100 kV images. Since basis materials in the image do not change under X-ray at different energy levels, accurate materialdecomposition maps can provide CT images under any spectrum if combined with the projection view acquired using the corresponding spectrum (e.g. 120 kV image from 120 kV projection, 150 kV image from 150 kV projection). In our experiments, the models were trained and tested using DECT images under 100 kV/Sn 140 kV scanning protocol. Therefore, the generated "materialdecomposition" maps are likely to be optimized for this specific protocol and may not be applied to images acquired using a different protocol (such as 80 kV/Sn 140 kV). An example is shown in Fig. 16. However, if the model was trained using images acquired from several different spectra, the "material-decomposition" maps would be much closer to the real ones and the model would have the potential to generate images under different spectra without re-training the MD-CNN. Also, regularizers may further be introduced and applied to matrix A during network training to enhance the robustness of the model.

There is a denoising-CNN included in the flowchart of FLESH-DECT. We use this network to reduce the impact of image noise and the results show its effectiveness. However, the denoising network is not mandatory and it can be replaced with other image denoising techniques [START_REF] Ma | Low-dose computed tomography image restoration using previous normal-dose scan[END_REF], such as non-local mean (NLM) [START_REF] Zhang | Iterative reconstruction for xray computed tomography using prior-image induced nonlocal regularization[END_REF], block-matching and 3D filtering (BM3D) [START_REF] Salehjahromi | A spectral ct denoising algorithm based on weighted block matching 3d filtering[END_REF], without see-ing severe degradation in performance. The denoising step may also be removed when inputting extremely-high quality images.

In this paper, the proposed method was tested on the Siemens DECT dataset as well as Philips spectral CT dataset. For Siemens DECT, both low-and high-energy projections can be attained from a CT system while in Philips spectral CT dataset, the projection views of those virtual monoenergetic images are unable to acquire from real CT machines. Here, we are using those simulated data to demonstrate that the proposed method has the potential to deal with different energy-level settings.

Despite all the advantages and potentials mentioned above, the proposed method has limitations. FLESH-DECT relies highly on the deep neural network to perform the material-decomposition-like operation. Since the domain knowledge learned by the network highly depends on training data, it is unlikely to provide reasonable results when there is a huge difference between training and testing data. For example, models trained under 100 kV/Sn 140 kV protocol may not generate correct results when inputting low-energy images scanned under 70 kV protocol. However, these limitations can be solved by training different models for different protocols.

Conclusion

In this paper, we proposed a deep learning approach to perform DECT imaging using a low-energy image and a single-view high-energy projection. Compared to standard DECT imaging, the approach can provide superior material-specific images with significantly reduced noise. It also has the potential to simplify the system design and reduce the radiation dose and allows us to perform high-quality DECT imaging without the conventional hardware-based DECT solutions. The approach may significantly extend the usage of the widespread standard SECT scanners by providing advanced DECT clinical applications, such as urinary stone characterization and differentiating intracerebral hemorrhage from iodinated contrast, and thus lead to a new paradigm of SECT imaging.

Figure 1 :

 1 Figure 1: The workflow of the proposed fully-sampled low-energy and single-view high-energy DECT imaging approach. During the training phase, the denoised DECT images together with the single-view dual-energy projections are used to train the projection domain convolutional neural network (CNN) and the MD-CNN. In the testing phase, the trained networks use the input low-energy images and the single-view dualenergy projections to infer the corresponding high-energy images.

Figure 2 :

 2 Figure 2: The architecture of the fully convolutional network for image denoising. A plain structure that encompasses 13 convolution layers is applied to learn the residual between the input image and the denoised image.

Figure 3 :

 3 Figure 3: Structure of the MD-CNN. A much simplified UNet-like structure with 14 convolution layers is used here to estimate the percentages of each corresponding "basis material" from the denoised low-energy image. Numbers under each block show the channel number of the multichannel feature maps.

Figure 4 :

 4 Figure 4: Structure of the projection-domain CNN. A concise 4-layer network is employed to slightly refine the inputting projection difference. Numbers under each block show the number of channels of the multichannel feature maps.

Figure 5 :

 5 Figure 5: Validation loss under different m settings. Best validation loss comes with m = 10.

Figure 6 :

 6 Figure 6: Illustration of the equiangular fan-beam geometry.

Figure 7 :

 7 Figure 7: Illustration of the equispaced fan-beam geometry.

Figure 9 :

 9 Figure 9: VNC images and Iodine maps reconstructed using original DECT, DL-DECT, and FLESH-DECT images. The first, third and fifth row are the VNC images on axial, sagittal and coronal view, respectively, while the second, fourth and sixth row are the corresponding Iodine maps. The kidney lesion in the transverse images is marked out with red arrows for better comparison.

Figure 10 :

 10 Figure 10: Results on three testing slices for different geometry models. From left to right are original 100 kV images, original 140 kV images, parallel-beam results, equispaced fan-beam results and equiangular fan-beam results, respectively.

Figure 11 :

 11 Figure 11: Mean MSE (standard deviation) in each crossvalidation group.

Figure 12 :

 12 Figure 12: Clinical evaluation results.

Figure 13 :

 13 Figure13: Example results on a Philips dataset testing slice and the difference with respect to the corresponding high-energy image. The CT images are displayed with a window width=300 HU and center=50 HU while the difference images are displayed under window width=300 HU and center=0 HU.

Figure 14 :

 14 Figure 14: MSE values on test patient1 when there is shift between projections at different energy-levels. (a) Shift in axial direction. (b) Shift in horizontal directions.

Table 1 :

 1 Quantitative comparisons between the predicted and the real high-energy CT images for the testing patients.

	MSE	PSNR SSIM	HU error Time(s)
	Patient1 875.06 36.80 0.8702 1.65±1.03	2.65
	Patient2 887.60 36.99 0.8789 2.09±1.48	2.31
	Patient3 927.47 36.89 0.8745 1.50±0.82	3.34

Table 2 :

 2 Quantitative comparisons the VNC image and the iodine maps reconstructed using original DECT and FLESH-DECT images.

		VNC	Iodine Map
		MSE	PSNR	MSE	PSNR
	Patient1	DL-DECT FLESH-DECT 3808.97 29.11 0.0308 23.73 3882.34 29.03 0.0313 23.65
	Patient2	DL-DECT FLESH-DECT 4128.01 29.04 0.0335 23.87 4216.36 28.95 0.0342 23.78
	Patient3	DL-DECT FLESH-DECT 3431.13 29.81 0.0306 24.95 3508.35 29.71 0.0313 24.85

Table 3 :

 3 Quantitative Noise Level Comparison on VNCs and Iodine Maps.

	Standard Deviation on ROIs VNC Iodine Map
	Patient1	Real DL-DECT	81.58 29.14	0.3363 0.1618
	FLESH-DECT	28.09	0.1499
		Real	74.69	0.2837
	Patient2	DL-DECT	24.81	0.1017
	FLESH-DECT	24.12	0.0998
		Real	76.95	0.2875
	Patient3	DL-DECT	25.97	0.0946
	FLESH-DECT	25.67	0.0891
	in ROIs. More than 500 ROIs were selected randomly
	within homogeneous areas. The mean standard deviations
	in ROIs on each testing patient are provided and com-
	pared in Table

Table 4 :

 4 Quantitative Results for Different Geometries.

			MSE	PSNR SSIM
		DL-DECT	2466.07 32.37 0.8558
	Patient1	Parallel Equiangular 875.06 866.48	36.84 0.8707 36.80 0.8702
		equispaced	871.65	36.82 0.8703
		DL-DECT	1530.75 34.61 0.8723
	Patient2	Parallel Equiangular 887.60 876.68	37.03 0.8793 36.99 0.8789
		equispaced	882.76	37.01 0.8789
		DL-DECT	1578.13 34.61 0.8660
	Patient3	Parallel Equiangular 932.40 927.47	36.89 0.8745 36.88 0.8741
		equispaced	927.15	36.90 0.8747

Table 5 :

 5 Computation Time for Different Methods in Seconds.

		Patient1 Patient2 Patient3
	Slices	308	265	387
	DL-DECT	5.43	4.68	6.86
	Parallel	2.33	2.04	2.94
	Equiangular	2.65	2.31	3.34
	equispaced	2.53	2.18	3.23

Table 6 :

 6 Quantitative Results During Cross-validation.

	MSE 823.62±141.00
	PSNR 37.36±0.7227
	SSIM 0.8862±0.0184

Table 7 :

 7 Quantitative Results Philips Dataset.

	Source(keV) Target(keV)	MSE	PSNR	SSIM
	70	110	13.45 ± 1.94 55.84 ± 0.64 0.9985 ± 1.43 × 10 -4
	70	150	22.27 ± 3.46 53.06 ± 0.66 0.9980 ± 2.65 × 10 -4
	110	150	1.64 ± 0.21 64.22 ± 0.60 0.9999 ± 1.47 × 10 -5

This work was supported by the State's Key Project of Research and Development Plan [grant numbers 2017YFA0104302, 2017YFC0109202 and 2017YFC0107900], by the National Natural Science Foundation [grant 81530060 and 61871117], by Science and Technology Program of Guangdong [grant number 2018B030333001], by NIH grant numbers 1R01CA223667 and R01CA227713] and a Faculty Research Award from Google Inc.


domain: spi phys

text:
Olivier Beltramo-Martin
email: olivier.beltramo-martin@lam.fr

Sam Ragland

Romain Fétick

Carlos Correia

Trent Dupuy

Giulian Fiorentino

Thierry Fusco

Laurent Jolissaint

Sebastian Kamann

Antonino Marasco

D Massari

B Neichel

L Schreiber

P Wizinowich

Scao Slao

Glao Moao

Ltao Mcao

Prime Beltramo-Martin

Airopa -Ciurlo

Martin 2016

Gendron

Ragland

Olguin

Clénet

Véran

Wagner

Fetick

Superstar -Marasco

Pampelmuse -Kamann

Fourier Correia

Neichel

Rigaut

Fusco

Oomao -Correia

COMPASS / Rocket -Ferreira

Dasp -Basden

Passata -Agapito

Blanco & Mugnier

Deconvolution

Review of PSF reconstruction methods and application to post-processing

Keywords:

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

Review of PSF reconstruction methods and application to post-processing

The determination of the point spread function (PSF) is a must in astronomical data processing. Indeed, knowledge of the PSF is critical to estimate the astrophysical signal (photometry/astrometry, extended objects surface topology, or galaxy kinematics, ...) through a post-processing pipeline, that relies on model-fitting or deconvolution for instance. Several issues must be faced to achieve an accurate PSF estimation, in particular field crowding, poor S/N or lack of PSF calibrators, that are strengthened in presence of Adaptive-optics (AO) due to spatio-temporal variability of the PSF. Moreover, the emergence of integral field spectrographs calls for considering 3D tools that also model the PSF spectral diversity. In order to overcome existing limits in the AO PSF determination, a constant effort has been deployed for decades in PSF reconstruction (PSF-R) approaches that aim to estimate the PSF from AO control loop data, e.g. AO wavefront measurements and commands. So far, PSF-R lacks of operational uses and science verification owing to dedicated expertise needed and data flow and storage management. A step further must be achieved to enable PSF-R, especially for Giant Segmented Mirrors Telescopes (GSMT). Our first goal is to provide a standardized categorization of the techniques that use focal-plane or pupil-plane data or both jointly combined (hybrid) so as to gain insight into the advantages and shortcomings of such techniques, with regard to post-processing requirements in terms of PSF knowledge. Secondly, I will summarize the ongoing and upcoming actions that are taken by the community for enabling the implementation of PSF-R facilities.

EXISTING PSF DETERMINATION TECHNIQUES 2.1 Landscape

In the following, we will use the keyword PSF determination as the ensemble of techniques that aim to providing a 2D representation of the PSF from any contextual data. Among them, we distinguish two top-level families as summarized in Fig. 1 :

• Focal-plane based techniques: all methods that use astronomical images to identify the PSF, that is that the input is a PSF or an extended object convolved by a PSF. Extraction and best-fitting methods are the most common techniques, thanks to their flexibility and the absence of AO expertise needed. [START_REF] Diolaiti | StarFinder: A code for stellar field analysis[END_REF][START_REF] Bertin | SExtractor: Software for source extraction[END_REF][START_REF] Stetson | DAOPHOT -A computer program for crowded-field stellar photometry[END_REF] They rely usually on a four-part process that consists of (i) correct frames from dark, backgrounds, flat field and distorsions [START_REF] Bernard | Optimal correction of distortion for high-angular-resolution images: Application to GeMS data[END_REF] (ii) detect stars and extract PSFs from the entire field or sub-fields to account for PSF spatial variations (iii) adjust a parametric model over extracted images or stack them to estimate a PSF function (iv) use the PSF model for science post-processing (deconvolution, photometry/astrometry, ...). However, these methods are sensitive to the image noise and source confusion and require the presence of a point-like source within the scientific field, [START_REF] Massari | GeMS MCAO observations of the Galactic globular cluster NGC 2808: the absolute age[END_REF][START_REF] Yelda | Improving Galactic Center Astrometry by Reducing the Effects of Geometric Distortion[END_REF][START_REF] Schödel | Accurate photometry with adaptive optics in the presence of anisoplanatic effects with a sparsely sampled PSF. The Galactic center as an example of a challenging target for accurate AO photometry[END_REF] calling for consideration of model-fitting instead.

-PSF fitting/interpolation: extracted point sources images and calibrating a PSF model to either deblend sources that overlap or provide a PSF at any position in the field/spectrum.

-PSF marginal estimation and deconvolution: retrieve the PSF from the image of an extended object assuming some priors on the object (and the PSF).

-PSF classification: estimate the PSF from classification methods, such as PCA, supervised classification, and machine-learning/deep-learning -based approaches generally.

• Pupil-plane based techniques : all methods that use wavefront measurements and/or models of wavefront propagation. The input is a time-series of wavefront measurements or integrated parameters derived from wavefront measurements.

-PSF simulation: end-to-end simulations based on integrated atmospheric and system parameters. They allow to include any non-linear effects and rely on wavefront propagation theory. Powerful tools have been developed such as OOMAO, 8 COMPASS, 9 DASP, 10 PASSATA, 11 SOAPY, 12 YAO, 13 OCTOPUS 14 and many others.

-PSF analytical models: From integrated parameters derived from the AO telemetry, one may use analytical models [START_REF] Jolissaint | Synthetic modeling of astronomical closed loop adaptive optics[END_REF][START_REF] Neichel | Tomographic reconstruction for wide-field adaptive optics systems: Fourier domain analysis and fundamental limitations[END_REF][START_REF] Rigaut | Analytical model for Shack-Hartmann-based adaptive optics systems[END_REF] to infer the PSF as well as the PSF spatial variations due to the anisoplanatism for instance [START_REF] Beltramo-Martin | Off-axis point spread function characterization in laser guide star adaptive optics systems[END_REF][START_REF] Britton | The Anisoplanatic Point-Spread Function in Adaptive Optics[END_REF][START_REF] Fried | Limiting Resolution Looking Down Through the Atmosphere[END_REF] -PSF reconstruction: Estimation of the PSF using AO control loop data or integrated parameters estimated from those. PSFR involves several approximations to derive the PSF from the residual phase power spectrum density.

• Hybrid techniques: all methods that uses jointly astronomical images, wavefront measurements and analytical models

Pupil-plane based techniques can only retrieve the spatial frequencies that are effectively measured by the AO wavefront-sensors (WFS). All non-common path aberrations (NCPA) or field-dependent static aberrations must be calibrated independently, which call for dedicated engineering time. On the contrary, focal-plane based techniques permit to recover a noisy representation of the real PSF including atmospheric and instrumental aberrations, but at specific positions in the field and/or for a particular wavelength. In order to be able to characterize the PSF in a whole field or over broadband observation, model of spatial and spectral PSF variabilities are needed. The PSF determination technique that works best the exploitation of a given astronomical image is obviously science-case dependent. We give in the current section an overview on recent approaches that were developed and successfully validated on-sky, which are the analytical PSF model presented in, [START_REF] Beltramo-Martin | Joint estimation of atmospheric and instrumental defects using a parsimonious point spread function model[END_REF][START_REF] Fétick | Physics-based model of the adaptive-optics-corrected point spread function. Applications to the SPHERE/ZIMPOL and MUSE instruments[END_REF] the PSF reconstruction [START_REF] Fusco | Reconstruction of the ground-layer adaptive-optics point spread function for MUSE wide field mode observations[END_REF][START_REF] Ragland | Status of point spread function determination for Keck adaptive optics[END_REF][START_REF] Ragland | Point spread function determination for Keck adaptive optics[END_REF][START_REF] Jolissaint | Adaptive Optics Point Spread Function Reconstruction at W. M. Keck Observatory in Laser Natural Guide Star Modes : Final Developments[END_REF] and the hybrid PSF reconstruction. [START_REF] Beltramo-Martin | PRIME: PSF Reconstruction and Identification for Multiple-source characterization Enhancement -application to Keck NIRC2 imager[END_REF][START_REF] Beltramo-Martin | Pushing Point spread function reconstruction to the next level. Application to SPHERE/ZIMPOL[END_REF]

Scope on advanced PSF models and reconstruction techniques

In advanced PSF models and PSF reconstruction, the PSF h(α) is commonly approximated 29 as a convolution of an instrumental contribution htel with a atmospheric residual term hAO as follows:

h(α) = F -1 htel (ρ/λ) hAO (ρ/λ) (1)
where

F -1 [x] is the inverse 2D Fourier transform of x.
The first term htel in Eq. 1 is obtained by the autocorrelation of the static phase (the amplitude is assumed to be uniform across the pupil) within the telescope pupil. As a consequence, the pupil-plane based approaches requires the calibration of the pupil masks and the static aberrations (e.g., primary mirror cophasing errors, NCPA, field-dependent aberrations for instance), contrary to focal-plane methods for which all these information are encoded into the data. Then, we already see that achieving pupil-plane-based PSF determination is a challenge and requires an expertise on the telescope, AO system and the post-AO instruments.

The calibration of instrumental aberrations, through phase diversity approaches, [START_REF] Jolissaint | Retrieving the telescope and instrument static wavefront aberration with a phase diversity procedure using on-sky adaptive optics corrected images[END_REF][START_REF] Mugnier | On-Line Long-Exposure Phase Diversity: a Powerful Tool for Sensing Quasi-Static Aberrations of Extreme Adaptive Optics Imaging Systems[END_REF][START_REF] Blanc | Calibration of NAOS and CONICA static aberrations. Application of the phase diversity technique[END_REF] focal plane sharpening [START_REF] Lamb | Quantifying telescope phase discontinuities external to adaptive optics systems by use of phase diversity and focal plane sharpening[END_REF] techniques, or even Deep-learning strategies, [START_REF] Andersen | Image-based wavefront sensing for astronomy using neural networks[END_REF][START_REF] Paine | Machine learning for avoiding stagnation in image-based wavefront sensing[END_REF] is particularly challenging. First of all, sampling the instrumental aberrations across the fov is feasible by the use of internal fibers within the system, [START_REF] Ciurlo | PSF reconstruction for integral field spectrographs[END_REF][START_REF] Sitarski | Modeling instrumental field-dependent aberrations in the nirc2 instrument on the keck ii telescope[END_REF] but can be potentially time consuming and requires direct control of the system. In other words, this must be anticipated as part of a PSF reconstruction project and can not be easily scheduled at short times scales. Another issue with static aberrations occurs with on-sky calibration as the PSF combines the instrumental part plus the residual atmospheric contribution that blurs the instrumental signature in the data. Current on-sky tests usually relies on the acquisition of long-exposure pairs of in and out-of-focus images, which can be particularly telescope time consuming, [START_REF] Ragland | Point spread function determination for Keck adaptive optics[END_REF][START_REF] Jolissaint | Adaptive Optics Point Spread Function Reconstruction at W. M. Keck Observatory in Laser Natural Guide Star Modes : Final Developments[END_REF] especially if one wants to map the instrumental aberrations with respect to the telescope elevation. Regardless the pupil-plane-based PSF determination approach, the calibration of static aberrations is mandatory to achieve very accurate (better than a few percents of accuracy on PSF figures of merits) PSF determination. Moving to pupil-plane based techniques is a major change of paradigm: it is far beyond from the sole processing of focal-plane images and involves multi-disciplinary expertise and hands-on work on the system to achieve successful PSF determination.

The second term in Eq. 1 is systematically calculated from 29

hAO (ρ/λ) = exp 2π λ 2 × (B AO (ρ) -B AO (0)) , (2)
where B AO is the covariance function of the residual optical path difference downstream AO correction. There are multiple branches that go from this point to the PSF and we draw here an non-exhaustive list.

• Fourier model-based approach. This consists of modeling how spatial frequencies of atmospheric phase disturbances propagate from turbulent layers to the detector focal plane. This requires to describe analytically each component of the AO system so as to understand how a spatial frequency is measured and reconstructed through the wavefront-sensing process. 15-17, 38, 39 The covariance function is then obtained through the Wiener-Khintchine theorem that assures that if the residual phase is stationary, we have:

B AO (ρ) = F -1 [W AO (k)] (3)
where W AO is the AO residual phase Power Spectrum Density (PSD). This latter is split into multiple terms that account for different physical limitations, such as the DM cut-off frequency that creates the seeinglimited PSF halo, or the servo-lag error that models how the system latency combined with the atmosphere temporal characteristics create a PSF elongation. The Fourier calculation assumes intrinsically that the residual phase is spatially stationary, meaning its variance remains identical regardless of the position within the pupil plane. which does not apply for modeling the cone effect, although some techniques can be deployed to mitigate this issue. [START_REF] Neichel | Tomographic reconstruction for wide-field adaptive optics systems: Fourier domain analysis and fundamental limitations[END_REF] The Fourier-based approach is very fast to compute (less than 0.5 s on a standard CPU with 200 pixels to describe a 8 m pupil) and requires only a few parameters (seeing, C 2 n profile, DM conjugation altitude, LGS asterisms,...), which is very useful for predicting AO performance with respect to the system design and site observing conditions.

This technique is also a good representation of the post-AO PSF in poor correction regime, as it has been applied on MUSE WFM. [START_REF] Fusco | Reconstruction of the ground-layer adaptive-optics point spread function for MUSE wide field mode observations[END_REF] For such systems that address the AO-assisted observations of very large fov (¿1') and rely on a GLAO system, the Fourier-based model has been shown to produce an estimation of PSF figures of merits (FWHM, ensquared energy) at the level of 10-20% as presented in Fig. 2. At this level of correction (factor 2 improvement of the seeing), the calibration of static aberrations and complex models of the AO system were not necessary to meet the science requirements, and consequently, regarding the this method is an analytical framework, the Fourier approach offers the feasibility to design an operational and robust pipeline from the minimal amount of information from the telemetry.

• Covariance model-based approach. This is a similar approach that involves the same hypothesis on atmosphere statistics. Similarly to the Fourier method, the AO residual is split into error terms (actually the same physical processes) but focuses on modeling the auto-and cross-correlation of those. [START_REF] Ferreira | Numerical estimation of wavefront error breakdown in adaptive optics[END_REF][START_REF] Martin | PSF reconstruction validated using on-sky CANARY data in MOAO mode[END_REF][START_REF] Gendron | A novel fast and accurate pseudo-analytical simulation approach for MOAO[END_REF][START_REF] Villecroze | PSF reconstruction for MUSE in wide field mode[END_REF] In the case of a point-wise description of the residual phase over n × n samples, this covariance is calculated over n 2 × n 2 pixels, which means that we derive the covariance of each phase sample with all others. If the phase is spatially stationary, this covariance matrix would have a Toeplitz structure with only n × n non-redundant elements, and we would not lose information by averaging redundant baselines to get a covariance map, which is connected to the PSD through the Wiener-Khintchine theorem. Compared to the PSD approach, the covariance model-based calculation allows a simpler calculation of the cross-terms, separates properly the tip-tilt modes from the other modes, and accounts correctly for the LGS cone effect, DM spatial filtering and phase spatial non-stationarity in general, which reveal to be particularly relevant in the case of PSF reconstruction for tomographic systems. • Reconstruction approach. The techniques presented above degrade the information contained in the telemetry down to few integrated parameters only. With the reconstruction approach, originally introduced by Véran, [START_REF] Veran | Estimation of the adaptive optics long-exposure point-spread function using control loop data[END_REF] B AO in Eq. 3 is derived from the empirical covariance of wavefront measurements, which contains more information (second order statistical moments) than the sole integrated parameters. This strategy generally allows improved PSF reconstruction accuracy despite the similar mathematical formalism. On-sky tests have been achieved through the two past decades, especially on PUEO, [START_REF] Veran | Estimation of the adaptive optics long-exposure point-spread function using control loop data[END_REF] NACO, [START_REF] Clénet | Tests of the PSF reconstruction algorithm for NACO/VLT[END_REF] Keck 24, 26 or CANARY, [START_REF] Martin | PSF reconstruction validated using on-sky CANARY data in MOAO mode[END_REF] that showed that reconstructing the PSF at the few percent level is doable. W.M. Keck Observatory has particularly pushed for enabling an accurate PSF reconstruction for the Keck SCAO-assisted observation, with estimation errors pushed down to few percents as presented in Fig. 3.

Major limitations seem to be the calibration of static or quasi-static aberrations (NCPA, segment differential piston,...) as well as optical gains variations (quad cell) that must be deeply characterized to break through this limit. There is currently no such PSF reconstruction verification on-sky for LTAO and MCAO systems. This is hypothetically due to lack of access to on-sky AO telemetry data from such operating systems (MUSE NFM, GeMS) and the need of complex recipes to model the PSF across the field. [START_REF] Gilles | Simulation model based approach for long exposure atmospheric point spread function reconstruction for laser guide star multiconjugate adaptive optics[END_REF] Note that this strategy is also envisioned for solar AO systems * • PSFAO19 model A recent analytical description has been proposed [START_REF] Beltramo-Martin | Joint estimation of atmospheric and instrumental defects using a parsimonious point spread function model[END_REF][START_REF] Fétick | Physics-based model of the adaptive-optics-corrected point spread function. Applications to the SPHERE/ZIMPOL and MUSE instruments[END_REF] to represent the PSD/PSF from physical-based parameters. This model jointly describes both the corrected and uncorrected parts of the PSF, with successful recent applications. [START_REF] Fétick | Closing the gap between Earth-based and interplanetary mission observations: Vesta seen by VLT/SPHERE[END_REF] In comparison of the previous discussed techniques, this model simplifies the PSF determination problem thanks to his parsimony and can be used as a parametric model (focal-plane-based approach) or be calibrated from the AO control loop data (pupil-plane-based approach). The parsimony is a major advantage here: with only 7 parameters (no static aberrations accounted for), the model shown excellent match over on-sky-observations as presented in Figure 4 that presents a comparison of on-sky PSFs obtained on 7 different AO-assisted instruments with corresponding fitted models. The model adapts to any kind of AO correction; the differences are at the level of a few percents among all systems regardless the imaging wavelength (Vis./near infra-red) and the pixel scale (undersampling/oversampling). For over 4,800 PSFs obtained on these instruments, the model succeeds in reproducing the PSF metrics (Strehlratio, FWHM) at less than 4% or error. [START_REF] Beltramo-Martin | Joint estimation of atmospheric and instrumental defects using a parsimonious point spread function model[END_REF] Note that this focal-plane-based/pupil-plane-based dual aspect could be also true for the Fourier-based model; we could distinguish a few parameters in this model to be retrieved on the image. Nevertheless, building such a model already requires AO expertise to describe the propagation of spatial frequencies, while the PSFAO19 model is a concise and parsimonious representation of the known AO PSD behavior, and still guarantee an excellent level of accuracy. The only drawback of the current implementation is that secondary parameters are not yet clearly connected to the AO control loop data, and seems to vary across the field and wavelength as well. This is an on-going analysis that is currently done in the framework of the ANR APPLY.

Hybrid reconstruction

Hybrid PSF determination techniques gather both focal-plane, pupil-plane measurements and analytical models of atmospheric/instrumental defects for achieving the best reconstruction possible. The presence of a star in the field does not guarantee our capability to identify accurately the PSF shape at any position in the field, especially regarding field PSF variations, source crowding and noise contamination. Consequently, constraining the PSF shape from a model that inherits from the AO telemetry allows to calibrate the PSF reconstruction process. For instance, accounting for the spatial decorrelation of phase disturbances so-called anisoplanatism effect, [START_REF] Fried | Anisoplanatism in adaptive optics[END_REF] and calibrating field-dependent static aberrations, 37 one may constrain the PSF shape at a point B in the field, knowing the PSF at a point A. Starting from this idea, a new Starfinder-based algorithm, so-called AIROPA, [START_REF] Ciurlo | PSF reconstruction for integral field spectrographs[END_REF][START_REF] Turri | Testing PSF reconstruction on sky with AIROPA[END_REF][START_REF] Witzel | The AIROPA software package: milestones for testing general relativity in the strong gravity regime with AO[END_REF] relies on a PSF extraction and characterization of PSF spatial variations through anisoplanatism models and calibration of field-dependent static aberrations. In a sense, it fits the hybrid reconstruction class regarding the level of expertise needed in AO models and calibration.

In the same spirit, the algorithm called PRIME [START_REF] Beltramo-Martin | PRIME: PSF Reconstruction and Identification for Multiple-source characterization Enhancement -application to Keck NIRC2 imager[END_REF] combines PSF fitting and PSF reconstruction techniques to calibrate the PSF model. Contrary to PSF reconstruction, atmospheric and system parameters on which the model relies are no longer fixed; they can be retrieved on the focal-plane image directly, which guarantees an optimal use of the AO telemetry for fitting the PSF For instance, the seeing value is determined from the PSF wings. On data acquired with SPHERE/ZIMPOL in V and NIRC2 at Keck in H/K, a 1% of error on the 2D PSF estimation has been obtained. [START_REF] Beltramo-Martin | PRIME: PSF Reconstruction and Identification for Multiple-source characterization Enhancement -application to Keck NIRC2 imager[END_REF][START_REF] Beltramo-Martin | Pushing Point spread function reconstruction to the next level. Application to SPHERE/ZIMPOL[END_REF] Hybrid techniques are particularly well suited to science cases that offer access to stars in the field.

Marginal PSF estimation (& deconvolution) from images of extended objects

Deconvolution refers to technique that aim to restoring high-spatial frequencies that are filtered out through the optical system. [START_REF] Fétick | Blind deconvolution in astronomy with adaptive optics: the parametric marginal approach[END_REF][START_REF] Blanco | Marginal blind deconvolution of adaptive optics retinal images[END_REF][START_REF] Mugnier | MISTRAL: a myopic edge-preserving image restoration method, with application to astronomical adaptive-optics-corrected long-exposure images[END_REF][START_REF] Mugnier | Myopic deconvolution from wave-front sensing[END_REF] In the isoplanetic domain, the imaging equation is

i = o h + n ( 4
)
where i is the image, o is the object, h is the PSF and n is an additive noise. Deconvolution is an estimation ô of the object that produced the image i and performed by the minimization of the following criterion

J (o, h) = ||W(i -o h)|| 2 + µ o R o (o) + µ h R h (h) (5)
where the first term is the data-attached term that is the quadratic distance from the model to the data and weighted by the noise variance, [START_REF] Fétick | Blind deconvolution in astronomy with adaptive optics: the parametric marginal approach[END_REF][START_REF] Mugnier | MISTRAL: a myopic edge-preserving image restoration method, with application to astronomical adaptive-optics-corrected long-exposure images[END_REF] and where R o (o) and R h (h) are penalty terms on respectively the object and the PSF. In practice, estimating jointly the object and the PSF (myopic deconvolution) is a highly undertermined problem and not feasible and the PSF can converge towards the Dirac function. [START_REF] Thiébaut | Strict a priori constraints for maximum-likelihood blind deconvolution[END_REF] Using a parametric model to describing the PSF and adjust jointly the object and PSF parameters also leads to algorithm degeneracy. [START_REF] Fétick | Blind deconvolution in astronomy with adaptive optics: the parametric marginal approach[END_REF][START_REF] Blanco | Marginal blind deconvolution of adaptive optics retinal images[END_REF] The sucessful approaches so far consisted in providing a strong prior on the PSF, such as

R h (h) = 0.5 f | h - hm | 2 /S h (f )
, with hm the mean OTF and S h (f ) the associated spatial PSD, which must be estimated empirically.

To circumvent this need of particular PSF observations, a technique is currently being developed to estimate the PSF directly from the image of a resolved object. This technique is called "marginal PSF estimation" and has successfully been applied in the past for biomedical images. [START_REF] Blanco | Marginal blind deconvolution of adaptive optics retinal images[END_REF] The approach has been adapted to astronomical images, [START_REF] Fétick | Blind deconvolution in astronomy with adaptive optics: the parametric marginal approach[END_REF] with direct applications to asteroid imaging for example. The basic principle of the method is to state that the observed image contains as much information on the object as on the PSF itself. Using a statistical approach, the object is integrated -marginalised -out of the problem to focus on the PSF estimation. A parametric model of the PSF is used so it reduces the number of unknowns to estimate by the marginal technique. Once the PSF parameters are known, we directly retrieve the PSF shape. Then this PSF is used to perform deconvolution of the image to restore the details of the observed object. In general, deconvolution processes are ill-posed problems highly sensitive to the PSF shape. The marginal PSF estimation provides an accurate PSF shape, usable for this difficult problem of deconvolution. The marginally estimated PSF may also be applied in any post-processing technique requiring the knowledge of the PSF.

Discussion on uncertainties

There is no real consensus on the estimation of uncertainties on the PSF through the PSF reconstruction process. Some work have been conducted in the past to address this question. [START_REF] Jolissaint | Exploring the impact of PSF reconstruction errors on the reduction of astronomical adaptive optics based data[END_REF][START_REF] Exposito | Estimation of errors on the PSF reconstruction process for myopic deconvolution[END_REF] In the case where the PSF is determined from an analytical model of the residual phase PSD (Fourier, covariance model-based), the PSF uncertainties are directly driven by the uncertainties on the model inputs parameters (seeing, C 2 n (h), L 0 ,...) and one may easily establish this connection. Note that biases are (must be) compensated by comparing the PSF estimator with on-sky observations. [START_REF] Fusco | Reconstruction of the ground-layer adaptive-optics point spread function for MUSE wide field mode observations[END_REF] Fig. 5 illustrates the photometric bias and precision obtained on simulated HARMONI PSF in SCAO mode, including the real ELT pupil and the cophasing error map. The PSF were obtain from the FAST simulator [START_REF] Neichel | A NEW TOOL TO EFFICIENTLY PREDICT YOUR FAVORITE AO PSF[END_REF] We have defined a distribution of seeing S = N T (µ s , σ s ), where N T is a truncated normal distribution with µ s = 0.65 the mean-value of the distribution and the reference seeing value, σ s the precision we set, and s > 0 ∀s ∈ S. We measured the distribution on the photometric error through a two-folds process: (i) compute PSF σ from a selection s ∈ S with σ s defined from 1% to 50% of µ s , and (ii) perform a linear regression to obtain a = argmin

a ||PSF 0 -a × PSF σ | | 2
, where PSF 0 is the ground truth, and measure the photometric error as Figure 5: Left: H-band photometric bias and precision (1-σ) obtained on simulated HARMONI PSFs with respect to the precision on the seeing estimates (bias on seeing estimate is assumed to be compensated). Right: histograms of the seeing and photometric error distributions for σs = 0.1µs. ∆m = -2.5log 10 ( a). First of all, the distribution of the photometric error is N (µ m , σ m ), a Normal distribution as well, with a mean µ m and standard-deviation σ s whose estimates are function of σ s as illustrated in Fig. 5.

For the usual 10% precision we obtain on seeing measurements, [START_REF] Andrade | Estimation of atmospheric turbulence parameters from Shack-Hartmann wavefront sensor measurements[END_REF][START_REF] Jolissaint | Determination of the optical turbulence parameters from the adaptive optics telemetry: critical analysis and on-sky validation[END_REF] we can obtain a marginally bias photometry measurements with a precision of 0.06 mag at 1 -σ, accounting on the sole precision on the seeing estimation (no noise) and assuming that the PSF is perfectly described by the Fourier model. However, if WFS time series are used (PSF reconstruction), the seeing estimate impacts the PSF wings model only (we stick to the isoplanatic case), while in the present analysis, the seeing shapes the PSF cores as well. Considering that the Fourier model is a perfect representation of the PSF wings (actually the PSF wings are modeled using this approach regardless the pupil-plane determination technique), this methodology reveals the upper-bound on the photometric error due to seeing estimates precision.

If the reconstruction process is more complex, the estimation of the uncertainties is not completely straightforward, but still feasible if we identify what parameters are critical, such as the optical gains for instance or the seeing estimates. Other aspect, such as accuracy of the temporal transfer function model or even DM influence function shape, matter but might not be the main cause of the uncertainty, although this deserves a more complete analysis that will be conducted in the future.

APPLICATIONS TO POST-PROCESSING

Characterization of brown dwarfs binary systems with PSF-R using NIRC2@Keck

One of the two science cases identified for the first phase of science verification is low-mass binaries. The binaries used for our science verification tests are sourced from a sample that has been regularly monitored with Keck AO, both NGS and LGS, for essentially the entire history of these Keck facilities (e.g.,; 60 ; 6162 ). We used the most recent orbital analysis from [START_REF] Dupuy | Individual Dynamical Masses of Ultracool Dwarfs[END_REF] to select binaries with precisely known orbits that can provide both on-and off-axis tests with NGS or LGS AO, although we have focused on testing on-axis PSF-R to date. Gl 569Bab (M8.5+M9) and HD 130948BC (L4+L4) can be observed as on-axis NGS targets because they are both in hierarchical triple systems with bright primary stars only 4.9" and 2.6" away, respectively. LP 349-25AB (M7+M8) and 2MASS J0746+2000AB (L0+L1.5) are an on-axis LGS targets. Of all these systems, Gl 569Bab is perhaps the most useful because its primary star is only 4 mag brighter than the binary components, so unsaturated images of all three components can be obtained, with the primary providing an empirical PSF reference only 4.9" off-axis.

Without PSF-R, fitting binary images requires either assuming an analytical PSF model (e.g., multi-component, elliptical Gaussians) or deriving an empirical PSF from the image itself (e.g., StarFinder; 1 ). Generally, the former is necessary for the tightest binaries where the two components are barely distinguishable from a single elongated PSF, while for wider binaries with more cleanly separated component PSFs the StarFinder algorithm can successfully identify two stars and estimate the PSF that they share in common. The median astrometric error of tight (separations of 40 -80 mas) binaries from 63 that were predominantly fit using the analytic PSF model is 0.8 mas. In contrast, wider binaries with 100-200 mas separations analyzed with StarFinder yield typical errors of 0.1 mas (1% of a NIRC2 pixel). In terms of photometry, StarFinder yields flux ratios with a typical precision of 0.01 mag as judged from the repeatability of the measurements across multiple epochs.

In principle, PSF-R should not significantly improve on StarFinder when the latter performs properly, regarding StarFinder capabilities to optimally determine the PSF without any theoretical assumptions about the atmosphere and telescope optics. Therefore, a reasonable goal of PSF-R would be to achieve 0.1 mas errors for the tightest binaries, i.e., accuracy comparable to what is possible when the PSF can be determined directly from the imaging data. The reconstructed PSFs are used in the science data analysis, in conjunction with the standard methods, to assess the performance. Essentially, binary-fitting analysis is performed for different cases such as using the reconstructed PSFs (PSF-R), the on-sky primary star PSF (GL569A), and the 2D Gaussian model. The resultant binary model parameters are compared with the reference values estimated from binary orbit fit in Tab. 1: Binary parameters estimated from different methods using 40 exposures of NIRC2 in H-band. These measurements, taken at just one binary separation, position angle and delta magnitude are not necessarily statistically representative of performance in general, but they do provide chance to look for possible issues. For example, the agreement between a simple Gaussian fit and the orbit astrometry is surprisingly good, but the disagreement in the magnitude difference shows that there is more to the story (e.g., perhaps the systematic errors in the relative astrometry cancel out in this fit better than those in the relative photometry). Such agreement between a simple 2D Gaussian model and true parameters are rarely seen in simulations of such binaries and thus require further investigations. The PSF-R results agree quite well with both the orbit baseline truth, and the fit using Gl 569A as the template. In particular, the magnitude difference is better reproduced by the PSF-R rather than the Gl 569A, indicating that the on-axis PSF-R may indeed be superior to the off-axis observed PSF as it includes the spatial variations of the PSF due to anisoplanatism and field-dependent static aberrations. [START_REF] Sitarski | Modeling instrumental field-dependent aberrations in the nirc2 instrument on the keck ii telescope[END_REF] Two other science cases, namely, the morphology of QSO host galaxies and microlensing, in addition to low-mass binaries, are considered for the next phase of science verification.

Improving astrometry/photometry measurements on SPHERE/ZIMPOL

As part of a technical proposal following the ESO Calibration Workshop 2017 † , Massari et. al observed the core of the Galactic globular cluster NGC 6121 with the AO assisted camera SPHERE/ZIMPOL 64 in V-band, and under the program ID 60.A-9801(S). More details of the observation are given in [START_REF] Massari | Precise photometry and astrometry of the globular cluster NGC 6121 using PSF-R techniques[END_REF][START_REF] Massari | Successful application of PSF-R techniques to the case of the globular cluster NGC 6121 (M 4)[END_REF] The observed field illustrated in Fig. 6 contains only five stars in a 3.5"×3.5" field of view, where the brightest star has a Vband magnitude of V'=10.9 mag, and has been used as the natural guide star for the AO system. Photometric measurements of the central star were obtained using four different techniques over 24 successive acquisitions, (i) aperture photometry, (ii) DAOPHOT-II, [START_REF] Stetson | DAOPHOT -A computer program for crowded-field stellar photometry[END_REF] (iii) PSF-fitting using the PSFAO19 model, [START_REF] Fétick | Physics-based model of the adaptive-optics-corrected point spread function. Applications to the SPHERE/ZIMPOL and MUSE instruments[END_REF] and (iv) PSF-fitting using PRIME. The analysis illustrates that numerical PSF models, such as PRIME or PSFAO19, reach more precise measurements than the a-posteriori techniques, both photometrically (by a factor of 10) and astrometrically (by a factor of 4), advocating for the potential of innovative PSF-modeling methods as a powerful tool for analyzing future AO-assisted observations.

Towards an upgrade of Starfinder illustrated on GeMS/GSAOI

Another illustration of the potential of advanced PSF models for photometry/astrometry measurements using PSF-fitting is illustrated in Fig. 7, that shows astrometric error histograms and color-magnitude diagrams obtained by deploying Starfinder 1, 67, 68 on R136 simulated images based on a GeMS/GSAOI configuration with a pixel scale of 20 mas/pixel in H and K bands. The PSF were simulated thanks to the Fourier PSF simulator [START_REF] Neichel | Tomographic reconstruction for wide-field adaptive optics systems: Fourier domain analysis and fundamental limitations[END_REF] which is broadly used in the community, in the HSIM simulator for HARMONI [START_REF] Zieleniewski | HSIM: a simulation pipeline for the HARMONI integral field spectrograph on the European ELT[END_REF] for instance. Fig. 7 illustrates that the use of the numerical model PSFAO19, [START_REF] Fétick | Physics-based model of the adaptive-optics-corrected point spread function. Applications to the SPHERE/ZIMPOL and MUSE instruments[END_REF] that is preliminary adjusted on simulated PSF, allows to unbias astrometry measurements and improve precision by a factor 4, while photometric measurement accuracies are enhanced by up to a factor 10. A current work is focusing on an upgrade of Starfinder in order to include the analytical expression of the PSFAO19 model to improve the whole processing chain performed in Starfinder.

Exploitation of MUSE NFM data with PampelMuse

Advanced PSF models can also be employed in the analysis of integral field spectroscopy, for example obtained using MUSE [START_REF] Bacon | The MUSE second-generation VLT instrument[END_REF] in the narrow field mode (NFM). MUSE NFM data of several globular clusters were obtained as part of the guaranteed time observations, in the star clusters program. [START_REF] Kamann | PampelMuse: Crowded-field 3D spectroscopy[END_REF] In order to extract individual stellar spectra from the data, the reduced cubes are processed with the PampelMuse software. [START_REF] Kamann | Resolving stellar populations with crowded field 3D spectroscopy[END_REF] PampelMuse first tries to determine the PSF of the observation as well as the coordinates of the resolved sources as a function of wavelength. In a second step, this information is used to perform an optimal extraction of the spectrum of each resolved star. As illustrated in Fig. 2, the default PSF model used by PampelMuse is a Moffat profile. While this profile provides a good representation of the actual PSF of MUSE data in the wide fied mode (WFM), it cannot account for the complexity of the NFM PSF. Therefore, the PSFAO19 22 model was implemented in PampelMuse. In Fig. 8, we illustrate the improvement in performance this model provides relative to a Moffat PSF in a MUSE NFM observation of the globular cluster NGC 5139. The central panel of Fig. 8 shows that even when the Moffat profile is complemented by a Gaussian core in order to represent the diffraction-limited part of the PSF, strong fit residuals around the bright stars remain. Instead, using the dedicated PSFAO19 model strongly suppresses these residuals.

Deconvolution of asteroids images obtained with SPHERE/ZIMPOL

The biggest asteroids of the Main Belt have been observed ‡ at high resolution with the VLT SPHERE/ZIMPOL imager. The SPHERE instrument benefits from the SAXO extreme adaptive optics system. Figure 9 (middle image) shows the observation 46 of the asteroid Vesta with SPHERE/ZIMPOL on June 8th, 2018. This image represents the state-of-the-art ground based instruments with current extreme AO systems, on which the overall shape is captured but finest details are attenuated. The result of the deconvolution for Vesta can be seen on Fig. 9 (right). It exhibits sharp edges, and finer details such as local albedo variations or craters down to 20 -30 km when observed from Earth at more than 170 million kilometers away. More particularly, PSF reconstruction refers to techniques that estimate or constrain the post-AO PSF from AO control loop data, e.g. time series of WFS/DM measurements and associated calibration. Multiple implementations (Fourier, covariance-based, reconstruction) may be envisioned, depending on the needs of PSF accuracy, RTC flexibility, data flow and data storage capabilities. On the one hand, there is the historical PSF reconstruction approach proposed by Véran et. al that derives the PSF, with a potential accuracy of 1% on PSF metrics or PSF mean square error, from the time-history of WFS measurements, which may be time and computational-power-demanding and constraining for the RTC data flow control. On the other hand, there are analytical and parsimonious techniques that require access to integrated parameters only (seeing, C 2 n (h) profile, photons return,...) and allow to model the PSF with a few up to 10% of accuracy on PSF metrics.

From the past work that invoked PSF reconstruction for the processing of AO-corrected images, we distinguish several scenarios:

• Very accurate PSF determination for analysis of crowded regions. This seems to be the most demanding science case in terms of needed accuracy on the PSF. This is also a case for which stars are present in the field and the most efficient way to benefit both PSF reconstruction and phased-shifted (wavelength-shifted) representation of the PSF in the focal plane is to go for hybrid PSF reconstruction, which is expected to beat the 1% of accuracy on the PSF recovery in the field (and in wavelength as well).

• Moderate ( 10%) up to very accurate PSF determination for the processing of images that contains stars. If point sources are observed in the field, an alternative way to calibrate the PSF model in the field is to rely on a parametric model of the PSF that also offers the possibility to take advantages of the spatial and spectral representation of the PSF. Such an approach is particularly meaningful when AO control loop data are not available, and also for tomographic AO systems for which an extension of the PSF reconstruction framework is still missing.

• Moderate PSF determination for the processing of stars-free astronomical images. In this situation, one must rely on AO control loop and use, either the classical PSF reconstruction recipe, or an analytical model. In both cases, the estimation process must be calibrated through a sample of on-sky observed PSF.

It is also crucial to keep in mind that PSF reconstruction can only reconstruct spatial frequencies to which the WFS is sensitive enough to deliver a linearly-dependent wavefront representation; non-linear wavefront sensors, high-order spatial frequencies and unseen frequencies (NCPA for instance) must be derived from alternative techniques. One of the advantage of parametric PSF models is that any fitting procedure is sensitive (noisedependent) to all spatial frequencies that feature the PSF shape and can be rendered thanks to an appropriate calibration. For instance, accurate photometry requires to model the AO PSF wings formed by uncompensated spatial frequencies; the use of the time-history WFS measurements is not needed in this situation as all pupilbased PSF determination techniques will rely on the same model; the PSF wings are derived from the Kolmogorov PSD model of the atmospheric phase that depends on the seeing estimation only, that maybe obtained from the AO control data or from the focal-plane image directly.

Finally, any decision to deploy such or another PSF reconstruction approach must absolutely be driven by science and user interface requirements before getting into technical implementation.

Figure 1 :

 1 Figure 1: Summary of existing PSF determination techniques and associated references.

Figure 2 :

 2 Figure2: Estimated FWHM from the Fourier reconstruction versus measurements obtained over 94 MUSE WFM PSFs. Three spectrum channels are considered here and color indicates the seeing conditions.[START_REF] Fusco | Reconstruction of the ground-layer adaptive-optics point spread function for MUSE wide field mode observations[END_REF]

Figure 3 :

 3 Figure 3: The reconstructed NGS PSF (middle) is compared with the target PSF (left) for a high Strehl ratio case (50%). The images were taken at 1.655 µm. Right: a scan across the PSF through the peak is shown in logarithmic scale 25

Figure 4 :

 4 Figure 4: Two dimensional comparison in log scale and over 64×64 pixels of (top) observed PSFs, (middle) fitted model for multiple types of AO systems and instruments working in either visible or NIR wavelengths. From left to right: SPHERE/ZIMPOL, SPHERE/IRDIS, Keck AO/NIRC2, SOUL/LUCI, CANARY/CAMICAZ in MOAO mode, GALACSI/MUSE NFM, and GEMS/GSAOI. All PSFs have been normalized to the max.

28

 28

Figure 6 :

 6 Figure 6: Left: Average of the SPHERE/ZIMPOL observations of NGC 6121. The five detected stars are marked with red circles. Right: Photometric (top panel) and astrometric (bottom panel) precision achieved on the guide star when using the aperture method (first column), DAOPHOT-II (second column), the PSFAO19 model (third column), and PRIME (last column).

Figure 7 :

 7 Figure 7: Left: astrometric error obtained with Starfinder applied to simulated images of R136 observed with GSAOI at Gemini, using the native PSF identification routine or the PSFAO19 model preliminary adjusted on simulated PSFs. Right: Corresponding color-magnitude diagrams. The empirical PSF is a median value of several PSF extracted from the image. In order to account for spatial variations, the image is segmented over 7x7 regions.

Figure 8 :

 8 Figure 8: Left : Image of NGC5139 obtained with MUSE NFM at 850 nm. Middle: residual map after stars detection and PSF subtraction using a Moffat-Gauss model for describing the PSF. Right: Residual map obtained using the PSFAO19 model. The mean residual error is improved by a factor 3 thanks to the new model definition.

Figure 9 :

 9 Figure 9: Left: Synthetic high-resolution image of VESTA produced by OASIS (L.Jorda, LAM) from the NASA/Dawn data. Middle : SPHERE/ZIMPOL observation in R band. Right : deconvolution using the marginal approach.50

50

 50 Figure 9: Left: Synthetic high-resolution image of VESTA produced by OASIS (L.Jorda, LAM) from the NASA/Dawn data. Middle : SPHERE/ZIMPOL observation in R band. Right : deconvolution using the marginal approach.50

  4. SUMMARY AND CONCLUSIONSWe have classified the PSF determination techniques into three classes: (i) focal-plane-based techniques, (ii) pupil-plane-based techniques, and (iii) hybrid techniques that combine observed image and models of models of atmospheric/instrumental aberrations.

  1.

		Binary separation (mas) Position angle (degrees) Delta magnitude (mag)
	Orbit	67.18 ± 0.24	169.2 ± 0.4	0.47 ± 0.04
	PSF-R	66.7 ± 0.8	169.4 ± 0.8	0.46 ± 0.03
	Gl 569 A	66.9± 0.6	169.9 ± 0.5	0.53 ± 0.03
	Gaussian	67.1 ± 1.0	169.6 ± 0.6	0.64 ± 0.05
	Table

* See L. Jolissaint presentation at AO4ASTRO: https://nuage.osupytheas.fr/s/9eS7WSXJ3ww57t5

† www.eso.org/sci/meetings/2017/calibration2017.html

‡ ESO Large Progam ID 199.C-0074 (PI: P. Vernazza)

ACKNOWLEDGMENTS

This work has been partially funded by the French National Research Agency (ANR) program APPLY -ANR-19-CE31-0011. This work also benefited from the support of the WOLF project ANR-18-CE31-0018 of the and the OPTICON H2020 (2017-2020) Work Package 1.