Qiyu Jin

Gabriele Facciolo
email: facciolo@cmla.ens-cachan.fr

Jean-Michel Morel
email: moreljeanmichel@gmail.com

A Review of an Old Dilemma: Demosaicking First, or Denoising First?

Image denoising and demosaicking are the most important early stages in digital camera pipelines. They constitute a severely ill-posed problem that aims at reconstructing a full color image from a noisy color filter array (CFA) image. In most of the literature, denoising and demosaicking are treated as two independent problems, without considering their interaction, or asking which should be applied first. Several recent works have started addressing them jointly in works that involve heavy weight CNNs, thus incompatible with low power portable imaging devices. Hence, the question of how to combine denoising and demosaicking to reconstruct full color images remains very relevant: Is denoising to be applied first, or should that be demosaicking first? In this paper, we review the main variants of these strategies and carry-out an extensive evaluation to find the best way to reconstruct full color images from a noisy mosaic. We conclude that demosaicking should applied first, followed by denoising. Yet we prove that this requires an adaptation of classic denoising algorithms to demosaicked noise, which we justify and specify.

Introduction

Most digital cameras capture image data by using a single sensor coupled with a color filter array (CFA). At each pixel in the array, only one color component is recorded and the resulting image is called a mosaic. The most common CFA is the Bayer color array [START_REF] Bayer | Color imaging array[END_REF], in which two out of four pixels measure the green component, one measures the red and one the blue. The process of completing the missing red, green and blue values at each pixel is called demosaicking. Noise is inevitable, especially in low light conditions and for small camera sensors like those used in mobile phones. The conventional approach in image restoration pipelines for processing noisy raw sensor data has long been to apply denoising and demosaicking as two independent steps [START_REF] Paliy | Demosaicing of noisy data: spatially adaptive approach[END_REF]. Furthermore, the immense majority of image processing papers addressing one of both operations do not address its combination with the other one. All classic denoising algorithms have been designed for color or grey level images with white noise added. Yet the realistic data are different: either a mosaic with white noise or a demosaicked image with structured noise.

Joint denoising/demosaicking methods. This has led several recent works to propose joint demosaicking-de-noising methods [START_REF] Chatterjee | Noise suppression in low-light images through joint denoising and demosaicing[END_REF][START_REF] Gharbi | Deep joint demosaicking and denoising[END_REF][START_REF] Hirakawa | Joint demosaicing and denoising[END_REF][START_REF] Khashabi | Joint demosaicing and denoising via learned nonparametric random fields[END_REF]. For example [START_REF] Goossens | An overview of state-ofthe-art denoising and demosaicking techniques: toward a unified framework for handling artifacts during image reconstruction[END_REF] proposed a variational model to solve jointly demosaicking, denoising and deblurring. It uses a sparsifying prior based on wavelet packets and applied on decorrelated color channels. More detail about the technicalities of this sophisticated method can be found in [START_REF] Aelterman | Computationally efficient locally adaptive demosaicing of color filter array images using the dual-tree complex wavelet packet transform[END_REF]. Life has become far easier for joint denoising/demosaicking with the emergence of machine learning methods. It is, indeed, easy to simulate as many learning data as needed. This methodology can be used to obtain ground breaking demosaicking algorithms such as [START_REF] Syu | Learning deep convolutional networks for demosaicing[END_REF]. This paper proposed in 2018 a demosaicking CNN outperforming the best handcrafted algorithms including ARI [START_REF] Monno | Adaptive residual interpolation for color and multispectral image demosaicking[END_REF] by nearly 2 decibels. In [START_REF] Khashabi | Joint demosaicing and denoising via learned nonparametric random fields[END_REF] a public ground truth dataset was introduced and used for one of the first joint demosaicking and denoising methods based on machine learning. In a rapid succession, two state of the art denoising+demosaicking methods based on deep learning were proposed: [START_REF] Gharbi | Deep joint demosaicking and denoising[END_REF] and [START_REF] Klatzer | Learning joint demosaicing and denoising based on sequential energy minimization[END_REF]. This last paper performs joint denoising and demosaicking by a customized neural network presented as a cascade of energy minimization methods tuned by learning. The outstanding results beat the claimed anterior best method [START_REF] Heide | Flexisp: A flexible camera image processing framework[END_REF] by 1 decibel. Then in 2018 we have two still better performing methods, [START_REF] Dong | Joint demosaicing and denoising with perceptual optimization on a generative adversarial network[END_REF] involving a GAN, which compares favorably to [START_REF] Gharbi | Deep joint demosaicking and denoising[END_REF] and [START_REF] Tan | Color image demosaicking via deep residual learning[END_REF]. The method recently proposed in [START_REF] Kokkinos | Iterative joint image demosaicking and denoising using a residual denoising network[END_REF] performs joint denoising and demosaicking by inserting many residual denoising layers in a CNN. This complex method is claimed to beat [START_REF] Gharbi | Deep joint demosaicking and denoising[END_REF] and [START_REF] Klatzer | Learning joint demosaicing and denoising based on sequential energy minimization[END_REF] by a good margin. Lastly, in 2019, [START_REF] Ehret | Joint demosaicing and denoising by overfitting of bursts of raw images[END_REF] introduced a "mosaic-to mosaic" training strategy analog to the noiseto-noise [START_REF] Lehtinen | Noise2Noise: Learning Image Restoration without Clean Data[END_REF] and frame-to-frame [START_REF] Ehret | Model-blind Video Denoising Via Frame-to-frame Training[END_REF] frameworks to handle noisy mosaicked raw data, and training both demosaicking and joint denoising and demosaicking networks without requiring ground truth. The method starts from pairs or bursts of raw images of the same scene. It registers them and learns to predict the missing colors.

Yet the question of how to combine denoising and demosaicking algorithms conceived as independent blocks remains very relevant, especially in the context of low power or portable devices, and given the fact that the main effort in denoising and demosaicking has addressed them independently. A big argument in favour of performing denoising before demosaicking is that most existing demosaicking algorithms have been developed under the unrealistic assumption of noise-free data [8, 10, 20, 23, 24, 28, 33-36, 38, 42, 45, 49, 53-56, 60]. Yet the perfor-mance of these algorithms can degrade dramatically when the noise level increases on the CFA raw image. Therefore, a previous denoising step is implicitly required by these algorithms.

In this paper we focus on the early CFA processing in the imaging pipeline (operating in linear space). We assume that the noise in the raw mosaic is additive white Gaussian (AWGN) and that its variance is known. This is realistic because, first, a variance stabilizing transform (VST) [START_REF] Anscombe | The Transformation of Poisson, Binomial and Negative-Binomial Data[END_REF] applied to a raw image results in a nearly AWG noise and, second, because an accurate noise model is often known or can be estimated [START_REF] Colom | Analysis and Extension of the Ponomarenko et al. Method, Estimating a Noise Curve from a Single Image[END_REF][START_REF] Ponomarenko | An automatic approach to lossy compression of AVIRIS images[END_REF]. In general, image denoising methods can be grouped into two major categories, the model based methods such as nonlocal means [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF][START_REF] Jin | Nonlocal means and optimal weights for noise removal[END_REF][START_REF] Jin | Convergence theorems for the non-local means filter[END_REF], nlBayes [START_REF] Lebrun | A nonlocal bayesian image denoising algorithm[END_REF], CBM3D [START_REF] Dabov | Image denoising with block-matching and 3d filtering[END_REF] and WNNM [START_REF] Gu | Weighted nuclear norm minimization with application to image denoising[END_REF], and deep learning methods such as [START_REF] Jain | Natural image denoising with convolutional networks[END_REF][START_REF] Zhang | Beyond a gaussian denoiser: Residual learning of deep cnn for image denoising[END_REF]. The ensuing CNNs are sometimes flexible in handling denoising problems with various noise levels.

Our goal here is to determine which strategy is more advantageous for coupling demosaicking and denoising : Is it applying denoising and then demosaicking (which we will denote DN &DM : DN and DM indicate denoising and demosaicking respectively), or is it better to apply first demosaicking and then denoising (DM &DN )?

DN&DM methods (i.e. denoising then demosaicking): advantages and drawbacks. Many state of the art works [START_REF] Kalevo | Noise reduction techniques for bayer-matrix images[END_REF][START_REF] Paliy | Demosaicing of noisy data: spatially adaptive approach[END_REF][START_REF] Sung | A case for denoising before demosaicking color filter array data[END_REF][START_REF] Zhang | Pca-based spatially adaptive denoising of cfa images for single-sensor digital cameras[END_REF] support the opinion that DN &DM outperforms DM &DN . Their first convincing argument is that after demosaicking noise becomes correlated, thus losing its independent identically distributed (i.i.d.) white Gaussian property. This increases the difficulty of applying efficient denoising and actually seems to discard all classic algorithms, that mostly rely on the AGWN assumption. A second obvious argument is that the best demosaicking algorithms have been designed with noisefree images.

For example, Park et al. [START_REF] Sung | A case for denoising before demosaicking color filter array data[END_REF] considered the classic Hamilton-Adams (HA) [START_REF] Jr | Adaptive color plan interpolation in single sensor color electronic camera[END_REF] and [START_REF] Dubois | Frequency-domain methods for demosaicking of bayer-sampled color images[END_REF] for demosaicking, combined with two denoising methods, BLS-GSM [START_REF] Portilla | Image denoising using scale mixtures of gaussians in the wavelet domain[END_REF] and CBM3D [START_REF] Dabov | Color image denoising via sparse 3d collaborative filtering with grouping constraint in luminance-chrominance space[END_REF]. This combination raises the question of adapting CBM3D to a CFA. To do so, the authors apply a sparsifying 4D color transform to the 4-channel image formed by rearranging the Bayer pixels, apply BM3D to each channel, then apply the inverse color transform.

In the very same vein, in the BM3D-CFA method [START_REF] Danielyan | Cross-color bm3d filtering of noisy raw data[END_REF] BM3D is applied directly on the CFA color array. To do so, "only blocks having the same CFA configuration are being compared to build the 3D blocks. This is the only modification of the original BM3D". A little thought leads to the conclusion that this amounts to denoise four different mosaics of the same image before aggregating the four values obtained for each pixel. The authors compare two denoising algorithms with two different setups: a) filtering CFA as a single image and b) splitting the CFA into four color components, filtering them separately, and recombining back the denoised CFA image. This paper showed a systematic improvement over [START_REF] Zhang | Pca-based spatially adaptive denoising of cfa images for single-sensor digital cameras[END_REF]. They use Zhang-Wu [START_REF] Zhang | Color demosaicking via directional linear minimum mean square-error estimation[END_REF] as demosaicking method for their comparison of result after demosaicking. In our comparisons the method of [START_REF] Danielyan | Cross-color bm3d filtering of noisy raw data[END_REF] will be mentioned every time we consider the DN &DM setup with BM3D. We will nevertheless replace the demosaicking of [START_REF] Zhang | Color demosaicking via directional linear minimum mean square-error estimation[END_REF] by RCNN [START_REF] Tan | Color image demosaicking via deep residual learning[END_REF] or RI [START_REF] Kiku | Residual interpolation for color image demosaicking[END_REF], which clearly outperform it.

Similarly in [START_REF] Chatterjee | Noise suppression in low-light images through joint denoising and demosaicing[END_REF] denoising is performed by an adaptation of NL-means to the Bayer pattern, where only patches with the same CFA configuration are being matched. This paper formulates the demosaicking as a super-resolution problem, assuming that the observed values are actually averages of four values in the high resolution image. It then guides this super-resolution problem by the NL-means weights. The method is compared with [START_REF] Menon | Joint demosaicking and denoisingwith space-varying filters[END_REF] and [START_REF] Zhang | Pca-based spatially adaptive denoising of cfa images for single-sensor digital cameras[END_REF]. The authors of [START_REF] Zhang | Pca-based spatially adaptive denoising of cfa images for single-sensor digital cameras[END_REF] also propose an DN &DM method, where the demosaicking method is [START_REF] Zhang | Color demosaicking via directional linear minimum mean square-error estimation[END_REF] and the the denoising method is an adaptation of nl-Bayes [START_REF] Lebrun | A nonlocal bayesian image denoising algorithm[END_REF] to a Bayer pattern. First, the method extracts blocks with similar configuration in the Bayer array and groups them by similarity, then it applies to them PCA and a Wiener denoising procedure which can be also interpreted as an LMMSE. In our experiments, this PCA method [START_REF] Zhang | Pca-based spatially adaptive denoising of cfa images for single-sensor digital cameras[END_REF] will be considered every time we evaluate the DN &DM scheme (but combined with a more recent demosaicking algorithm such as RCNN [START_REF] Tan | Color image demosaicking via deep residual learning[END_REF]). The more recent paper [START_REF] Zhang | Joint denoising and demosaicking of noisy cfa images based on inter-color correlation[END_REF] involves similar arguments. This paper uses [START_REF] Alleysson | Linear demosaicing inspired by the human visual system[END_REF], a linear filter to extract the luminance from the CFA. Then it remarks that this luminance is correlated, so it applies a variant of NL-means that attempts to decorrelate the noise. The same method is applied to each downsampled color channel and the high frequency of the grey level is transported back to the color channels. This method under-performs with respect to others considered here, so we shall not include it to our final comparison tables. Nevertheless, it remains of interest as a fast method compatible with low power cameras. The paper proves that it has a performance very close to a combination of [START_REF] Zhang | Pca-based spatially adaptive denoising of cfa images for single-sensor digital cameras[END_REF] and [START_REF] Hirakawa | Joint demosaicing and denoising[END_REF].

The paper [START_REF] Patil | Poisson noise removal for image demosaicing[END_REF] is another method promoting denoising before demosaicking, involving dictionary learning methods to remove the Poisson noise from the single channel images prior to demosaicing. Experimental results on simulated noisy images as well as real camera acquisitions, show the advantage of these methods over approaches that remove noise subsequent to demosaicing. The paper nevertheless uses [START_REF] Henrique S Malvar | Highquality linear interpolation for demosaicing of bayerpatterned color images[END_REF] which is a historic but outdated demosaicking method.

To summarize, in the DN &DM strategy all classic denoising algorithms such as CBM3D, nlBayes, nlMeans have been adapted to handle a noisy mosaic where only one of R, G or B is known at each pixel. Several of them [START_REF] Kalevo | Noise reduction techniques for bayer-matrix images[END_REF][START_REF] Paliy | Demosaicing of noisy data: spatially adaptive approach[END_REF][START_REF] Sung | A case for denoising before demosaicking color filter array data[END_REF][START_REF] Zhang | Pca-based spatially adaptive denoising of cfa images for single-sensor digital cameras[END_REF] address this realistic case by processing the noisy CFA images as a half-size 4-channel color image (with one red, two green and one blue channels) and then apply a multichannel denoising algorithm to it. The advantage of the denoising step of DN &DM is that the Poisson noise can be led back by the classic Anscombe transform to the case of i.i.d. white Gaussian, and the disadvantage is that the resolution of the image is reduced and, as a result, some details might be lost after denoising. Another issue of this strategy is that the spatial relative positions of the R, G, and B pixels are lost by handling the image as a four channel half size image.

In this paper, we address the above mentioned issues. We shall first delve into the advantages and disadvantages of DN &DM and DM &DN approaches. We shall then analyze noise properties after demosaicking and adjust two existing classic denoising algorithms (CBM3D and nlBayes) to accommodate them to this type of noise. Then, we shall perform a thorough experimental evaluation that will lead us to conclude that DM &DN (with an adjusted noise parameter) is superior to DN &DM . This result is opposite to the conclusion of [START_REF] Kalevo | Noise reduction techniques for bayer-matrix images[END_REF][START_REF] Paliy | Demosaicing of noisy data: spatially adaptive approach[END_REF][START_REF] Sung | A case for denoising before demosaicking color filter array data[END_REF][START_REF] Zhang | Pca-based spatially adaptive denoising of cfa images for single-sensor digital cameras[END_REF]. The advantages of DM &DN seem to be linked to the fact that this scheme does not handle a half size 4-channels color image; it therefore uses the classic denoising methods directly on a full resolution color image; this results in more details being preserved and avoids checkerboard effects.

The rest of the paper is structured as follows. In Sec-

The demosaicking and denoising framework

In a single-sensor camera equipped with a color filter array (CFA) [START_REF] Bayer | Color imaging array[END_REF], only one pixel value among the three RGB values is recorded at each pixel. Consider a CFA block as shown in Fig. 1. The raw Bayer CFA images are scalar mosaic matrices with noise. Obtaining high quality color images requires completing the missing color channels and removing the noise. As mentioned in the introduction, for this task we will consider two main schemes: DM &DN (demosaicking then denoising) and DN &DM (denoising then demosaicking). Park et al. [START_REF] Sung | A case for denoising before demosaicking color filter array data[END_REF] argued that demosaicking introduces chromatic and spatial correlations to the noisy CFA image. Then the noise is no longer i.i.d. white Gaussian, which makes it harder to remove. In [START_REF] Kalevo | Noise reduction techniques for bayer-matrix images[END_REF], some experiments were done to show that DN &DM schemes are more efficient to suppress noise than DM &DN schemes. Based on this argument several denoising methods [START_REF] Akiyama | Pseudo four-channel image denoising for noisy cfa raw data[END_REF][START_REF] Lee | Denoising algorithm for cfa image sensors considering inter-channel correlation[END_REF][START_REF] Sung | A case for denoising before demosaicking color filter array data[END_REF][START_REF] Zhang | Pca-based spatially adaptive denoising of cfa images for single-sensor digital cameras[END_REF] for raw CFA images before demosaicking were introduced. Other denoising methods that are not explicitly designed to handle raw CFA images (such as CBM3D and nlBayes) can also be applied on noisy CFA images by rearranging the CFA image into a halfsize four-channels image with two greens on which the denoising algorithm is applied [START_REF] Sung | A case for denoising before demosaicking color filter array data[END_REF]. The denoised CFA is then recovered by undoing the pixel rearrangement. How- In each experiment: left, detail of the demosaicked and denoised image; right, the difference with original that should contain mainly noise. DN : CBM3D denoising; DM : demosaicking (here we use RCNN). 1.5DN means that if noise level is σ 0 , the input noise level parameter of denoising method DN is σ = 1.5σ 0 ; DN &DM : uses the BM3D-CFA framework [START_REF] Danielyan | Cross-color bm3d filtering of noisy raw data[END_REF] for denoising.

ever, this strategy reduces the resolution of the image seen by the denoiser, and we observed checkerboard effects resulting from chromatic aberrations in the two green channels after denoising. To address this issue, Danielyan et al. [START_REF] Danielyan | Cross-color bm3d filtering of noisy raw data[END_REF] proposed BM3D-CFA which amounts to denoise four different mosaics of the same image before aggregating the four values obtained for each pixel.

Modeling demosaicking noise. In order to solve the above two problems, we shall revisit the DM &DN scheme. Compared to the DN &DM scheme, the advantage of DM &DN is that it does not halve the image size. This is a way around the above mentioned problems. A serious drawback, though, is that chromatic and spatial correlations have been introduced by the demosaicking in the noise of the CFA image. The result is that the noise is no longer white. We next analyze some properties of the demosaicked noise.

Definition Given a ground truth color image (R, G, B) we define the demosaicked noise associated with a demosaicking method DM in the following way: first the image is mosaicked so that only one value of either R, G, DM is kept at each pixel, according to a fixed Bayer pattern. Then white noise with standard deviation σ 0 is added to the mosaicked image, and the resulting noisy mosaic is demosaicked by DM , hence giving a noisy image ( R, G, B). We call demosaicked noise the difference

( R -R, G -G, B -B).
In short, it is the difference between the demosaicked version of a noisy image and its underlying ground truth.

The model of the demosaicked noise depends on the choice of the demosaicking algorithm DM . For the demosaicking step we will evaluate the following state of the art methods, which have an increasing complexity: HA [START_REF] Jr | Adaptive color plan interpolation in single sensor color electronic camera[END_REF], RI [START_REF] Kiku | Residual interpolation for color image demosaicking[END_REF], MLRI [START_REF] Kiku | Minimized-laplacian residual interpolation for color image demosaicking[END_REF], ARI [START_REF] Monno | Adaptive residual interpolation for color and multispectral image demosaicking[END_REF], LSSC [START_REF] Mairal | Non-local sparse models for image restoration[END_REF], RCNN [START_REF] Tan | Color image demosaicking via deep residual learning[END_REF] and JCNN [START_REF] Gharbi | Deep joint demosaicking and denoising[END_REF]. We are interested in algorithms with low or moderate power; only HA, RI, MLRI and RCNN have a reasonable complexity in this context. For the denoising step we shall likewise consider two classic hand-crafted algorithms, CBM3D and nlBayes.

Fig. 2 (c1) and (c2) shows an example where noisy CFA images with noise of standard deviation σ 0 were first demosaicked by RCNN and then restored by CBM3D assuming a noise parameter σ = σ 0 . The output of CBM3D with σ = σ 0 has a strong residual noise. Similar results are also obtained with nlBayes (see the supplementary material). To understand empirically the right noise model to adopt after demosaicking, we simulated this DM &DN pipeline for different levels of noise σ 0 , and applied CBM3D after demosaicking with a noise parameter corresponding to σ 0 multiplied by different factors (1.0, 1.1, • • • , 1.9).

Table 1:

Denoising after demosaicking DM &DN , where DN is CBM3D [START_REF] Dabov | Color image denoising via sparse 3d collaborative filtering with grouping constraint in luminance-chrominance space[END_REF] with noise parameter equal to C σ 0 , while noise in the raw image has standard deviation σ 0 = 20. Each row shows the CPSNR result for C ranging from 1.0 to 1.9. Each column corresponds to a different demosaicking method DM . The best result of each column is marked with a box . The best result of each line is in red and the second best one is in green. The best factor C for all methods is C 1.5, the same is true for different values of σ 0 as well (see supplementary material). The results are shown in Table 1, where the classic color peak signal-to-noise ratio (CPSNR) [START_REF] Alleysson | Linear demosaicing inspired by the human visual system[END_REF] is adopted as a logarithmic measure of the performance of the algorithms. It is defined by

CPSNR = 10 log 10 255 2 X=R,G,B MSE(X)/3 , with MSE(X) = 1 |Ω| (i,j)∈Ω ( X(i, j) -X(i, j)) 2 ,
where X denotes the ground truth image and X is the estimated color image. From 1.0 to 1.9, the CPSNR increases first and then decreases. We can see that the best values are distributed on the lines with factors from 1.4 to 1.7. A similar behavior was also observed using nlBayes for denoising as well as for other levels of noise (see the supplementary material).

This does not mean that the overall noise standard deviation has increased after demosaicking. Let us consider the noise standard deviation estimated as the mean RMSE of the demosaicked images from the Imax [START_REF] Zhang | Color demosaicking by local directional interpolation and nonlocal adaptive thresholding[END_REF] dataset with different noise levels, given in Table 2. We observe that for low noise (σ 0 = 1) there is a serious demosaicking error, of about 4, not caused by the noise, but by the demosaicking itself. However, for σ 0 > 10 we see that the RMSE of the demosaicked image tends to roughly 3/4 of the initial noise standard deviation.

At first sight, this 3/4 factor seems to contradict the observation that denoising with a parameter 1.5σ 0 yields better results. This leads us to analyzing the structure of the residual noise. Fig. 3 shows an image contaminated with AWG noise with standard deviation σ 0 = 20 and its resulting demosaicked noise for respectively HA, MLRI, RCNN. In the last row of the figure, one can observe the color spaces (in standard (R,G,B) Cartesian coordinates) of each of these noises, each cloud being presented in its projection with maximal area. As expected, the AWG color space is isotropic and has an apparent diameter proportional to 4σ 0 80. The color space of the demosaicked noise is instead elongated in the luminance direction Y = R+G+B √ 3 to about 6σ 0 120 and squeezed in the others. This amounts to an increased noise standard deviation for Y after demosaicking, and much less noise in the chromatic directions. This is confirmed by Table 3 that shows variances and covariances of (R, G, B) and (Y, U, V ) respectively for an AWG noise with σ 0 = 20, and then for the demo- Cartesian coordinates) of each noise, presented in their projection with maximal area. As expected, the AWG color space is isotropic, while the color space after demosaicking is elongated in the luminance direction Y and squeezed in the others. This amounts to an increased noise standard deviation for Y after demosaicking, and less noise in the chromatic directions.

saicked noise obtained from it after demosaiciking with RI, MLRI and RCNN. In Table 3 (a) these statistics are computed on a pure white noise image with σ = 20.

Hence the variance of Y is 400, as the (R, G, B) → (Y, U, V ) transform is implemented as an isometry of R 3 . The variance of Y is a growing sequence for the demosaicked noise obtained by increasingly sophisticated demosaicking: 715 for RI, 772 for MLRI, 972 for RCNN. In contrast, the demosaicked noise is reduced in the U and V axes, with its variance passing from 400 for AWGN to 168 and 94 for RI, and even down to 43 and 55 for RCNN. Hence, the noise standard deviation on U or V has been divided by a factor between 2 and 3. But Table 3 also shows that the residual noise on U and V is strongly spatially correlated, it is therefore a low frequency noise, that will require stronger filtering than white noise to be removed. This table also shows that the Y component of the demosaicked noise remains almost white. This leads to a simple conclusion: since image denois-ing algorithms are guided by the Y component [START_REF] Dabov | Color image denoising via sparse 3d collaborative filtering with grouping constraint in luminance-chrominance space[END_REF][START_REF] Lebrun | A nonlocal bayesian image denoising algorithm[END_REF], we can denoise with methods designed for white noise, but with a noise parameter adapted to the increased variance of Y .

To understand why the variance of Y is far larger than the AWG noise it comes from, let us study in Table 4 the correlation between the three channels (R, G, B) in the demosaicked noise of RI, MLRI and RCNN. We observe a strong (R, G, B) correlation ranging from 0.6 for RI to 0.89 for RCNN, which is caused by the "tendency to grey" of all demosaicking algorithms. Assuming that the demosaicked noisy pixel components (denoted R , G , B ) have a correlation coefficient close to 1 then we have

Y = R + G + B √ 3 ∼ √ 3 N (0, σ 0 ).
This factor of about 1.7 corresponds to the case with maximum correlation. Our empirical observation of an optimal factor near 1.5 responds to a lower correlation between the colors.

Experimental evaluation

To evaluate the proposed framework for denoising and demosaicking, we shall use two classic noise free color image datasets: Kodak and Imax. The Imax dataset [START_REF] Zhang | Color demosaicking by local directional interpolation and nonlocal adaptive thresholding[END_REF] consists of 18 images of 500 × 500 pixels, cropped from high-resolution images of size: 2310 × 1814. The Kodak dataset consists of 25 images of 768 × 512 pixels released by the Kodak Corporation for unrestricted research usage 1 . We also evaluated it on a set of 14 real raw images from the SIDD dataset [START_REF] Abdelhamed | A High-Quality Denoising Dataset for Smartphone Cameras[END_REF], which comes with ground truth acquisitions.

Evaluation of DN&DM and DM&DN strategies. We performed simulations with the schemes: DN &DM and DM &DN . The considered demosaicking methods range from classic to very modern: HA [START_REF] Jr | Adaptive color plan interpolation in single sensor color electronic camera[END_REF], RI [START_REF] Kiku | Residual interpolation for color image demosaicking[END_REF], MLRI [START_REF] Kiku | Minimized-laplacian residual interpolation for color image demosaicking[END_REF], ARI [START_REF] Monno | Adaptive residual interpolation for color and multispectral image demosaicking[END_REF], and RCNN [START_REF] Tan | Color image demosaicking via deep residual learning[END_REF]. For the denoising stage two classic hand-crafted patch-based denoising algorithms were considered: CBM3D [START_REF] Dabov | Color image denoising via sparse 3d collaborative filtering with grouping constraint in luminance-chrominance space[END_REF] and nlBayes [START_REF] Lebrun | A nonlocal bayesian image denoising algorithm[END_REF]. As commented in the introduction, both methods can be adapted to handle mosaics (in the DN &DM setting). In 1 Image source: http://r0k.us/graphics/kodak

Table 3:

Variance and covariance of (R, G, B) and (Y, U, V ) (each first row) and the corresponding correlations (each second row) between pixels (i, j) and (i + s, j + t), s, t = 0, 1, 2 first for AWGN (a) with standard deviation σ = 20, then for its demosaicked versions by RI (b), MLRI (c) and RCNN (d). the case of CBM3D this amounts to applying the method by Danielyan et al. [START_REF] Danielyan | Cross-color bm3d filtering of noisy raw data[END_REF], while for nlBayes this is done by denoising the 4-channel image associated to the mosaic. The denoising and demosaicking schemes with the above mentioned demosaicking algorithms and denoising methods were applied to the mosaic images of the Imax image dataset corrupted by additive white Gaussian noise with standard deviations σ 0 = 1, 3, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60. Due to space constraints, in Table 5 we only report the results corresponding to one noise level σ 0 = 20. Results corresponding to other noise levels are in the supplementary material. From Table 5, we can see that DM &DN with parameter σ = σ 0 is not better than DN &DM , but DM &1.5DN (which denotes denoising DN with parameter σ = 1.5σ 0 ) beats clearly DN &DM .This might explain why many researchers think that the scheme DN &DM was superior to the scheme DM &DN . Table 5: Comparison in CPSNR(dB) of average restoration performance between DN &DM and DM &DN for a fixed level of noise σ 0 = 20. We test two denoisers DN namely CBM3D, and nlBayes, and 1.5DN means that if noise level is σ 0 , the noise level parameter for the denoising method DN is σ = 1.5σ 0 . Both denoisers can be adapted to handle mosaics in the DN &DM schemes (see in the text). The best result of each column is marked with a box . The best result of each line is in red and the second best one is in green. In addition to the good CPSNR results, one important advantage of the DM &DN schemes is the high visual quality of the final restored images. Fig. 2 demonstrates the differences between the various solutions (based on BM3D) obtained on the test image number 3 of the Imax dataset with σ 0 = 20. To save space, only crops of the full-color results and corresponding differences with the ground truth are shown here.

(i,j) (i,j+1) (i,j+2) (i+1,j) (i+1,j+1) (i+1,j+2) (i+2,j) (i+2,j+1) (i+2,
The DN &DM scheme shown in Fig. 2 (b1) and (b2) uses BM3D-CFA [START_REF] Danielyan | Cross-color bm3d filtering of noisy raw data[END_REF] for denoising; we can observe some minor checkerboard artifacts. From Fig. 2 (c1) and (c2), we can deduce that there is no checkerboard effect but that much noise remains in the restored image by DM &DN schemes with parameter 1.0σ 0 . The result of DM &1.5DN (Fig. 2 (d1) and (d2)) are smooth without checkerboard effects. Fig. 2 (e1) and (e2) correspond to the outputs of the CNN joint denoising and demosaicking method JCNN [START_REF] Gharbi | Deep joint demosaicking and denoising[END_REF].

One can observe thin structures in the upper left corner of In short, it appears that the DM &DN scheme with an appropriate parameter (namely DM &1.5DN ) outperforms the competition in terms of visual quality. This is due to the fact that it efficiently uses spatial and spectral image characteristics to remove noise, preserve edges and fine details. Indeed, contrary to the DN &DM schemes, DM &1.5DN does not reduce the resolution of the noisy image. Using an DN &DM scheme ends up over-smoothing the result. It comes to no surprise that JCNN performs slightly better than the other methods; however, it is much more computationally demanding and only works for σ ≤ 20.

In a systematic comparison between the schemes involving CBM3D and nlBayes, schemes with CBM3D proved to perform slightly better. Furthermore, the DM &DN schemes with CBM3D are about four times faster than nlBayes. Hence, the following experiments are more focused on CBM3D.

Comparison with methods from the literature. To complete this comparison we went back to all DN &DM schemes proposed in the literature, and performed a systematic comparison for the two classic Kodak and Imax datasets. These databases are always used in demosaicking evaluations, because they illustrate different challenges of the demosaicking problem, Imax being difficult by its color contrast, and Kodak challenging for the recovery of fine structure. In Tables 6 and7 we compare representative DN &DM methods from the literature with the best DM &DN methods identified above (all of them DM &1.5DN ):

-The two best performing demosaicking before denoising methods (DM &1.5DN ) from on Table 5 are considered. Namely, RCNN for demosaicking followed by CBM3D (denoted RCNN+CBM3D) or nl-Bayes (RCNN+nlBayes) for denoising.

-We also consider a "low-cost" DM &1.5DN combination using MLRI [START_REF] Kiku | Minimized-laplacian residual interpolation for color image demosaicking[END_REF] for demosaicking and CBM3D for denoising (MLRI+CBM3D).

The considered DN &DM methods from the literature are:

-The BM3D-CFA filter was proposed in [START_REF] Danielyan | Cross-color bm3d filtering of noisy raw data[END_REF] to avoid the checkerboard effects resulting from independently We compare an DN &DM scheme BM3D+RCNN [START_REF] Danielyan | Cross-color bm3d filtering of noisy raw data[END_REF], with three DM &1.5DN RCNN+CBM3D, RCNN+nlBayes and MLRI+BM3D. As a reference we also include the result of a joint CNN method JCNN [START_REF] Gharbi | Deep joint demosaicking and denoising[END_REF]. But its results are only available for noise with σ ≤ 20 because the network is not trained beyond that level.

applying BM3D to the color phases of CFA images. We evaluate BM3D-CFA [START_REF] Danielyan | Cross-color bm3d filtering of noisy raw data[END_REF] followed by Hamilton Adams demosaicking (denoted BM3D+HA), as well as followed by a state-of-the-art RCNN demosaicking [START_REF] Tan | Color image demosaicking via deep residual learning[END_REF] (BM3D+RCNN).

-The CFA denoising framework of Park et al. [START_REF] Sung | A case for denoising before demosaicking color filter array data[END_REF] effectively compacts the signal energy while the noise is distributed equally in all dimensions by using a color representation from the principal components analysis of the pixel RGB values in the Kodak dataset and then removes noise in each channel by BM3D. This preprocessing is advantageous for the Kodak image set, but inadequate for the Imax image set. We evaluate this framework [START_REF] Sung | A case for denoising before demosaicking color filter array data[END_REF] with BM3D [START_REF] Dabov | Image denoising with block-matching and 3d filtering[END_REF] followed by the RCNN demosaicking [START_REF] Tan | Color image demosaicking via deep residual learning[END_REF] (Park+RCNN).

-The PCA-CFA filter proposed in [START_REF] Zhang | Pca-based spatially adaptive denoising of cfa images for single-sensor digital cameras[END_REF] is a spatiallyadaptive denoising based on principal component analysis (PCA) that exploits the spatial and spectral correlations of CFA images to preserve color edges and details. We evaluate PCA-CFA [START_REF] Zhang | Pca-based spatially adaptive denoising of cfa images for single-sensor digital cameras[END_REF] followed by DLMM demosaiciking [START_REF] Zhang | Color demosaicking via directional linear minimum mean square-error estimation[END_REF] (PCA+DLMM) and RCNN demosaicking [START_REF] Tan | Color image demosaicking via deep residual learning[END_REF] (PCA+RCNN).

-Finally, as a reference, we include the deep learning based joint denoising and demosaicking (JCNN) of [START_REF] Ehret | A Study of Two CNN Demosaicking Algorithms[END_REF][START_REF] Gharbi | Deep joint demosaicking and denoising[END_REF]. But its results are only available for noise with σ ≤ 20 because the network is not trained beyond that level.

From Tables 6 and7 we see that the DM &DN method RCNN+CBM3D as well as RCNN+nlBayes yield the best results on the Kodak dataset, and the margin with respect to the best DN &DM method (BM3D+RCNN, i.e. BM3D-CFA [START_REF] Dabov | Color image denoising via sparse 3d collaborative filtering with grouping constraint in luminance-chrominance space[END_REF] with RCNN [START_REF] Tan | Color image demosaicking via deep residual learning[END_REF]) is quite large: more than 1.5dB on average. In Fig. 3 we compare some results obtained on an image from the Kodak dataset. From the upper-left extract we can see that textures are better restored with RCNN+CBM3D and MLRI+CBM3D, while JCNN introduces some defects. From the extract we see that the DM &1.5DN methods preserve much more details than BM3D+RCNN, and the result is comparable to

JCNN.

On the Imax database RCNN+CBM3D has the highest CPSNRs on high noise levels, by a small gap though. For low noise levels BM3D+RCNN is better, but the difference with RCNN+CBM3D is very small. The joint denoising-demosaicking network JCNN [START_REF] Gharbi | Deep joint demosaicking and denoising[END_REF] yield the best results on the Imax dataset for σ ≤ 20 (not trained above those levels) yet, the margin with respect to RCNN+CBM3D is again small. Overall, by looking at the average CPSNR we can say that the DM &1.5DN scheme RCNN+CBM3D is indeed much more robust than BM3D+RCNN.

Evaluation on real images. We evaluated on a set of 14 raw images taken from the Small SIDD dataset [START_REF] Abdelhamed | A High-Quality Denoising Dataset for Smartphone Cameras[END_REF]. For simplicity, the selected images correspond to phones from the same manufacturer. We adopted the simple pipeline proposed by the authors, which yields photo finished images that can be compared with the ground truth. The considered methods (RCNN+CBM3D, CBM3D+RCNN, and JCNN) were applied at the demosaicking stage (in linear space). Before any denoising step (DN ) we applied a VST (squared root [START_REF] Anscombe | The Transformation of Poisson, Binomial and Negative-Binomial Data[END_REF]), which whitens the noise, and invert it afterwards. The noise level was estimated using [START_REF] Colom | Analysis and Extension of the Ponomarenko et al. Method, Estimating a Noise Curve from a Single Image[END_REF] and provided to the denoising algorithms and JCNN. Table 8 reports the average CPSNR obtained on these images and the average of the estimated noise levels (after whitening). These values are consistent with the simulated results obtained on the Kodak database (Table 7). The result in Fig. 5, and the supplementary material, support the case in favor of the DM &1.5DN schemes (RCNN+CBM3D).

Conclusions

This paper analyzed the advantages and disadvantages of denoising before demosaicking (DN &DM ) schemes, versus demosaicking before denoising (DM &DN ), to recover high quality full-color images. We showed that for the DM &DN schemes a very simple change of the noise parameter of the denoiser DN coped with the structure of demosaicked noise, and led to efficient denoising after demosaicking. We found that, this allowed to preserve fine structures that are often smoothed out by the DN &DM schemes. Our best performing combination in terms of quality and speed is a DM &1.5DN scheme, where demosaicking DM is done by a fast algorithm RCNN [START_REF] Tan | Color image demosaicking via deep residual learning[END_REF] followed by CBM3D denoising 1.5DN with noise parameter equal to 1.5σ 0 .

Nevertheless it seems ineluctable to see deep learning win the end game when solutions will be found to have more compact or more rapid joint demoisaickingdenoising algorithms.

Figure 1 :

 1 Figure 1: Bayer color filter array, CFA, which is used by most cameras.

  (a1) Ground truth (a2) Ground truth (b1)DN &DM /26.92dB (b2)DN &DM /26.92dB (c1) DM &DN /25.38dB (c2) DM &DN /25.38dB (d1) DM &1.5A /26.95dB (d2) DM &1.5A /26.95dB (e1) JCNN /27.46dB (e2) JCNN /27.46dB

Figure 2 :

 2 Figure 2: Comparison of different denoising and demosaicking schemes with noise σ 0 = 20.In each experiment: left, detail of the demosaicked and denoised image; right, the difference with original that should contain mainly noise. DN : CBM3D denoising; DM : demosaicking (here we use RCNN). 1.5DN means that if noise level is σ 0 , the input noise level parameter of denoising method DN is σ = 1.5σ 0 ; DN &DM : uses the BM3D-CFA framework[START_REF] Danielyan | Cross-color bm3d filtering of noisy raw data[END_REF] for denoising.

Figure 3 :

 3 Figure 3: AWG noise image and demosaicking noise with standard deviation σ = 20 for respectively HA, MLRI, RCNN. Last row: the color spaces (in standard (R,G,B)Cartesian coordinates) of each noise, presented in their projection with maximal area. As expected, the AWG color space is isotropic, while the color space after demosaicking is elongated in the luminance direction Y and squeezed in the others. This amounts to an increased noise standard deviation for Y after demosaicking, and less noise in the chromatic directions.

Fig 2 (

 2 a1), but they disappear in the restored image by DN &DM . The proposed DM &1.5DN scheme restores them. The second column of Fig 2 illustrates a similar situation in which thin details are recovered by DM &DN and DM &1.5DN but not in the others.

Figure 4 :

 4 Figure4: Demosaicking and denoising results on an image from the Kodak dataset with σ = 20. We compare an DN &DM scheme BM3D+RCNN[START_REF] Danielyan | Cross-color bm3d filtering of noisy raw data[END_REF], with three DM &1.5DN RCNN+CBM3D, RCNN+nlBayes and MLRI+BM3D. As a reference we also include the result of a joint CNN method JCNN[START_REF] Gharbi | Deep joint demosaicking and denoising[END_REF]. But its results are only available for noise with σ ≤ 20 because the network is not trained beyond that level.

Figure 5 :

 5 Figure 5: Details of a real images (enhanced contrast) from the SIDD [1] dataset. From left to right: noisy input (demosaicked), BM3D+RCNN , and RCNN+CBM3D.

C

  HA GBTF RI MLRI ARI LSSC RCNN 1.0 28.15 27.58 28.46 27.95 28.70 27.19 27.28 1.1 28.56 28.15 28.83 28.44 28.98 27.89 28.05 1.2 28.85 28.55 29.08 28.80 29.18 28.43 28.67 1.3 29.05 28.81 29.23 29.03 29.29 28.78 29.09 1.4 29.18 28.96 29.31 29.17 29.35 29.00 29.34 1.5 29.23 29.00 29.32 29.22 29.35 29.06 29.41 1.6 29.25 29.01 29.30 29.23 29.33 29.06 29.41 1.7 29.25 28.97 29.26 29.20 29.29 29.02 29.36 1.8 29.22 28.92 29.20 29.15 29.23 28.95 29.28 1.9 29.17 28.85 29.13 29.08 29.17 28.88 29.20

Table 2 :

 2 RMSE between original and demosaicked image for different demosaicking algorithms in presence of noise of standard deviation σ 0 .

	σ0 HA GBTF RI MLRI ARI LSSC RCNN
	1 5.04 5.10 4.17 4.06 3.72 4.40	3.21
	5 6.78 6.87 6.12 6.10 5.74 6.36	5.59
	10 10.18 10.27 9.53 9.74 9.09 9.96	9.65
	20 17.75 17.83 16.77 17.56 16.06 18.16 18.04
	40 32.67 32.76 30.77 32.64 29.36 33.68 33.98
	60 46.14 46.35 43.43 46.11 41.44 48.11 47.95

Table 4 :

 4 Covariances (each first row) and correlations (each second row) of the three color channels (R, G, B) of the demosaicked noise, when the initial CFA white noise satisfies σ 0 = 20

		R	G	B		R	G	B
	R	336.44 206.29 175.01 1.0000 0.6542 0.5097	R	361.42 224.39 201.41 1.0000 0.6826 0.5501
	G	206.29 295.54 200.96 0.6542 1.0000 0.6244	G	224.39 298.94 216.86 0.6826 1.0000 0.6512
	B	175.01 200.96 350.46 0.5097 0.6244 1.0000	B	201.41 216.86 370.92 0.5501 0.6512 1.0000
			(a) RI				(b) MLRI
		R	G	B		R	G	B
	R	359.90 320.44 302.85 1.0000 0.8967 0.8461	R	334.84 297.31 275.28 1.0000 0.8675 0.8181
	G	320.44 354.83 299.85 0.8967 1.0000 0.8437	G	297.31 350.81 270.32 0.8675 1.0000 0.7848
	B	302.85 299.85 355.99 0.8461 0.8437 1.0000	B	275.28 270.32 338.17 0.8181 0.7848 1.0000
			(c) RCNN				(d) JCNN

  .[START_REF] Monno | Adaptive residual interpolation for color and multispectral image demosaicking[END_REF] 27.97 28.69 27.27 DM &DN 28.15 28.46 27.95 28.70 27.28 DM &1.5DN 29.24 29.32 29.22 29.36 29.41 nlBayes DN &DM 28.17 28.17 28.17 28.18 28.28 DM &DN 28.67 28.99 28.57 29.21 28.02 DM &1.5DN 29.29 29.26 29.22 29.31 29.36

	DN Algorithm	HA	RI MLRI ARI RCNN
	CBM3D	DN &DM	28.11 28

Table 6 :

 6 Comparison of the results (CPSNR in dB) between different denoising and demosaicking methods for the Imax image set. The best result of each line is in red, the second best one is in green and the third best one is in blue. 34.63 38.53 35.37 33.99 37.52 38.36 38.42 36.52 38.59 5 33.43 35.62 32.86 32.69 34.87 35.39 35.29 34.60 33.48 10 31.84 32.92 30.06 30.73 31.89 32.75 32.59 32.36 33.09 20 29.22 29.55 26.86 27.57 27.99 29.41 29.25 29.22 29.79 40 25.50 25.51 23.86 23.50 23.57 25.52 25.09 25.39 -60 21.55 21.34 21.75 20.89 20.89 22.78 22.31 22.63 -Av 28.09 28.88 26.89 26.71 27.53 28.99 28.72 28.58 -

				DN &DM				DM &1.5DN
		BM3D BM3D Park PCA PCA RCNN RCNN MLRI
	σ	+	+	+	+	+	+	+	+	JCNN
		HA RCNN RCNN DLMM RCNN CBM3D nlBayes CBM3D

1

Table 7 :

 7 Comparison of the results (CPSNR in dB) between different denoising and demosaicking methods for the Kodak image set. The best result of each line is in red, the second best one is in green and the third best one is in blue. 34.70 40.55 40.36 38.19 39.12 40.98 40.98 38.52 41.15 5 32.84 34.89 34.87 34.99 35.42 36.55 36.42 35.71 34.13 10 30.34 30.93 30.85 31.83 32.01 33.36 33.18 32.94 33.27 20 27.59 27.70 27.42 28.11 28.14 29.98 29.87 29.70 29.95 40 24.79 24.78 24.88 24.15 24.08 26.71 26.29 26.44 -60 22.58 22.55 23.19 21.77 21.70 24.42 23.93 24.16 -Av 27.47 28.35 28.36 27.96 28.09 30.19 29.93 29.64 -

				DN &DM				DM &1.5DN
		BM3D BM3D Park PCA PCA RCNN RCNN MLRI
	σ	+	+	+	+	+	+	+	+	JCNN
		HA RCNN RCNN DLMM RCNN CBM3D nlBayes CBM3D

1

Table 8 :

 8 Average CPSNR over 14 raw images taken from the Small SIDD dataset[START_REF] Abdelhamed | A High-Quality Denoising Dataset for Smartphone Cameras[END_REF]. The reported average noise level is scaled to the range 0-255.

	mean σ	CBM3D+RCNN RCNN+CBM3D	JCNN
	7.65	38.19	39.64	38.54

Acknowledgments : Work partly financed by Office of Naval research grant N00014-17-1-2552 and DGA Astrid project n • ANR-17-ASTR-0013-01.


Eloi Zalczer

François-Xavier Thomas

Laurent Chanas

Gabriele Facciolo

Frédéric Guichard

Depth Map Quality Evaluation for Photographic Applications

As depth imaging is integrated into more and more consumer devices, manufacturers have to tackle new challenges. Applications such as computational bokeh and augmented reality require dense and precisely segmented depth maps to achieve good results. Modern devices use a multitude of different technologies to estimate depth maps, such as time-of-flight sensors, stereoscopic cameras, structured light sensors, phase-detect pixels or a combination thereof. Therefore, there is a need to evaluate the quality of the depth maps, regardless of the technology used to produce them. The aim of our work is to propose an end-result evaluation method based on a single scene, using a specifically designed chart. We consider the depth maps embedded in the photographs, which are not visible to the user but are used by specialized software, in association with the RGB pictures. Some of the aspects considered are spatial alignment between RGB and depth, depth consistency, and robustness to texture variations. This work also provides a comparison of perceptual and automatic evaluations.

Introduction

Computational Bokeh, and various other depth-based photographic features, have become major selling points for flagship smartphones in the last few years. These applications have different quality requirements than most depth imaging technologies. Bokeh simulation requires depth maps to be very precisely aligned with an RGB picture, and suffers largely from artifacts and segmentation errors. On the other hand, distance estimation is generally less important, and the relative depth levels of the elements of the scene is what matters. Additionally, those applications use various types of depth imaging technologies, ranging from single-image deep learning estimation to time-of-flight sensors.

In this article, we present a quality evaluation framework for depth maps as a means towards photographic applications. The different aspects of quality that we evaluate can be classified in two main categories: segmentation and depth consistency. Examples of segmentation problems are offset edges or incomplete detection of small elements. Under depth consistency, we regroup aspects such as sensitivity to texture variations, depth estimation and depth regularity.

Our approach is based on a single scene, designed to highlight those defects and make them measurable. The scene contains a 2-dimensional chart and a background. Our measurement algorithms output several metrics, mostly local. Rather than a single global evaluation metric, we propose feature-level measurements to enable manufacturers to evaluate precisely the strengths and weaknesses of their devices.

The rest of the article is organized as follows. The next section presents the chart that we designed as the main element of our test setup. The following section introduces the experimen-Figure 1. Picture of the bokeh setup proposed in [START_REF] Hauser | Image quality benchmark of computational bokeh[END_REF] and the associated depth map captured with a smartphone. Depth map defects are obviously visible in this example, mostly around the head and in the edges. Some background elements start to appear, and some holes of the palettes are blended in with the foreground.

tal setup and test protocol used for our analyses. Afterwards, we present the metrics we propose as well as the implemented algorithms, followed by the analysis of the results and a comparison with perceptual evaluation.

Related Work

A number of papers on depth map quality evaluation have been published over the last few years. Several trends can be distinguished. Haddad [START_REF] Haddad | Non-reference depth map quality evaluation in immersive video applications[END_REF] and Kim et al. [START_REF] Kim | Depth map quality metric for three-dimensional video[END_REF] focus on depth maps for application such as Free Viewpoint Television or 3D Television, whereas Koch et al. [START_REF] Koch | Evaluation of cnn-based single-image depth estimation methods[END_REF] proposes a set of metrics to evaluate the output of CNN-based single-image depth estimation. Several papers [START_REF] Banitalebi-Dehkordi | A study on the relationship between depth map quality and the overall 3d video quality of experience[END_REF][START_REF] You | Perceptual quality assessment for stereoscopic images based on 2d image quality metrics and disparity analysis[END_REF] propose evaluations based on well-known 2D quality metrics (PSNR, SSIM, VIF, etc.) applied to stereoscopic cases. Others, such as Benoit et al. [START_REF] Benoit | Quality assessment of stereoscopic images[END_REF], propose local quality metrics for disparity distortion. Finally, many papers [START_REF] Eigen | Depth map prediction from a single image using a multi-scale deep network[END_REF][START_REF] Liu | Deep convolutional neural fields for depth estimation from a single image[END_REF][START_REF] Li | Depth and surface normal estimation from monocular images using regression on deep features and hierarchical crfs[END_REF][START_REF] Wang | Towards unified depth and semantic prediction from a single image[END_REF][START_REF] Roy | Monocular depth estimation using neural regression forest[END_REF] also propose neural network architectures to extract depth maps from single images, using various quality metrics for training.

Those works have in common that they do not take advantage of a complete prior knowledge of the scene. They are meant to be usable for image-to-image comparisons with a ground truth image, which makes them easily testable against datasets such as KITTI [START_REF] Geiger | Vision meets robotics: The kitti dataset[END_REF], NYU v2 [START_REF] Silberman | Indoor segmentation and support inference from rgbd images[END_REF] or Middlebury [START_REF] Scharstein | A taxonomy and evaluation of dense two-frame stereo correspondence algorithms[END_REF]. The novelty of our work is to use a precisely known laboratory setup to compute local metrics, targeted at specific known defects.

Chart design

The main part of our experimental setup is the chart. Our goal was to design the elements of the chart to mimic some real life situations, and highlight some specific defects. Therefore, its conception was very much based on observation and experimental prototyping. As starting point we used the Bokeh chart proposed in [START_REF] Hauser | Image quality benchmark of computational bokeh[END_REF] shown in Fig. 1. This work was designed for perceptual evaluation, and was not fit for automated measurements because of the imperfect knowledge of the geometry of the chart. However, it enabled us to observe and classify many defects.

From those observations, we established a list of interesting elements to put in our new chart. Starting with textures, it was obvious that strong texture variations caused some devices to infer depth variations, even if there was none. This is likely caused by some kind of guided filtering [START_REF] Hosni | Fast cost-volume filtering for visual correspondence and beyond[END_REF], which causes background details to appear in the depth maps. Furthermore, highly repetitive patterns are tricky for devices that use stereoscopic vision, and low-reflectivity surfaces usually cause depth estimation issues for time-of-flight or structured light sensors.

Segmentation issues were most visible around elements that are very close together, such as fingers. Artifacts were also visible near the end of thin elements, where the resolution of the depth sensors may not be sufficient and filtering algorithms come into play. In low-contrast areas, the difference between foreground and background was sometimes not detected, highlighting the importance of the background.

Knowing that face detection is almost universally used in the process of depth estimation in smartphones, a simulated head is an essential element to put in our setup. It is likely that some devices will infer the shape of a body under the head, introducing some interesting defects. Furthermore, the face is the part of the body where human perception focuses the most and where any small defect can result in a very unrealistic portrait.

The decision was made to create a two-dimensional chart, shown in Fig. 2, because measurements can be carried out much more easily than on a three-dimensional one. It is easier to know the geometry with sufficient precision, and avoids occlusion and perspective issues. The legend for the regions of interest in Fig. 2 are given in Tab. 1. The bottom part of the chart is simply used as a support, so that it can be held without wasting any useful space. The chart was printed on DIBOND R aluminium composite material, which makes it rigid and reasonably flat.

Experimental setup

Our experimental setup consists of the chart and a background, providing two layers of depth. As previously stated, the importance of the background is capital. The contrast and the patterns of the background have a big impact on the end result, especially parts that are connected or very close to the foreground elements. Therefore, the choice of the background to be used was None, used to check the printer alignment Table 1. Elements of the chart non-trivial. Our first idea was to use a uniform white backdrop, but this was unfair to stereoscopic devices, and did not really resemble any real life situation. In the end, our choice was to use the same background that is used for the bokeh setup of [START_REF] Hauser | Image quality benchmark of computational bokeh[END_REF].

This background is quite challenging because it has lots of small elements, however it is fair for every kind of technology, which is what we are mainly looking for. It is printed on a backdrop and is therefore flat. In order to ensure that results are comparable among devices, the same exact framing is used in every case. During a shooting session, both the distance from the chart to the background and the distance from the device to the chart need to be adjusted depending on the focal length of the device, using the following formula

d 2 = f 2 f 1 d 1 ,
where d 2 is the new distance between the device and the element, d 1 the previous distance, and f 1 and f 2 are the focal lengths of each device.

We use two lighting conditions during our shooting sessions: D65 at 1000 lux, simulated by a fluorescent illuminant with a color temperature of 6500K, and A at 50 lux, a tungsten illuminant at 2856K. Those two cover the most common illuminants and offer widely different spectra. In both cases, our shooting protocol consists of five pictures taken on a tripod, forcing the autofocus to run again between pictures. The chart is placed parallel to the sensor of the device or at low angles (≤ 10 • ).

Metrics and algorithms

The goal of our work was to propose a number of local metrics rather than one global metric. In every case, we started from a known defect, and tried to find a metric to measure it. All of our metrics take advantage of a near-perfect knowledge of the scene to compute relevant values.

The linearity issue

As a starting point, we need a precise estimate of the ground truth of the scene. By considering that the chart and the background are both flat, the depth becomes easy to model. However, the interpretation of the depth maps is more problematic. There is currently no standard for depth map encoding, and the formats used by the manufacturers are usually not documented. Since the absolute distance is not important in this use case, depth maps are not directly associated to real distances. The conversion from a distance to a gray level is done by an unknown function, which we assume strictly monotonic and continuous. In order to measure the accuracy, we need to project the ground truth in the depth map space, and therefore estimate this function.

The four markers in the corners of the chart enable us to know its position in the RGB image with an error of less than two pixels, even in noisy conditions. It is therefore possible to extract the foreground pixels of the depth map and fit a parametric model to the surface. However, a trade-off needs to be found : an overly complicated model would risk over-fitting, and we would not be able to detect defects on the surface ; on the other hand, another extreme solution would be to consider that the chart is perfectly parallel to the sensor and use a single depth value for the chart, but this would leave too much room for human error during the shooting sessions.

A reasonable solution would have been to try several parametric models (e.g. linear, logarithmic, inverse...), but we observed that the variations were too small and noisy to extract any meaningful tendency since the chart was almost parallel to the sensor. As an illustration, Fig. 3 displays a plot of the gray levels along the width of the chart. Therefore, we chose a simple linear model for this fitting.

To reduce the impact of outliers, we use the robust Huber regression distance [START_REF] Huber | Robust estimation of a location parameter[END_REF]. This parametric model enables us to infer an expected depth value for every foreground pixel. That way, we can compute a difference between the expected background and foreground values for every pixel, and use this difference to adapt the thresholds in our measures. We will refer to this difference as the pixel dynamic. When computing metrics on areas, we use the maximal pixel dynamic in the area and refer to it as area dynamic. Depending on the encoding, the value of the foreground pixels may be lower (darker) than the background, or the opposite : henceforth, we shall assume that the lowest depth values correspond to positions closer to the camera.

The only global metric of this work is called error areas, and corresponds to the proportion of the pixels where the difference between the expected value and the actual value is superior to 10% of the pixel dynamic. It is computed similarly for the background and the foreground, on a cropped area to avoid errors due to framing. The second output of this measurement is an error map (Fig. 4), which enables the user to see at a glance where the problematic areas are.

In addition to this metric, we introduce the planarity error. It is basically the same metric, applied to a specific part of the chart of Fig. 2. Those areas are the face, the checkerboard (seen in Fig. 5) and the two textured strips, enabling us to quickly see which part is most challenging for the device.

Edge segmentation

We define an ideal edge as a sharp transition between two depth levels corresponding to the background and foreground. This transition should happen on a single pixel and be perfectly aligned with the ideal position of the edge. The first step of our process is to apply an homography on the depth map, using bilinear interpolation, to re-align the four markers of the target, which means that the edges of the chart are perfectly parallel to the edges of the image. We chose to define twelve regions of interest along the edges of the chart (shown in Fig. 6), three on each side.

Gradient-based measurements. The first half of our metrics are gradient-based metrics, aimed at spatial characteristics of the edges. The preliminary step to these metrics is to convolve the depth map of the region of interest with a Sobel [START_REF] Kanopoulos | Design of an image edge detection filter using the sobel operator[END_REF] filter. This will allow to precisely compute the position of the edges and the value of the edge gradient. Starting from there, we define the following three metrics. In all the following formulas, n stands for the length of the edge and x edge is the position of the depth map edge along the perpendicular axis. This position should ideally be constant, because the ideal edges are aligned with the edges of the image. x t is the ideal position of the edge, and is constant.

The Pixel Shift (PS) metric is the average distance between the detected edge and the ideal edge. It is expressed as PS = 1 n+1 ∑ n i=0 x edge (i) -x t . This metric, being an average, does not tell us if the shift is global or if the average is biased by a few outliers. This is the purpose of our next metric, which consists of the standard deviation of the edge position (σ edge ), defined as σ edge =

1 n+1 ∑ n i=0 x edge (i) -x edge 2 .
Finally, our third spatial metric is the Relative Gradient (RG). It is a measurement of the edge gradient value, rather than its position. As explained before, the expected behavior is a sharp transition between the background value and the foreground value. The gradient at each point of the edge should then be equal to the pixel dynamic. The value computed is an average result over the region of interest, expressed as a proportion of the pixel dynamic. Fig. 7 gives a visual representation of some of the proposed metrics.

Histogram-based measurements. The second half of our metrics are based on the distribution of depth values in the region of interest of the edges, previously defined. They complement the first three ones with a fine analysis of the depth levels in the entire region of interest, rather than just along the main edge. Their purpose is to detect defects such as background details appearing in the depth map or depth steps effects.

The first metric we define is the Pixel Repartition (PR) metric. We first find the two main peaks of the histogram of depth values, corresponding to the background and foreground values, and then multiply the histogram with a bimodal gaussian distribution centered on the peaks. The PR metric is the sum of this weighted histogram, divided by the total number of pixels in the area. The gaussian windows have a height of 1, which means that the metric cannot exceed 1, and their standard deviations are normalized proportionally to the dynamic. A visual representation is given in Fig. 8.

The use of a Gaussian window makes it so that only a perfect device can obtain a perfect PR of 1, but devices that show slight depth variations will not be too penalized. However, the larger the variations, the lower the PR will be.

This metric works well as a single indicative value, however it does not show which part of the depth map is most problematic. As a complement, we also compute the background and foreground standard deviation in the region of interest, using the average of the two main peaks as the separation threshold.

The last metric is the base-2 entropy of the depth histogram, defined as H = -∑ k p k • log 2 (p k ), where p k is the probability corresponding to the k-th bin of the normalized histogram. This has the advantage of showing whether all pixels have the same or different depth values, thus enabling the user to differentiate a step effect and a gradient effect.

Holes detection and dynamic

The main metric regarding the holes of the chart is the Dynamic Proportion (DP). Since the background of our experimental setup is flat, the expected behavior is that the value in every hole should be equal to the background value around the chart. However, due to filtering issues and imprecision of depth sensors, this is often not the case in real depth map images. This metric is computed as follows on a given region: DP = max i, j D i, j -Dfg Dbg -Dfg , with D i, j the value of the depth map at coordinates i, j, Dfg the maximum estimated depth of the chart using the parametric model defined earlier, and Dbg the depth of the background. Using the parametric model makes the metric robust to defects that could appear inside or around the hole. In some cases, holes are detected closer rather than farther by devices, in which case the metric may be negative.

For gradual holes (see Tab. 1), we implement the Contour Coverage (CC), which measures the detected area of the hole as a proportion of the ideal area. Depending on the size of the hole and the background elements behind it, it is frequent that devices only detect some part of the holes. The threshold used is the average between the expected background value and foreground values. The metric is the proportion of pixels in the hole area whose value is greater than the threshold. An example is given in Fig. 9.

For the two triangular holes, this metric would not be as interesting because their area is much larger. The most interesting aspect of those is that their tip is often wrongly filtered out, as shown in the example in Fig. 10. Hence, we define the Visible Height metric, which is computed as the proportion of visible pixels along the bisector of a triangle.

Tab. 2 summarizes the proposed metrics along with the targeted defects. All of the metrics are repeatable, with no stochastic aspect. They are also invariant to the resolution of the depth maps. The metrics that are based on the dynamic are normalized, making them invariant to the encoding of the depth. Filtering issues and hardware limitations Table 2. Implemented metrics and targeted defects.

Validation and results

Computing the correlation coefficient between our metrics and a perceptual evaluation is difficult because our metrics are local and focus on very specific features of the depth maps. Our work provides building blocks, upon which aggregated metrics could be designed to match human perception. In its current state, this is more intended to be a tool for the tuning devices rather than a perceptual quality evaluation. Nevertheless, the obtained results are very promising. The strengths and weaknesses of the different technologies are easily visible, as in the example of Fig. 5, which was created by a stereoscopic device. Many of the artifacts caused by filtering algorithms are also easily measurable. Examples are given in Tab. 3.

Our measures also work when the chart is not perfectly parallel to the sensor. Our tests included pictures with a target tilted up to 10 • around the vertical axis. On some devices, we observed a smoothing of depth levels on tilted surfaces, revealing a heavy filtering. We also see different segmentation issues and artifacts.

The only metric in our set which was easily correlated with perceptual measurements was the error areas, because it is global. We used the aggregated results of the perceptual analysis performed on the Bokeh chart as reference. Using seven uniformly distributed devices in our database, we find a 78 % correlation. The monotonicity is overall preserved, and we can expect a top device to have a low value for the error areas. However, our metric is more sensitive to small depth variations than perceptual measurements, which explains some outliers. Moreover, the perceptual evaluation takes into account the Bokeh simulation output, which does not entirely rely on the depth map.

Conclusion

We have presented a full measurement protocol and objective metrics to evaluate the quality of depth maps for photographic applications. Our method is based on a single scene, containing a specifically designed chart and a backdrop. It is designed to be challenging for all commonly used technologies and aims to reproduce and measure some common defects. Eleven metrics Linearity is very good, while RG is lower because the transition is unsharp.

Obvious stepping effect, causing a high σ edge and PS, and a low RG.

The depth transition is barely visible, causing a very low RG. However, linearity and alignment are fairly good.

The defect here is poorly detected because the maximal gradient is linear even though the depth level varies. This is an edge case. Table 3 are used (mostly local) to characterize precisely defined regions of interest. Those regions of interest are located on the edges of the chart, on flat surfaces or around holes. In total, 148 values are computed for each image, providing a very fine-grained analysis.

Considering that we work on the depth maps, which are not visible to the user, our work is oriented to provide a tool for tuning of and for selecting hardware and software components. In the case of computational bokeh, the background blur simulation rely on different parameters that affect the quality of the bokeh independently of the depth map. This explains the difficulty to correlate our metrics with perceptual evaluation. However, with respect to the depth maps, the results are relevant.

The proposed metrics are designed to be easily understandable, and could be used as building blocks in future works to extract higher level metrics. In the future, the possible adoption of an uniformized depth map format [START_REF]Encoding depth and confidence[END_REF] would facilitate its comparison across devices.

Figure 2 .

 2 Figure 2. Annotated depth map chart.ElementExpected defects 1. Triangular holes Filtering issues, wrong detection of thin parts 2. Gradual holes No detection or mixed detection of thin holes 3. Checkerboard Depth estimation issues due to repetitive pattern 4. Textured strip Depth inconsistency due to texture variation 5. Black and white strip Depth inconsistency due to reflectivity variation 6. Head Depth inference caused by face detection 7. Letters Filtering issues because of very contrasted elements 8. Markers None, used for chart detection 9. Crown Filtering issues, wrong detection of thin parts 10. Alignment marker None, used to check the printer alignment Table 1. Elements of the chart

Figure 3 .

 3 Figure 3. Evolution of the gray levels along the width of the chart. We can see that the levels are very noisy and it is difficult to extract a clear tendency.A linear regression is given for reference.

Figure 4 .

 4 Figure 4. Error map displaying the error for each pixel, relative to the pixel dynamic.

Figure 5 .

 5 Figure 5. Artifact due to a high frequency pattern.

Figure 6 .

 6 Figure 6. Regions of interest for the measurements on the edges (in red).

Figure 7 .

 7 Figure 7. Graph of the gray levels along a line orthogonal to an edge of the chart, with visual representations of the relative gradient and pixel shift metrics.

Figure 8 .

 8 Figure 8. Histogram of a region of interest (in orange) and the Gaussian windows (in blue). We can see a small peak on the left, caused by some artifact.

Figure 9 .

 9 Figure 9. Visualisation of the gradual holes of the chart, with visible depth estimation defects. The values for Dynamic Proportion (DP) and Contour Coverage (CC) are given in the table.

Figure 10 .

 10 Figure 10. Example of incomplete triangular hole, with the ideal shape and its bisector.

  estimation inconsistencies, on surfaces, artifacts due to high frequency patterns or texture variations Pixel shift (PS) Misalignment of depth map and image, wrong detection of edges Edge standard deviation Alteration of the shape of contours, due to background elements or filtering Relative gradient edges, sometimes added by devices to blur potential defects Dynamic proportion (DP) Wrong depth estimation in small discontinuities Contour coverage (CC) Incomplete or no detection of small discontinuities Visible height

  Region of interest Results σ edge = 1.72 px RG = 82.1 % PR = 85.0 % σ edge = 0.52 px RG = 51.1% PR = 87.2% H = 4.33 bits σ edge = 20.84 px RG = 36.0% PS = 18.43 px σ edge = 3.45 px RG = 19.8% PS = 2.45 σ edge = 3.12 px RG = 78.4% PS = 1.66 px PR = 56.9 % Remarks Sharp transition. The σ edge is slightly higher because of shifts at the top and bottom ends of the edge.

  . Examples of results obtained for multiple regions of interest of the edges. Only the most relevant metrics are displayed. The regions of interest considered are the areas enclosed in a white rectangle, the corresponding areas in the reference image are enclosed in magenta.


Roland Akiki

Roger Marí

Carlo De Franchis

Jean-Michel Morel

Gabriele Facciolo

Robust Rational Polynomial Camera Modelling for SAR and Pushbroom Imaging

The Rational Polynomial Camera (RPC) model can be used to describe a variety of image acquisition systems in remote sensing, notably optical and Synthetic Aperture Radar (SAR) sensors. RPC functions relate 3D to 2D coordinates and vice versa, regardless of physical sensor specificities, which has made them an essential tool to harness satellite images in a generic way. This article describes a terrain-independent algorithm to accurately derive a RPC model from a set of 3D-2D point correspondences based on a regularized least squares fit. The performance of the method is assessed by varying the point correspondences and the size of the area that they cover. We test the algorithm on SAR and optical data, to derive RPCs from physical sensor models or from other RPC models after composition with corrective functions.

Introduction

Developing a remote sensing application requires a set of tools, one of which is geolocation. Geolocation relates the 3D world coordinates to the 2D image. This is represented by means of a projection function P : R 3 → R 2 , that maps 3D points to the image plane, and its inverse, the localization function L : R 2 × R → R 3 . When all the physical phenomena and components involved in the acquisition process are known, the geolocation functions can be defined by a chain of operations that model such factors, in what is known as a physical or rigorous sensor model.

Pushbroom scanners are the most common optical satellite image acquisition system, typically consisting of a single line of pixel sensors mounted on a platform that captures each line of the image at a different moment in time. As a result, the exterior orientation parameters, i.e. the perspective center and the attitude angles, change from line to line. The intrinsic parameters (e.g. pixel size, focal length, lens distortion), related to the physical design of the sensor, are constant across the image [START_REF] Grodecki | IKONOS stereo feature extraction -RPC approach[END_REF]. A detailed description of a simplified physical sensor model for pushbroom scanners can be found in [START_REF] De Franchis | Attitude refinement for orbiting pushbroom cameras: a simple polynomial fitting method[END_REF]. In the case of Synthetic Aperture Radar (SAR) images, the most used physical sensor model is the Range-1 arXiv:2102.13423v1 [eess.IV] 26 Feb 2021 Doppler model detailed in [START_REF] Curlander | Location of spaceborne SAR imagery[END_REF]. SAR satellites send an electromagnetic wave that is reflected on the ground. The image is acquired line by line (similar to a pushbroom system), and the position of a ground patch in the image is related to its distance to the sensor, known as the range. The Range-Doppler model is constructed based on ephemeris data (time, position and velocity samples along the orbit) and the acquisition timing information. The ephemeris data needs to be interpolated to obtain continuous geolocation functions along the orbit.

Image vendors have adopted the generic Rational Polynomial Camera1 (RPC) model to save customers from having to deal with the complex specificities of rigorous sensor models. The RPC model is independent of physical properties and offers flexibility to work with different coordinate systems. RPCs have become essential metadata to process satellite images in a generic way, from different sources and for multiple tasks, e.g. photogrammetry and radargrammetry based 3D reconstruction or image ortho-rectification and coregistration.

In this article, we describe a terrain-independent algorithm to fit a RPC model from a physical sensor model or any other geolocation model. Our contributions are:

-An open-source implementation of the method as an easy-to-use Python package, which is available at https://github.com/cmla/rpcfit.

-An evaluation of the algorithm's precision and robustness based on real scenarios. We test our method using Sentinel-1 and WorldView-3 images, to fit a SAR physical sensor model or correct an existing RPC model by composing it with a complementary transformation.

Mathematical formulation of the RPC model

The RPC model defines the projection function P as

r n = a (X n , Y n , Z n ) b (X n , Y n , Z n ) c n = e (X n , Y n , Z n ) f (X n , Y n , Z n ) , (1)
where a, b, e, f are cubic polynomials. X, Y, Z represent the longitude, latitude and height of a 3D point; and r, c are the row and column of its projection on the image plane. Equation 1 uses normalized coordinates for better numerical stability, hence the subscript n. Normalized values are in the range [-1, 1] and they are obtained using two scalars, an offset and a scale factor: X n = (X -X offset )/X scale , where X could be r, c, X, Y or Z from Equation 1.

Each RPC polynomial p is defined by 20 coefficients as

p (X, Y, Z) = p 0 + p 1 Z + p 2 Y + p 3 X + p 4 ZY + p 5 ZX +p 6 Y X + p 7 X 2 + p 8 Y 2 + p 9 Z 2 + p 10 ZY X +p 11 Z 2 Y + p 12 Z 2 X + p 13 Y 2 Z + p 14 Y 2 X +p 15 ZX 2 + p 16 Y X 2 + p 17 Z 3 + p 18 Y 3 + p 19 X 3 , (2)
where p i is the i-th coefficient of p. Since we set p 0 = 1 for the RPC denominator polynomials, a total of 78 coefficients need to be determined to define a, b, e and f in Equation 1.

Related work

RPCs have been used for high-resolution optical satellite imaging since the launch of Ikonos in 1999 [START_REF] Grodecki | IKONOS stereo feature extraction -RPC approach[END_REF][START_REF] Tao | A comprehensive study of the rational function model for photogrammetric processing[END_REF][START_REF] Fraser | Sensor orientation via RPCs[END_REF][START_REF] Long | RPC estimation via 1 -norm regularized least squares (L1LS)[END_REF]. In the last decade, they have been proven to be extremely accurate for SAR acquisition systems as well [START_REF] Zhang | Evaluation of the RPC model for spaceborne SAR imagery[END_REF][START_REF] Zhang | Rational function modeling for spaceborne SAR datasets[END_REF].

The RPC model of a satellite image can be constructed using a set of correspondences between image and object space coordinates. Depending on the nature of these correspondences, the literature can be classified into terrain-dependent or independent methods (or a combination of both).

Terraindependent strategies use Ground Control Points (GCPs), whose object and image coordinates are known in advance relying on manual labeling or onsite measurements. Oppositely, terrain-independent methods derive virtual sets of 2D-3D point correspondences from other geolocation functions, usually a physical sensor model. Once the point correspondences are available, least squares algorithms are typically used to estimate the RPC coefficients that min-imize the error between the projected 3D points and their image locations.

Several works have underlined the importance of using uniformly distributed points in sufficient amount, covering the different parts of the image and the whole altitude range of the scene [START_REF] Tao | A comprehensive study of the rational function model for photogrammetric processing[END_REF][START_REF] Long | RPC estimation via 1 -norm regularized least squares (L1LS)[END_REF]. As a result, regularized least squares methods have become widely used to gain robustness to different configurations and enforce well-conditioned normal equations [START_REF] Tao | A comprehensive study of the rational function model for photogrammetric processing[END_REF][START_REF] Long | RPC estimation via 1 -norm regularized least squares (L1LS)[END_REF][START_REF] Zhang | Rational function modeling for spaceborne SAR datasets[END_REF][START_REF] Wang | RPC estimation via feature points for urban areas[END_REF]. Additionally, terrain-dependent strategies have explored the selection of optimal and balanced subsets of GCPs, e.g. [START_REF] Tao | A comprehensive study of the rational function model for photogrammetric processing[END_REF] propose a bucketing strategy or [START_REF] Wang | RPC estimation via feature points for urban areas[END_REF] study the benefits of encouraging correspondences located at building edges in urban scenarios. In contrast, terrain-independent methods can arbitrarily generate regular sets of points, but require special care to the boundaries and density of the structure, e.g. [START_REF] Zhang | Rational function modeling for spaceborne SAR datasets[END_REF] investigate the impact of different number of elevation layers for flat and mountainous areas.

Method

We follow a terrain-independent approach similar to [START_REF] Tao | A comprehensive study of the rational function model for photogrammetric processing[END_REF] to fit an RPC model to another input geolocation model. The data used to fit the model consists of a 3D grid of uniformly distributed Control Points (CNPs) within some longitude, latitude and altitude boundaries (Fig. 1). The 2D image point of each CNP can be obtained by projecting it with the input geolocation model, so that each sample results in 5 normalized values, i.e. (X i , Y i , Z i , r i , c i ). For simplicity, we drop the subscript n of normalized coordinates and replace it by i to refer to the sample index.

Using N CNPs, Equation 1 can be rewritten as a system of equations, in matrix form, following the derivation of [START_REF] Tao | A comprehensive study of the rational function model for photogrammetric processing[END_REF]:

W T I -W G = 0, (3)
where In Equation 3, W is a weight matrix with shape 2N × 2N , b(X i ) denotes the RPC polynomial b evaluated with the 3D coordinates of the i-th CNP; T is the design matrix with shape 2N × 78; I is the solution vector with the 78 RPC coefficients necessary to determine a, b, e and f in Equation 1; and G is a vector of length 2N containing the CNPs image coordinates. Equation 3 can be solved by least squares minimization to estimate I, using the normal equation

W = diag 1 b(X 1 ) , ..., 1 b(X N ) , 1 f (X 1 ) , ..., 1 f (X N ) b(X i ) = b(X i , Y i , Z i ), T = block diag [M r , M c ] M r , M c ∈ R N ×39 M ri = 1, Z i , Y i , ..., X 3 i , -r i Z i , -r i Y i , ..., -r i X 3 i M ci = 1, Z i , Y i , ..., X 3 i , -c i Z i , -c i Y i , ..., -c i X 3 i , I = [a 0 , ...,
T T W 2 T I -T T W 2 G = 0. (4)
To increase numerical stability, ridge estimation regularization [START_REF] Zhang | Rational function modeling for spaceborne SAR datasets[END_REF][START_REF] Wang | RPC estimation via feature points for urban areas[END_REF] is often added so that the normal equation becomes

(T T W 2 T + h 2 E)I -T T W 2 G = 0, ( 5
)
where E is the identity matrix and h is a scalar controlling the regularization that is applied. To choose the best regularization factor h, the L-curve criterion was introduced in [START_REF] Zhang | Rational function modeling for spaceborne SAR datasets[END_REF]. This heuristic computes the log norm of the solution (log I h ) versus the log norm of the residual (log W T I -W G h ) across different values of h that extend from the minimal to the maximal singular value of T. This curve usually has a L-shape, in which the optimum corresponds to the maximum regularization parameter that achieves a small residual. The value corresponding to the corner of the curve, at the position of maximal curvature, is taken to set h automatically. For non-weighted regularized least squares (i.e. weights are set to the identity, W = E), the L-curve criterion is fast since the curvature can be computed with closed form expressions [START_REF] Hansen | The use of the L-curve in the regularization of discrete ill-posed problems[END_REF]. Therefore, we first set W (0) = E (where the superscript denotes the iteration number) and use the Lcurve criterion to determine the optimal h and an initial solution I (0) . Then, for i ≥ 1 , W (i) is determined from I (i-1) and is plugged in Equation 5 to solve for I (i) iteratively (the SVD least squares solver is used for stability). The iterations stop when the change in terms of RMSE between the RPC projected CNPs and their image coordinates becomes lower than a tolerance value. After convergence some final ICCV (Iteration by Correcting Characteristic Value [START_REF] Zhang | Rational function modeling for spaceborne SAR datasets[END_REF]) iterations are computed to remove possible biases introduced by the regularization. The same stopping criterion based on the RMSE improvement is used. Each ICCV iteration k can be expressed as

(T T (W (k) ) 2 T + E)I (k) = T T (W (k) ) 2 G + I (k-1) . (6)

Experiments

Data and use cases description

• SAR. 32 Sentinel-1 SAR images (Table 1). The data is in Interferometric Wide Swath (IW) mode, each product contains 3 subswaths, and each subswath contains multiple bursts that need to be stitched together to get a continuous image. We construct the Range-Doppler physical sensor model and use the method from Section 3 to fit an equivalent RPC model for each image of the dataset.

• Optical. 47 WorldView-3 panchromatic images (Table 1), from the 2016 IARPA Multi-View Stereo 3D Mapping Challenge [START_REF] Bosch | A multiple view stereo benchmark for satellite imagery[END_REF]. The original RPCs of the images exhibit small inaccuracies, mainly due to inexact knowledge of the satellite attitude angles, which cause a 3D point to be projected to non corresponding pixels across different images. Bundle adjustment (BA) algorithms are a well-known approach to correct RPC errors [START_REF] Fraser | Sensor orientation via RPCs[END_REF][START_REF] Beyer | The Ames Stereo Pipeline: NASA's open source software for deriving and processing terrain data[END_REF][START_REF] Marí | To bundle adjust or not: A comparison of relative geolocation correction strategies for satellite multi-view stereo[END_REF]. We apply a BA similar to [START_REF] Beyer | The Ames Stereo Pipeline: NASA's open source software for deriving and processing terrain data[END_REF] to correc the projection function P of each RPC into a new P BA , expressed as

P BA (X) = P(R(X -T -C) + C), (7)
where X is a 3D point. That is, each RPC is corrected by applying a translation T followed by a rotation R around an approximate camera center C, before applying the original projection P. C is derived by regressing a projective model from each RPC model. Our method from Section 3 is used to fit P BA from the composition of P with T, R, C.

Performance assessment

To apply the method described in Section 3, the altitude limits of the grid of CNPs are set using the altitude ranges in Table 1. To assess the RPC fitting, we use a grid of Check Points (CKPs), which are located in the middle of each pair of consecutive CNPs. The RMSE between the image coordinates obtained by projecting the CKPs with the output RPC and the image coordinates obtained using the input geolocation model, measured in pixels, is used as evaluation metric. We set a tolerance of 10 -10 for the stopping criterion based on the RMSE improvement, and a maximum of 20 iterations for the weighted least squares and the ICCV iterations. Two type of experiments were conducted to assess the performance and robustness of the method:

-Varying surface area. For each image, we fit different RPCs by gradually increasing the longitude and latitude limits of the grid of CNPs from a small square centered at the image center to a larger square including the entire image. The number of CNPs is fixed, with 50 × 50 samples in the longitude and latitude dimensions and 10 elevation layers (25000 points in total). -Varying grid length. For each image, we fit different RPCs, by increasing the number of CNPs in each elevation layer, i.e. n × n where n is the grid length. The longitude and latitude boundaries of the grid are fixed using the equivalent limits of the image plane.

Overall, the results presented in Fig. 2 show that the RPCs constructed with our method approximate the geolocation models with very high accuracy for the two datasets and the two scenarios outlined in Section 4.1. The RMSE is in the order of 10 -4 pixels or less in both dimensions of the image plane for the majority of configurations that were tested, which emphasizes the robustness of the method. The experiments with different grid lengths (Fig. 2b and2d) show that 10 samples in the longitude and latitude directions is already a good choice and increasing this number beyond 20 does not result in significant improvements. The experiments with varying surface area (Fig. 2a and2c) show that both the overall RMSE values and its variation across the different images increase with the size of the area being fitted. This is probably due to the fact that the SAR physical sensor model is less smooth for large neighborhoods. A similar behavior is obtained with the optical dataset, where the original RPCs are known to behave locally as an affine camera [START_REF] Fraser | Sensor orientation via RPCs[END_REF][START_REF] Marí | To bundle adjust or not: A comparison of relative geolocation correction strategies for satellite multi-view stereo[END_REF].

Conclusion

This article described an automatic algorithm to fit the RPC model of a satellite image in a terrainindependent manner. The inputs of the method are a regular grid of 3D points (CNPs), with multiple elevation layers, and the 2D locations of the points on the image plane. We evaluated the method on real scenarios using collections of SAR and optical satellite images, and assessed its performance by varying the CNPs configuration. Finally, we release an opensource implementation of the algorithm as an easyto-use Python package.

Figure 1 :

 1 Figure 1: The RPC model is derived using a grid of 3D CNPs (Control Points) and its projection onto the satellite image.

  a 19 , b 1 , ..., b 19 , e 0 , ..., e 19 , f 1 , ..., f 19 ] T , G = [r 0 , ..., r N , c 0 , ..., c N ]T .

Figure 2 :

 2 Figure 2: RPC fitting error varying grid length and surface area on the SAR (a-b) and optical data (c-d) described in Section 4.1. Each vertical bar corresponds to the [-σ/2, µ, σ/2] values of the root-mean-square error (RMSE) evaluated across the different images of the datasets, in each dimension of the image plane, where µ corresponds to the mean and σ to the standard deviation.

Table 1 :

 1 SAR and optical data used in the experiments. The altitude range of the area covered by each collection of images is defined using the [min, max] values from the corresponding SRTM digital elevation model ±500 m to consider tall buildings or fine irregularities beyond bare ground level.

		SAR dataset Optical dataset
	platform	S1 A/B	WorldView-3
	number of images	32	47
	geographic area	Albania	Argentina
	(lon, lat) center	(18.82,41.02)	(-58.61,-34.47)
	first acquisition date	2019-08-03	2014-11-15
	last acquisition date	2020-02-05	2016-01-13
	altitude range (m)	[-533, 2969]	[-513, 548]

Also referred to as Rational Polynomial Coefficients camera model or Rational Function Model in the literature.

Acknowledgements Work partly financed by IDEX Paris-Saclay IDI 2016, ANR-11-IDEX-0003-02, Office of Naval research grant N00014-17-1-2552 and N00014-20-S-B001, DGA Astrid project « filmer la Terre » n o ANR-17-ASTR-0013-01, MENRT, and by a grant from Région Île-de-France.