Frédéric Sur
email: frederic.sur@loria.fr

Michel Grédiac
email: michel.grediac@univ-bpclermont.fr

F Sur Is With Loria

† M Grédiac

Sensor Noise Modeling by Stacking Pseudo-Periodic Grid Images Affected by Vibrations

Keywords: Image noise modeling, generalized Anscombe transform, grid method

This letter addresses the problem of noise estimation in raw images from digital sensors. Assuming that a series of images of a static scene are available, a possibility is to characterize the noise at a given pixel by considering the random fluctuations of the gray level across the images. However, mechanical vibrations, even tiny ones, affect the experimental setup, making this approach ineffective. The contribution of this letter is twofold. It is shown that noise estimation in the presence of vibrations is actually biased. Focusing on images of a pseudo-periodic grid, an algorithm to discard their effect is also given. An application to the generalized Anscombe transform is discussed.

Introduction

Noise is inherent to any digital signal acquisition device. Concerning images, noise is often simply modeled as an independent, identically distributed (i.i.d.) Gaussian process because many processing tasks and several noise sources are involved. However, a more accurate model is available for the raw images from a linear camera (see, e.g., [START_REF]Standard 1288, standard for characterization of image sensors and cameras[END_REF][START_REF] Healey | Radiometric CCD camera calibration and noise estimation[END_REF]). In this case, the camera response is linear which means that, disregarding the analog-to-digital conversion, the gray level measured at a pixel is proportional to the quantity of light photons arriving at a CCD sensor cell, plus a so-called dark current and read noise. A simplified yet realistic Poisson-Gaussian model is often used [START_REF]Standard 1288, standard for characterization of image sensors and cameras[END_REF][START_REF] Boulanger | Patch-based nonlocal functional for denoising fluorescence microscopy image sequences[END_REF][START_REF] Faraji | CCD noise removal in digital images[END_REF][START_REF] Foi | Practical Poissonian-Gaussian noise modeling and fitting for single-image raw-data[END_REF]:

u(x, y) = g • η p(x,y) (x, y) + δ(x, y) (1)
where u is the observed image, g > 0 is the gain of the electronic system, the number of collected photo-electrons η p(x,y) (x, y) at a pixel (x, y) is a random variable following a Poisson distribution of parameter p(x, y) (spatially independently distributed), and the dark current / read noise term δ(x, y) is modeled as a Gaussian white noise of mean µ δ and variance σ 2 δ . The random variables η p(x,y) (x, y) and δ(x, y) are assumed to be independent. Within this model, a straightforward calculation gives:

E(u(x, y)) = gp(x, y) + µ δ (2)
Var(u(x, y)) = g 2 p(x, y)

+ σ 2 δ (3)
where E and Var denotes respectively the expectation and the variance of any random variable. With (2) and (3), Var(u(x, y)) = gE(u(x, y)) + σ 2 δ -gµ δ (4)

The variance of the noise linearly depends on the expected gray level value. Provided that estimations of E(u(x, y)) and Var(u(x, y)) are available, a linear regression on the scatter plot of the pairs (E(u(x, y)),Var(u(x, y))) gives a reliable estimation of the noise parameters [START_REF] Boulanger | Patch-based nonlocal functional for denoising fluorescence microscopy image sequences[END_REF][START_REF] Boulanger | An adaptive statistical method for 4d-fluorescence image sequence denoising with spatio-temporal discontinuities preserving[END_REF][START_REF] Delpretti | Multiframe sure-let denoising of timelapse fluorescence microscopy images[END_REF][START_REF] Ramani | Deconvolution of 3d fluorescence micrographs with automatic risk minimization[END_REF]. This letter discusses the estimation of E(u) and Var(u) where u is a pseudo-periodic grid image. Characterizing noise in such images is of increasing interest in experimental mechanics when strain components are to be determined (see Section 3). These quantities being tiny in general, noise potentially strongly affects the results. Characterizing this effect is of prime importance to assess the actual metrological performance of contactless displacement and strain measurement techniques. Although it would be possible to characterize noise with a dedicated setting [START_REF]Standard 1288, standard for characterization of image sensors and cameras[END_REF], we propose here a low-cost measurement with a series of images of the static grid.

Image stacking is indeed certainly the most direct route to estimate the expectation and the variance at a given pixel, by computing respectively the empirical mean and variance.

However, even though the camera and the imaged object are carefully installed and fixed, some mechanical vibrations of very small amplitude are unavoidable in practice. The resulting peak-to-peak amplitude can be up to 4 µm measured on the sensor, while the pixel size of, e.g., a PCO Sensicam QE camera is 6.45 µm [START_REF] Grédiac | Effect of sensor noise on the resolution and spatial resolution of displacement and strain maps estimated with the grid method[END_REF]. It is shown in Section 2 that these vibrations actually give a biased estimator of the variance where the image gradient is not negligible. The special case of a pseudo-periodic grid, described in Section 3, is addressed in this letter. Taking advantage of periodicity, an algorithm is proposed in Section 4 to discard the small oscillations due to vibrations in the temporal evolution of the gray level at a given pixel. Experiments in Section 5 prove that this is an effective procedure to estimate the noise parameters, which can be used afterward in the generalized Anscombe transform.

Several studies are close to the topic of this letter: in [START_REF] Boulanger | Patch-based nonlocal functional for denoising fluorescence microscopy image sequences[END_REF] and [START_REF] Foi | Practical Poissonian-Gaussian noise modeling and fitting for single-image raw-data[END_REF], noise is estimated on a single image, which permits to discard the vibration problem. However a complex segmentation step is required. Noise is also characterized by time-averaging in [START_REF] Healey | Radiometric CCD camera calibration and noise estimation[END_REF] or [START_REF] Foi | Noise measurement for raw-data of digital imaging sensors by automatic segmentation of non-uniform targets[END_REF], but camera shake is overlooked.

Modeling noisy images affected by vibrations

An ideal, noise-and vibration-free image s(x, y) being considered, the gray value u(x, y, t) at pixel (x, y) of the t-th image of a series is modeled, following (1), by

u(x, y, t) = s(x + α t , y + β t ) + n s(x+αt,y+βt) (x, y, t) (5)
where

• (α t , β t ) is the translation vector modeling the displacement due to vibrations between the reference image s and the t-th image. It is modeled as an i.i.d. 0-mean process;

• n s (x, y, t) denotes 0-mean random variables of variance as + b, spatially and temporally independently distributed, with a = g and b = σ 2 δ -gµ δ from (4).

We also assume that the amplitude of the vibrations is small enough so that it is justified to identify s(x + α t , y + β t ) and its first-order Taylor expansion:

s(x + α t , y + β t ) = s(x, y) + (α t , β t )∇s(x, y) (6)
where ∇s is the gradient of s.

Let X be the empirical mean 1/T T t=1 X(t) of any random process (X(t)) 1≤t≤T . By definition of the random variables α t , β t and n s , and from ( 6), E(u(x, y, t)) = s(x, y) and E(u(x, y)) = 1/T E(u(x, y, t)) = s(x, y). Moreover, with ( 5), ( 6) and the law of the total variance1 :

Var(u(x, y, t)) = E(Var(u(x, y, t)|α t , β t )) + Var(E(u(x, y, t)|α t , β t )) (7)
= E(as(x + α t , y

+ β t ) + b) + Var(s(x + α t , y + β t )) (8)
= as(x, y) + b

+ Var(α t ) ∂s ∂x 2 + Var(β t ) ∂s ∂y 2 + 2 Covar(α t , β t ) ∂s ∂x ∂s ∂y (9)
While the empirical mean is (to the first order) an unbiased estimator of E(s), the empirical (sample) variance of u(x, y, t) is a biased estimator of as(x, y) + b. The bias simply writes ∇s T Var(α, β)∇s [START_REF] Foi | Noise measurement for raw-data of digital imaging sensors by automatic segmentation of non-uniform targets[END_REF] where Var(α, β) is the variance-covariance matrix of the vibration random process. The intuition behind this result is that a vibration orthogonal to the gradient of the image does not disturb the local estimation. In the remainder of this letter, we focus on pseudo-periodic grid images, inspired by a problem from experimental solid mechanics.

The grid method in experimental mechanics

The grid method [START_REF] Badulescu | Investigation of the grid method for accurate in-plane strain measurement[END_REF] is one of the full-field techniques for measuring in-plane displacement and strain of a specimen subjected to a load, and consequently lightly deformed. It consists in first depositing a regular grid on the surface of a material and then taking high-resolution pictures of the grid before and after deformation. Figure 1 (left) shows a close-up of a grid deposited on a deformed specimen. Displacement and strain maps are derived from the local deformation of the grid. A grid image is modeled as a 2D pseudo-periodic function:

s(x, y) = A 2 2 + γℓ(2πf x + φ 1 (x, y)) + γℓ(2πf y + φ 2 (x, y)) (11)
where A > 0 is the average field illumination, γ ∈ [0, 1] is the contrast of the oscillatory pattern, the line profile ℓ is a 2π-periodic real function with peak-to-peak amplitude equal to 1 and average value 0, f is the frequency of the carrier (it depends on the experimental setup, e.g., f = 1/5 pixel -1 in Figure 1), and φ 1 (x, y) and φ 2 (x, y) are the maps of the carrier phase modulations along the x-and y-axes respectively. Note that the model does not depend on a specific ℓ. In practice, the derivatives of φ 1 and φ 2 are very small with respect to 2πf for most constitutive materials used in real structures (10 -2 -10 -4 vs. 2π/5 pixel -1 ). This is an important property which permits to go beyond standard approaches based on Fourier transform used in, e.g., fringe pattern analysis in interferometry [START_REF] Surrel | Photomechanics, ser. Topics in Applied Physics[END_REF]. See [START_REF] Sur | Towards deconvolution to enhance the grid method for inplane strain measurement[END_REF][START_REF] Grédiac | Using deconvolution to improve the metrological performance of the grid method[END_REF], and [START_REF] Sur | Enhancing with deconvolution the metrological performance of the grid method for in-plane strain measurement[END_REF] for a digest.

Removing the non-random part of gray level fluctuations in grid images

The aim of this section is to design an algorithm to remove the bias in the estimation of the empirical variance, taking benefit of the pseudo-periodicity of the grid images. Noting p = 1/f the grid pattern pitch, we can see that ℓ is p-periodic. However, s(x, y) is not perfectly p-periodic because of the phases φ 1 and φ 2 which are not p-periodic. In this section, it is shown that s(x, y) can still be locally considered as a periodic function, because φ 1 and φ 2 have tiny variations as noted above.

For every integer k, l, (11) gives:

s(x + kp, y + lp) -s(x, y) = γA 2 (ℓ (2πf x + φ 1 (x + kp, y + lp)) -ℓ (2πf x + φ 1 (x, y))) + (ℓ (2πf y + φ 2 (x + kp, y + lp)) -ℓ (2πf y + φ 2 (x, y))) (12)
With Taylor's theorem and using the chain rule, there exists (ε

1 , ε 2 , ε ′ 1 , ε ′ 2 ) such that ε 1 , ε 2 ∈ [x, x + kp], ε ′ 1 , ε ′ 2 ∈ [y, y + lp],
and:

s(x + kp, y + lp) -s(x, y) = γA 2 (kp, lp), ℓ ′ (2πf x + φ 1 (ε 1 , ε ′ 1 ))∇φ 1 (ε 1 , ε ′ 1 ) + (kp, lp), ℓ ′ (2πf y + φ 2 (ε 2 , ε ′ 2 ))∇φ 2 (ε 2 , ε ′ 2 ) ( 13
)
where ℓ ′ is the derivative of the 1D periodic function ℓ, and < •, • > denotes the 2D dot product.

Let M be the supremum of |ℓ ′ | and S i be the supremum of the Euclidean norm ||∇φ i ||. Cauchy-Schwartz' inequality gives

|s(x + kp, y + lp) -s(x, y)| ≤ γAM (S 1 + S 2 )p 2 k 2 + l 2 . (14)
A is the mean value of s, thus the relative error made by approximating s(x + kp, y + lp) by s(x, y) is bounded by γM (S 1 +S 2 )p 2 √ k 2 + l 2 . The contrast γ is bounded by 1, and typical values are p = 5 pixels, M = 1 if ℓ is, e.g., a perfect sine, and sup ||∇φ • || ≃ 10 -3 . It is thus licit to approximate s(x + kp, y + lp) by s(x, y) as long as k, l are small with respect to the variations of φ. For example, with these numerical values, the bound on the relative error is less than 1% if √ k 2 + l 2 ≤ 4. This simplification gives with the model of (5):

u(x + kp, y + lp, t) = s(x + α t , y + β t ) + n s(x+αt,y+βt) (x + kp, y + lp, t) (15)
Hence for every x, y, t:

1 N k,l u(x + kp, y + lp, t) = s(x + α t , y + β t ) + 1 N k,l n s(x+αt,y+βt) (x + kp, y + lp, t) (16)
where (k, l) spans a set of cardinality N .

According to the hypotheses on n, 1 N k,l n s(x+αt,y+βt) (x + kp, y + lp, t) has 0-mean and a variance decreasing to 0 as 1/N . The noise term in ( 16) is thus negligible as soon as N is large enough. This justifies that 1 N k,l u(x + kp, y + lp, t) is a good estimation of s(x + α t , y + β t ) for any x, y, t, which corresponds to the grid images impaired by vibrations with noise ruled out.

Estimating noise parameters. Subtracting this estimation to u(x, y, t) for every x, y, t, gives, from ( 5), an estimation of n s(x+αt,y+βt) (x, y, t) for every x, y, t. Since the noise is assumed to be independently distributed, a straightforward calculation gives its expected temporal empirical variance:

1 T T t=1 Var(n s(x+αt,y+βt) ) = a 1 T T t=1 s(x + α t , y + β t ) + b (17)
By approximating s(x+α t , y+β t ) by its first-order Taylor expansion s(x, y)+α t ∂s ∂x (x, y)+ β t ∂s ∂y (x, y) in accordance with [START_REF] Boulanger | An adaptive statistical method for 4d-fluorescence image sequence denoising with spatio-temporal discontinuities preserving[END_REF], the expected variance of ( 17) becomes: as(x, y) + b + (α, β)∇s which converges to as(x, y) + b when the number of averaged images grows, since (α t , β t ) is a 0-mean i.i.d. process.

Summary.

The following process gives an unbiased estimation of the noise variance from a stack of T grid images.

1) Spatially average the gray values over pixels distant from multiple of p around a given pixel (x, y) for a given t, as in [START_REF]Sensicam QE -1288 datasheet[END_REF], to obtain the noise-free component s(x + α t , y + β t ), then subtract the result to the grid images in order to obtain the noise component. In practice, the average is computed with Gaussian weights with respect to the distance to the considered pixel (x, y). The periodicity assumption on the grid image is indeed valid only locally because of slight illumination variations and of manufacturing constraints. Moreover, localization is required by [START_REF] Grédiac | Using deconvolution to improve the metrological performance of the grid method[END_REF] to ensure the validity of Taylor's approximation. However, the standard deviation σ of this Gaussian function should be large enough (as N should be large to neglect the noise in ( 16)); here σ = 29 pixels turns out to give good results.

2) The temporal average of the s(x + α t , y + β t ) is a first-order approximation of s(x, y), and the temporal empirical variance is a first-order approximation of as(x, y) + b as soon as T is large enough to neglect α, β. Hence a and b can be estimated by a linear regression on the plot of temporal average vs. temporal empirical variance (calculated at each pixel).

We call Step 1 Non-Random Signal Reduction (NRSR).

Experiments

A series of T = 200 grid images from a PCO Sensicam QE camera is studied. It is impaired by tiny mechanical vibrations, cf. Figure 1 (right). It should be noted that the 16-bit output TIFF files are obtained by multiplying the 12-bit camera output by 2 4 . The gray levels (resp. variances) given here should thus be divided by 2 4 (resp. 2 8 ) to correspond to 12-bit data.

The procedure described in Section 4 is applied. is the noise-free component s(x + α t , y + β t ), estimated at the same locations as the ones in Figure 1 (right). Plot b) is the random part of the signal, centered at the temporal mean value. This is the result of the NRSR algorithm. The remaining signal is actually randomly distributed, and its amplitude is larger for the green line than for the red and blue lines. This is consistent with the linear dependence between the noise variance and the gray level. Plots c) and d) depict the distribution of the empirical variance as a function of the empirical mean on the raw data (as in Figure 1) and after NRSR (as in plot b). The former plot shows that the effect of tiny vibrations actually gives a strong . This gives a = 2 4 × 0.519 = 8.304 and b = (1.33 -45.46 × 0.519) × 2 8 = -5699.5 which is very close to our estimation. The software for [START_REF] Foi | Practical Poissonian-Gaussian noise modeling and fitting for single-image raw-data[END_REF][START_REF] Foi | Noise measurement for raw-data of digital imaging sensors by automatic segmentation of non-uniform targets[END_REF] has been tested but it does not give consistent results, probably because grid images are not adapted to the segmentation process.

Figure 4 shows the estimation of a and b against the number T of stacked images. T ≥ 100 gives reliable results.

As an application (motivated by [START_REF] Grédiac | Effect of sensor noise on the resolution and spatial resolution of displacement and strain maps estimated with the grid method[END_REF]), we check whether the estimated parameters permit to correctly apply the generalized Anscombe transform (GAT, see [START_REF] Murthagh | Image restoration with noise suppression using a multiresolution support[END_REF], after [START_REF] Anscombe | The transformation of Poisson binomial and negative-binomial data[END_REF]) which stabilizes to 1 the variance of any Poisson-Gaussian variable like the one in [START_REF]Standard 1288, standard for characterization of image sensors and cameras[END_REF]. With our notations, the GAT writes:

GAT(s(x, y)) = 2 g gs(x, y) + 3 8 g 2 + σ 2 δ -gµ δ (18)
Estimating the parameters of the GAT from the noise model was discussed in [START_REF] Boulanger | An adaptive statistical method for 4d-fluorescence image sequence denoising with spatio-temporal discontinuities preserving[END_REF] and also [START_REF] Boulanger | Patch-based nonlocal functional for denoising fluorescence microscopy image sequences[END_REF][START_REF] Delpretti | Multiframe sure-let denoising of timelapse fluorescence microscopy images[END_REF][START_REF] Ramani | Deconvolution of 3d fluorescence micrographs with automatic risk minimization[END_REF]. Figure 3 a) shows empirical variance against empirical mean after NRSR and GAT. Noise is now homoscedastic, with a variance roughly constant, equal to 1. The benefit of the NRSR is graphically assessed by comparing the spatial distribution of the empirical standard deviations. Map b) (directly estimated from the raw data, without NRSR and GAT) clearly shows a grid pattern, due to the relation [START_REF] Grédiac | Effect of sensor noise on the resolution and spatial resolution of displacement and strain maps estimated with the grid method[END_REF] between the variance and the intensity. If the GAT is applied on the raw data (with parameters estimated by linear regression on the raw data, i.e., the red line in Figure 2, plot c), then Map c) still shows a grid pattern. This means that NRSR is required to discard the effect of vibrations. Map d), estimated after GAT with the data processed by NRSR, has indeed no visible pattern. The noise has now a variance normalized to 1, which illustrates the efficiency of our estimation procedure for a and b.

Conclusion

Although observing the fluctuation of the gray level against time is a way to characterize sensor noise, it is demonstrated here that the empirical variance is biased by vibrations, which are likely to occur even in a controlled experimental setting. In the case of a pseudoperiodic grid, it is still possible to discard the bias, giving reliable estimates of the noise parameters which permits us to successfully stabilize the noise variance. Such an approach makes it possible in experimental mechanics to use easily low-cost cameras whose noise characteristics are not publicly available.

Figure 1 :

 1 Figure1: Left image: Close-up of a grid. In this typical example, the deformations of the grid are barely visible to the naked eye. Right image: Typical time evolution of the gray level of three distant pixels lying between lines of the grid. In-phase oscillations are caused by mechanical vibrations and cannot be explained by an independently distributed digital noise. The grid or the camera is not rigorously static.

Figure 2 :

 2 Figure 2: a) Estimated non-random oscillations; b) Noise component around the expected gray level; c) Empirical variance vs. empirical mean gray level on raw data; d) Empirical variance vs. empirical mean gray level after NRSR (same axis as in c).

Figure 2 Figure 3 :

 23 Figure 3: a) Empirical variance vs. empirical mean gray level after generalized Anscombe transform; b) Map of the standard deviation of the noise on raw data; c) Map after GAT but without NRSR on raw data; d) Map after applying first NRSR then GAT.

Figure 4 :

 4 Figure 4: Estimation of the regression parameters a (blue curve, left axis) and b (green curve, right axis) as a function of the number of stacked images.

If X and Y are random variables, Var(X) = E(Var(X|Y )) + Var(E(X|Y )).