Thanh Hai Thai

Florent Retraint

Rémi Cogranne
email: remi.cogranne@utt.fr

Generalized Signal-Dependent Noise Model and Parameter Estimation for Natural Images

Keywords: Signal-Dependent Noise Model, Noise Measurement, Noise Parameter Estimation, Denoising

The goal of this paper is to propose a generalized signal-dependent noise model that is more appropriate to describe a natural image acquired by a digital camera than the conventional Additive White Gaussian Noise model widely used in image processing. This non-linear noise model takes into account effects in the image acquisition pipeline of a digital camera. In this paper, an algorithm for estimation of noise model parameters from a single image is designed.

Then the proposed noise model is applied with the Local Linear Minimum Mean Square Error filter to design an efficient image denoising method.

Introduction

Noise has been studied for decades in computer vision, image processing and statistical signal processing because of its impact in various applications such as image denoising, image segmentation or edge detection. To improve performance in those applications, it is important to identify noise characteristics. Noise models proposed in the literature can be roughly divided into two groups: signal-independent and signal-dependent. A typical model for the group of signal-independent noise is the Additive White Gaussian Noise (AWGN) that is widely used in image processing.

However, this signal-independent AWGN model is not relevant due to the dominant contribution of the Poisson noise corrupting a natural image acquired by imaging device [START_REF] Healey | Radiometric CCD camera calibration and noise estimation[END_REF][START_REF] Foi | Practical Poissonian-Gaussian noise modeling and fitting for single-image raw-data[END_REF]. While signal-independent noise models assume the stationarity of noise in the whole natural image, regardless original pixel intensity, signal-dependent noise models take into account the proportional dependence of noise variance on the original pixel intensity. Signal-dependent noise models include Poisson noise or film-grain noise [START_REF] Kuan | Adaptive noise smoothing filter for images with signal-dependent noise[END_REF], Poisson-Gaussian noise [START_REF] Luisier | Image denoising in mixed Poisson-Gaussian noise[END_REF][START_REF] Makitalo | Optimal inversion of the generalized Anscombe transformation for Poisson-Gaussian noise[END_REF], heteroscedastic noise model [START_REF] Foi | Practical Poissonian-Gaussian noise modeling and fitting for single-image raw-data[END_REF][START_REF] Thai | Camera model identification based on the heteroscedastic noise model[END_REF], and non-linear noise model [START_REF] Faraji | CCD noise removal in digital images[END_REF][START_REF] Liu | Estimation of signal dependent noise parameters from a single image[END_REF]. The signal-dependent noise model gives the noise variance as a function of pixel's expectation. This function can be linear [START_REF] Luisier | Image denoising in mixed Poisson-Gaussian noise[END_REF][START_REF] Makitalo | Optimal inversion of the generalized Anscombe transformation for Poisson-Gaussian noise[END_REF][START_REF] Foi | Practical Poissonian-Gaussian noise modeling and fitting for single-image raw-data[END_REF][START_REF] Thai | Camera model identification based on the heteroscedastic noise model[END_REF] or non-linear [START_REF] Faraji | CCD noise removal in digital images[END_REF][START_REF] Liu | Estimation of signal dependent noise parameters from a single image[END_REF]. To identify noise characteristics Nikon D70: Estimated data Nikon D70: Fitted data Nikon D200: Estimated data Nikon D200: Fitted data Figure 1: Scatter-plot of pixels' expectation and variance from RAW images acquired by Nikon D70 and Nikon D200 cameras [START_REF] Gloe | The Dresden image database for benchmarking digital image forensics[END_REF].

or attenuate noise impact in many image processing applications, it is desirable to design an algorithm that estimates noise model parameters accurately.

Estimation of noise model parameters can be performed from a single image or multiple images. The contribution of this paper is threefold. Firstly, by modeling the main steps of the image acquisition pipeline, this paper starts from the heteroscedastic noise model [START_REF] Foi | Practical Poissonian-Gaussian noise modeling and fitting for single-image raw-data[END_REF][START_REF] Thai | Camera model identification based on the heteroscedastic noise model[END_REF] that accurately characterizes a natural RAW image and takes into account the impact of gamma correction to develop a generalized non-linear noise model. This model has not been proposed yet in the literature. Secondly, an algorithm for estimation of noise model parameters from a single image is proposed. Finally, the Local Linear Minimum Mean Square Error (LLMMSE) filter that was proposed in [START_REF] Kuan | Adaptive noise smoothing filter for images with signal-dependent noise[END_REF] is combined with the proposed noise model to design an efficient image denoising method.

The paper is organized as follows. Section 2 presents the generalized signal-dependent noise model for a natural image. Next, Section 3 designs the algorithm for estimation of noise model parameters. Section 4 presents numerical results on synthetic images and real natural images. Finally, Section 5 concludes the paper.

Generalized Signal-Dependent Noise Model

Image processing pipeline involves several steps from light capturing to image storage performed in a digital camera, see details about image processing pipeline in [START_REF] Ramanath | Color image processing pipeline[END_REF][START_REF] Thai | Statistical model of natural images[END_REF][START_REF] Thai | Statistical model of quantized DCT coefficients: Application in the steganalysis of Jsteg algorithm[END_REF]. After measuring light intensity at each pixel, RAW image that contains exactly information recorded by the image sensor goes through some typical post-acquisition processes, e.g. demosaicing, white-balancing and gamma correction, to render a full-color high-quality output image, referred to as TIFF image.

RAW image can be modeled by considering noise sources that corrupt the image during its acquisition process [START_REF] Healey | Radiometric CCD camera calibration and noise estimation[END_REF][START_REF] Foi | Practical Poissonian-Gaussian noise modeling and fitting for single-image raw-data[END_REF][START_REF] Thai | Camera model identification based on the heteroscedastic noise model[END_REF]. Typically, the RAW image model consists of a Poissonian part that addresses the photon shot noise and dark current and a Gaussian part for the remaining stationary disturbances, e.g. read-out noise. For the sake of simplification, the Gaussian approximation of the Poisson distribution can be exploited because of a large number of incident photons, which leads the heteroscedastic noise model [START_REF] Thai | Camera model identification based on the heteroscedastic noise model[END_REF][START_REF] Foi | Practical Poissonian-Gaussian noise modeling and fitting for single-image raw-data[END_REF] x i ∼ N µ x i , aµ x i + b ,

where x i denotes a RAW pixel, i ∈ {1, . . . , N} and N is the number of pixels. By convention, µ X and σ 2 X denote expectation and variance of a random variable X respectively. The index of color channel is omitted for simplicity. This model gives pixel's variance σ 2

x i as a linear function of pixel's expectation µ x i . As discussed in [START_REF] Thai | Camera model identification based on the heteroscedastic noise model[END_REF][START_REF] Foi | Practical Poissonian-Gaussian noise modeling and fitting for single-image raw-data[END_REF], the heteroscedastic noise model characterizes a RAW image more accurately than the conventional AWGN model widely used in image processing. The heteroscedastic noise model in a natural RAW image is illustrated in Figure 1.

In order to study noise statistics in a TIFF image, it is necessary to consider the effects of post-acquisition operations. In this paper, the demosaicing and white-balancing operations are assumed to be linear [START_REF] Ramanath | Color image processing pipeline[END_REF][START_REF] Thai | Statistical model of quantized DCT coefficients: Application in the steganalysis of Jsteg algorithm[END_REF]. Therefore, a short calculation shows that the white-balanced pixel still follows the Gaussian distribution and the relation between its expectation and variance remains linear

y i ∼ N µ y i , ãµ y i + b , (2)
where y i denotes the white-balanced pixel and the parameters (ã, b) differ from the parameters (a, b) due to effects of demosaicing and white-balancing. The white-balanced pixel y i can be equivalently rewritten as

y i = µ y i + η y i with η y i ∼ N 0, ãµ y i + b , (3)
where η y i represents the zero-mean signal-dependent noise after white-balancing. Meanwhile, the gamma correction is defined as the following element-wise power-law expression

z i = y 1 γ i = (µ y i + η y i ) 1 γ = µ 1 γ y i 1 + η y i µ y i 1 γ , ( 4
)
where γ is the correction factor (typically, γ = 2.2) and z i denotes the gamma-corrected pixel. The first order of Taylor's series expansion of (1 + x)

1
γ at x = 0 leads to

z i = µ 1 γ y i + 1 γ µ 1 γ -1 y i η y i + o η y i µ y i ≈ µ z i + 1 γ µ 1-γ z i η y i , (5)
where

µ z i = µ 1 γ
y i is the expected value of the gamma-corrected pixel z i . Taking expectation and variance on the both sides of the equation ( 5), it follows that

σ 2 z i = 1 γ 2 µ 2-2γ z i σ 2 y i = 1 γ 2 µ 2-2γ z i (ãµ γ z i + b). (6)
Finally, the gamma-corrected image undergoes the quantization Q ∆ with step ∆ in the image acquisition pipeline.

Under mild assumptions [START_REF] Widrow | Statistical theory of quantization[END_REF], the quantization noise can be modeled as an additive noise that is uniformly distributed and uncorrelated with the input signal. Taking into account the variance of the quantization noise, the generalized noise model of a natural image is derived as

σ 2 z i f (µ z i ; ã, b, γ) = 1 γ 2 µ 2-2γ z i (ãµ γ z i + b) + ∆ 2 12 , (7)
where, to simplify the notations, z i is referred to as the final output pixel. For the sake of simplification, it is assumed that the quantization step is unitary, i.e. ∆ = 1. Since this generalized noise model accounts for heteroscedasticity of noise, it is more appropriate to characterize than existing non-linear models used in [START_REF] Liu | Estimation of signal dependent noise parameters from a single image[END_REF][START_REF] Faraji | CCD noise removal in digital images[END_REF]. The generalized noise model ( 7) is illustrated in Figure 2.

Estimation of Noise Model Parameters from a Single Image

The generalized noise model ( 7) is non-linear, which causes a difficulty of estimating the noise model parameters.

When the gamma factor γ is known in advance, an obvious approach is to invert the gamma correction for obtaining again the heteroscedastic relation [START_REF] Kuan | Adaptive noise smoothing filter for images with signal-dependent noise[END_REF], and then to perform the Weighted Least Squares (WLS) estimation as proposed in [START_REF] Thai | Camera model identification based on the heteroscedastic noise model[END_REF]. Unfortunately, this approach leads to many problems in practice [START_REF] Kleihorst | An adaptive order-statistic noise filter for gamma-corrected image sequences[END_REF]. Firstly, the value of γ can not be known in practice. One method is proposed in [START_REF] Farid | Blind inverse gamma correction[END_REF] to estimate γ blindly without calibration information or knowledge of imaging device. However the stability of this method on a large real image database is still questioned. Secondly, even when the value of γ is exactly known, the effect of the quantization Q ∆ makes the inversion of the gamma correction illconditioned. Finally, this non-linear inversion would introduce artefacts into the signal, which prevents from obtaining a subsequent good estimation of parameters. Therefore, the goal of this section is to develop an algorithm that works directly on the non-linear generalized noise model [START_REF] Faraji | CCD noise removal in digital images[END_REF] for estimating the noise model parameters.

This section presents an algorithm for estimation of noise model parameters from a single image. The proposed algorithm consists of three fundamental steps: homogeneous block detection, level-set segmentation, and Maximum Likelihood (ML) estimation of parameters. The first two steps aim at detecting homogeneous blocks and partition the image into non-overlapping level sets (or segments) in which the pixels are assumed to be independent and identically distributed. Thus local expectations and local variances in each segment can be calculated, allowing to estimate the noise model parameters simultaneously.

Homogeneous Block Detection and Level-Set Segmentation

Let Z be a two-dimensional matrix representing a natural image. Firstly, an estimation of image structure is performed using a denoising filter D: Z app = D(Z) where Z app denotes the approximate image structure. The residual image Z res , which is the difference between the noisy image Z and the denoised image Z app , is further used to estimate local noise variances. Since it is desirable that the proposed algorithm can be further applied on JPEG images, and JPEG compression works separately on each 8 × 8 block, it is proposed to decompose the image Z (accordingly Z app and Z res ) into 64 vectors of pixels z L = (z L,1 , . . . , z L,N b ), where L ∈ {1, . . . , 64} denotes the location index in the 8 × 8 grid and N b is the number of blocks. Therefore, the vector z L contains all the pixels at the same location of the 8 × 8 grid and the pixels (z 1,B , . . . , z 64,B ) are in the same block B.

In order to identify if a 8 × 8 block is homogeneous or contains an edge or discontinuity, it is proposed to calculate the standard deviation of each block and compare it with a threshold τ. The median of absolute deviations (MAD), which is considered as a robust estimator of standard deviation [START_REF] Rousseeuw | Alternatives to the median absolute deviation[END_REF], is employed to calculate the standard deviation of the block. Therefore, the standard deviation of block B is calculated in the DCT domain as follows ŝB = 1.4826

Here, the denoised image Z app is employed instead of the noisy image Z because the noise may severely contaminate the calculation of standard deviation. Moreover, only 63 AC coefficients are used in [START_REF] Liu | Estimation of signal dependent noise parameters from a single image[END_REF]. The DC coefficient is excluded. The block B is selected if the standard deviation ŝB is smaller than the threshold τ. Hence the set of homogeneous blocks is defined by

S = 1 ≤ B ≤ N b : ŝB ≤ τ . (9)
After detecting homogeneous blocks, it is proposed to use only a sub-image z L for partitioning into K nonoverlapping segments by dividing the dynamic range. Each segment S k , k ∈ {1, . . . , K} is defined by

S k = z L,B : z app L,B ∈ u k - ∆ k 2 , u k + ∆ k 2 , B ∈ S . (10)
The number of segments K is set to the number of quantization levels, e.g. K = 2 8 and ∆ k = 1. If the value of K is larger than the number of quantization levels, the result is finer, i.e. the pixels in each level can be more probably identically distributed, but the number of pixels is smaller, which could lead to a case that there is not enough data for the subsequent parameter estimation. The situation is opposite when the value of K is smaller than the number of quantization levels. For the sake of clarity, the pixel in each segment S k is now denoted z k,i , i ∈ {1, . . . , n k } where n k is the number of pixels in segment S k . Analogously, z app k,i and z res k,i denote respectively its denoised value and residual value.

Maximum Likelihood Estimation

Consequently, local expectation and local variance in each segment are given by μk

= 1 n k n k i=1 z app k,i (11)
σ2 k = 1 n k -1 n k i=1 (z res k,i -z res k ) 2 with z res k = 1 n k n k i=1 z res k,i . (12)
Because the local expectation μk is calculated as the average of all denoised value in each segment, it is assumed that its variance is negligible when the number of pixels is large, i.e. the local expectation μk is close to the true value µ k : μk µ k . Meanwhile, the variance of σ2 k is more crucial and needs to be treated carefully. In virtue of Lindeberg Central Limit Theorem (CLT) [17, theorem 11.2.5], for a very large number of pixels n k , the local variance σ2 k follows the Gaussian distribution

σ2 k ∼ N(σ 2 k , d k σ 4 k ) with d k = 2 n k , ( 13
)
where σ 2 k = f (µ k ; ã, b, γ) is the true variance with respect to µ k . The Figure 2 illustrates the scatter-plot of couples { μk , σ2 k } K k=1 and the generalized noise model [START_REF] Faraji | CCD noise removal in digital images[END_REF] in natural images in JPEG format acquired by Nikon D70 and Nikon D200 cameras [START_REF] Gloe | The Dresden image database for benchmarking digital image forensics[END_REF].

The ML approach is used to fit the global parametric model σ 2 k = f (µ k ; ã, b, γ) to the scatter-plot of couples { μk , σ2

k } K k=1 . The log-likelihood function of K segments is given by

L = - 1 2 K k=1 log 2πd k f 2 ( μk ; ã, b, γ) + σ2 k -f ( μk ; ã, b, γ) d k f 2 ( μk ; ã, b, γ) . (14)
Here, because the true value µ k is unknown in practice, µ k is replaced by μk in the log-likelihood function L. The ML estimates of (ã, b, γ) are obtained by maximizing the log-likelihood function L

( â, b, γ) = arg max (ã, b,γ) L(ã, b, γ). ( 15
)
Because there is no closed form for ML estimates, the problem ( 15) is proposed to be solved numerically by using the Nelder-Mead method [START_REF] Nelder | A simplex method for function minimization[END_REF].

Numerical Results

It can be noted that the accuracy of the homogeneous block detection and segmentation depends on the performance of the denoising filter D. We have conducted some denoising methods such as Gaussian filter, Wiener filter, wavelet-based filter [START_REF] Mihçak | Spatially adaptive statistical modeling of wavelet image coefficients and its application to denoising[END_REF] and Block-Matching and 3D (BM3D) filter [START_REF] Dabov | Image denoising by sparse 3-d transform-domain collaborative filtering[END_REF]. The usual denoising methods as Gaussian filter and Wiener filter provide very poor results. The results provided by the wavelet-based filter and BM3D filter are equivalent but the BM3D filter takes much processing time. The wavelet-based denoising filter is employed in this paper because of its relative accuracy and computational efficiency.

The proposed algorithm for estimation of noise model parameters requires an appropriate threshold τ such that we have sufficient statistics for estimation process. In this paper, the threshold τ is defined as the median of absolute deviations of all residual pixels

z res i τ = 1.4826 • MAD z res 1 , . . . , z res N . (16)
This threshold is simple and efficient for rejecting blocks with strong edges. Besides, the sub-image used segmentation and parameter estimation corresponds to the location (4, 4) of the 8×8 grid since compression error is higher for pixels near block boundaries, and especially high at block corners [START_REF] Robertson | DCT quantization noise in compressed images[END_REF].

Parameter Estimation on Synthetic Images

The reference images from TID2008 database [START_REF] Ponomarenko | TID2008 -a database for evaluation of full-reference visual quality assessment metrics[END_REF] are chosen to evaluate the accuracy of the proposed approach.

The parameters (ã, b, γ) are used to generate synthetics images according to the generalized noise model [START_REF] Faraji | CCD noise removal in digital images[END_REF]. Those parameters are estimated from natural JPEG images that are acquired by Nikon D70 and Nikon D200 cameras (see Figure 2). The synthetic images are then compressed with different quality factors {70, 80, 90}. The Table 1 shows the statistics of estimated parameters on non-compressed TIFF images and JPEG images with different quality factors. It can be noted that the estimated parameters are close to the ground truth. Moreover, the distortion caused by JPEG compression with moderate-to-high quality factors (Q ≥ 70) seems to weekly interfere in the estimation process. Since the generalized noise model [START_REF] Faraji | CCD noise removal in digital images[END_REF] has not been proposed yet in the literature, there is no existing algorithm for noise parameter estimation to compare with the proposed one.

Parameter Estimation on Natural Images

To highlight the relevance of the proposed approach, experiments are then conducted on a large image database.

The Dresden database [START_REF] Gloe | The Dresden image database for benchmarking digital image forensics[END_REF] that covers different camera devices, different imaged scenes, different camera settings and different environmental conditions is chosen for this experiment. All images of the database are acquired with the highest available JPEG quality setting and maximum available resolution. Figure 3 shows estimated parameters (ã, b) over 1000 JPEG images of different camera models. As expected, the estimated parameters of the same camera model are close to each other. Furthermore, it is desirable to compare estimated gamma provided by the proposed algorithm with the algorithm proposed by Farid [START_REF] Farid | Blind inverse gamma correction[END_REF]. Figure 4 shows estimated gamma of the two algorithms on the JPEG images taken from Nikon D200 camera model. It can be noted that the variability of gamma estimated by the proposed algorithm is considerably smaller than Farid's, thus the amount of gamma can be estimated more efficiently. 2 shows the averaged Peak Signalto-Noise Ratio (PSNR) of the extended LLMMSE filter compared with other denoising methods, e.g. some classical methods such as Wiener filter and wavelet-based filter [START_REF] Mihçak | Spatially adaptive statistical modeling of wavelet image coefficients and its application to denoising[END_REF], and some state-of-the-art methods such as BM3D filter [START_REF] Dabov | Image denoising by sparse 3-d transform-domain collaborative filtering[END_REF] and Shape-Adaptive DCT (SA-DCT) filter [START_REF] Foi | Pointwise shape-adaptive DCT for high-quality denoising and deblocking of grayscale and color images[END_REF]. It can be noted that the extended LLMMSE filter outperforms the Wiener, wavelet-based, BM3D, and is rather equivalent to the SA-DCT filter. This can be justified due to the fact that the SA-DCT is designed for signal-dependent noise whilst the others do not consider the non-stationarity of noise. Different denoising filters are illustrated in Figure 5.

Remark 1. One of the main contributions of the paper is to propose the generalized noise model that has not been provided yet in the literature. This work is accomplished by studying the main steps of the image processing pipeline inside a digital camera. The relevance of the proposed model is highlighted by applying on a real image database.

However, this approach involves a limitation, namely that the demosaicing and white balancing are not completely modeled. The proposed approach assumes that those operations are linear, which simplifies the statistical study and results in a more exploitable model, say the generalized noise model. It must be noted that the demosaicing always involves a convolution operation, which requires us to consider multivariate distribution. The resulting model could be more relevant but inexploitable in practice. Overall, the generalized noise model can be seen as a trade-off between the reality in image acquisition and the exploitability in practice.

Conclusion

This paper proposes a novel generalized signal-dependent noise model that is more relevant to characterize a natural image acquired by a digital camera. An algorithm for estimating noise model parameters accurately from a single image is also designed. The proposed algorithm can work on JPEG images with moderate-to-high quality factors. Another strength of the proposed algorithm is the ability to estimate the gamma factor more efficiently.

Moreover, the LLMMSE filter is extended by combining with the proposed generalized noise model. The proposed noise model could be useful in many applications. A first step is to exploit the parameters of the generalized noise model as camera fingerprint for camera model identification, as proposed in [START_REF] Thai | Camera model identification based on generalized noise model in natural images[END_REF].

  From a practical point of view, this paper mainly focuses on noise parameter estimation from a single image. Several methods have been proposed in the literature for estimation of signal-dependent noise parameters[START_REF] Foi | Practical Poissonian-Gaussian noise modeling and fitting for single-image raw-data[END_REF][START_REF] Thai | Camera model identification based on the heteroscedastic noise model[END_REF][START_REF] Faraji | CCD noise removal in digital images[END_REF][START_REF] Liu | Estimation of signal dependent noise parameters from a single image[END_REF]. They rely on similar basic steps but differ in details. The common methodology starts from obtaining local estimates of noise variance and image content, then performing the curve fitting to the scatter-plot based on the prior knowledge of noise model. The existing methods involve two main difficulties: influence of image content and spatial correlation of noise in a natural image. In fact, homogeneous regions where local means and variances are estimated are obtained by performing edge detection and image segmentation. However, the accuracy of those local estimates may be contaminated due to the presence of outliers (textures, details and edges) in the homogeneous regions. Moreover, because of the spatial correlation, the local estimates of noise variance can be overestimated. Overall, the two difficulties may result in inaccurate estimation of noise parameters.

Figure 2 :

 2 Figure 2: Scatter-plot of pixels' expectation and variance from JPEG images acquired by Nikon D70 and Nikon D200 cameras [9].

Figure 3 :Figure 4 :

 34 Figure 3: Estimated parameters (ã, b) on JPEG images issued from different camera models in Dresden image database [9].

Figure 5 :

 5 Figure 5: Illustration of different denoising filters for ã = 1, b = 10, γ = 0.85.

Table 1 :

 1 Parameter estimation on synthetic images

		TIFF images	JPEG images (Q = 90)	JPEG images (Q = 80)	JPEG images (Q = 70)
		Avg	Std Dev	Avg	Std Dev	Avg	Std Dev	Avg	Std Dev
	ã = -0.0012	-0.0012	2.2e -4	-0.0012	2.5e -4	-0.0013	2.8e -4	-0.00095	4.1e -4
	b = 0.11	0.114	0.0097	0.117	0.0102	0.126	0.015	0.132	0.025
	γ = 0.8	0.807	0.053	0.811	0.057	0.79	0.058	0.814	0.11
	ã = -0.0025	-0.0024	3.5e -4	-0.0024	4.2e -4	-0.0026	5.4e -4	-0.0023	5.8e -4
	b = 0.20	0.196	0.0084	0.191	0.0098	0.215	0.0117	0.191	0.021
	γ = 0.85	0.845	0.049	0.853	0.058	0.845	0.061	0.842	0.086

  • MAD DCT z

	app 1,B , . . . , z app 64,B .

Table 2 :

 2 PSNR of different denoising filters

		Simple Wiener filter	Wavelet-based filter [19]	BM3D filter [20]	SA-DCT filter [23]	Extended LLMMSE filter
	ã = -0.0012
	b = 0.11	30.06	41.33	42.85	49.75	50.97
	γ = 0.8
	ã = -0.0025
	b = 0.20	30.05	40.37	42.56	47.73	48.61
	γ = 0.85
	ã = 1
	b = 10	24.88	26.13	26.32	27.08	29.14
	γ = 0.85

Application to Image Denoising

To highlight the usefulness of the generalized noise model [START_REF] Faraji | CCD noise removal in digital images[END_REF], it is proposed to combine it with the LLMMSE filter [START_REF] Kuan | Adaptive noise smoothing filter for images with signal-dependent noise[END_REF] to design an efficient image denoising method. The LLMMSE filter is based on the non-stationary mean, non-stationary variance image model. From [START_REF] Makitalo | Optimal inversion of the generalized Anscombe transformation for Poisson-Gaussian noise[END_REF], the pixel z can be decomposed as

where u denotes the original pixel and η u is the zero-mean signal-dependent noise, σ 2 η u = ãµ γ u + b. The index of pixel is omitted for the sake of clarity. The original pixel u involves non-stationary mean and non-stationary variance.

The non-stationary mean describes the gross structure of an image and the non-stationary variance characterizes edge information of the image [START_REF] Kuan | Adaptive noise smoothing filter for images with signal-dependent noise[END_REF]. In the decomposition [START_REF] Lehmann | Testing Statistical Hypotheses, 3rd Edition[END_REF], the quantization noise is assumed to be negligible. As explained in [START_REF] Kuan | Adaptive noise smoothing filter for images with signal-dependent noise[END_REF], the LLMMSE filter for any signal-dependent noise model is formulated as

where µ z and σ 2 z are respectively local mean and local variance of the pixel z, µ u and σ 2 u are respectively local mean and local variance of the original pixel u. Since the noise η u is zero-mean, it follows from (17) that µ z = µ u . The LLMMSE filter [START_REF] Nelder | A simplex method for function minimization[END_REF] is the weighted sum of original signal mean µ u and the noisy observation z where the weight is determined as the ratio of the original signal variance and noise variance. A simple technique to obtain local statistics µ z and σ 2 z is to calculate over a sliding window of size (2r + 1) × (2c + 1)

Therefore, it remains to calculate the variance σ 2 u . From [START_REF] Lehmann | Testing Statistical Hypotheses, 3rd Edition[END_REF], the variance σ 2 z can be given by

By using the Taylor series expansion of u 2-2γ around µ u , the expression of E u 2-2γ can be simplified as

Combining ( 21) and ( 22), the expression of the variance σ 2 u is derived as

The LLMSSE filter [START_REF] Nelder | A simplex method for function minimization[END_REF] follows immediately.


Thanh Hai Thai
email: thanh_hai.thai@utt.fr

Florent Retraint

Rémi

Rémi Cogranne

Ro

Camera Model Identification Based on DCT Coefficient Statistics

Keywords: Digital Forensics, Camera Model Identification, Hypothesis Testing, Natural Image Model, Discrete Cosine Transform

Camera model identification based on DCT coefficient statistics

Introduction

The evolution of digital imaging technology and information technologies in the past decades has raised a number of information security challenges. Digital images can be easily edited, altered or falsified because of a large availability of low-cost image editing tools. Consequently, the credibility and trustworthiness of digital images have been eroded. This is more crucial when falsified images that are utilized as evidence in a courtroom could mislead the judgement and lead to either imprisonment for the innocent or freedom for the guilty. Digital image forensics has been emerged in response to the increasing need to verify the trustworthiness of digital images, see [START_REF] Stamm | Information forensics: An overview of the first decade[END_REF] and references therein for detailed introductions.

State of the Art

One of the key problems of digital image forensics is identification of image origin, which aims to verify whether a given image was acquired by a specific camera or determine camera models/brands as well as types of imaging mechanism (e.g. scanners, cell-phone cameras, or computer graphics). Basically, when an image is captured by a camera, it is stored with the metadata headers in the memory storage device. The metadata, e.g. Exchangeable Image File (EXIF) and JPEG headers, contain all recording and compression history. Therefore, a simplest way to determine the image's source is to read out directly from the metadata. However, such metadata headers are not always available in practice if the image is resaved in a different format or recompressed. Another problem is that the metadata headers are not reliable as they can be easily removed or modified using low-cost editing tools. Therefore, it is desirable for law enforcement agencies to build up a set of reliable forensic tools for image origin identification.

In general, there are two approaches to address the problem of image origin identification. Active forensics involves generating extrinsically security measures such as digital watermarks [START_REF] Podilchuk | Digital watermarking: algorithms and applications[END_REF] and digital signatures [START_REF] Friedman | The trustworthy digital camera: restoring credibility to the photographic image[END_REF], referred to as extrinsic fingerprints, and adding to the image file. However, active forensics is of limited extent due to many strict constraints in its protocols. In order to solve the problem of image origin identification in its entirety, passive forensics has been quickly evolved. In contrast to active forensics, passive forensics does not impose any constraint, nor require any prior information including the original reference image. Forensic analysts have only the suspect image at their disposal and must explore useful information from that image to gather forensic evidence and trace the acquisition device. The common philosophy in passive approach is to rely on inherent intrinsic fingerprints that the digital camera leaves in a given image. Passive forensics can be widely applied to millions of images that circulate daily on communication networks. This paper mainly addresses the origin identification of images acquired by digital cameras based on passive approach. It is important to distinguish the problem of camera instance identification and the problem of camera model/brand identification. More specifically, fingerprints used for camera instance identification should capture individuality, especially cameras coming from the same brand and model. For camera model/brand identification, it is necessary to exploit fingerprints that are shared between cameras of the same model/brand but discriminative for different camera models/brands. Passive forensic methods proposed for those problems can be divided into two fundamental categories.

Technically, any method proposed for image origin identification must respond to following questions:

1. Which fingerprints are utilized for identification? 2. How to extract these fingerprints accurately from a given image? 3. Under which frameworks is the method designed to exploit the discriminability of fingerprints extracted from images captured by different sources 1 and to calculate the similarity of fingerprints extracted from images captured by the same source? Every stage from real-world scene acquisition to image storage can provide intrinsic fingerprints for forensic analysis. Therefore to design a camera fingerprint, it is necessary to study the image processing pipeline of a digital camera. Although the image processing pipeline is common for most cameras, each processing step is performed according to manufacturers' own design. Thus the information left by each processing step is useful to trace down to the device source.

Several fingerprints have been proposed in the literature. Sensor Pattern Noise (SPN), which is caused by imperfections during the manufacturing process and non-uniformity of photo-electronic conversion due to inhomogeneity of silicon wafers, is used in [START_REF] Lukas | Digital camera identification from sensor pattern noise[END_REF] for camera instance identification. Two main components of the SPN are the Fixed Pattern Noise (FPN) and the Photo-Response Non-Uniformity (PRNU) noise. The FPN is also used for camera instance identification in [START_REF] Kurosawa | CCD fingerprint method-identification of a video camera from videotaped images[END_REF]. However it can be compensated by subtracting a dark frame from the output image, thus it is not a robust fingerprint and no longer used in later works. The PRNU is directly exploited in [START_REF] Chen | Determining image origin and integrity using sensor noise[END_REF][START_REF] Fridrich | Digital image forensics[END_REF][START_REF] Li | Color-decoupled photo response non-uniformity for digital image forensics[END_REF]. Moreover, PRNU can be also used for camera model identification as proposed in [START_REF] Filler | Using sensor pattern noise for camera model identification[END_REF] based on the assumption that the fingerprint obtained from images in the TIFF or JPEG format contains traces of post-acquisition processes (e.g. demosaicing)

that carry information about the camera model. Other fingerprints include lens aberration [START_REF] Choi | Source camera identification using footprints from lens aberration[END_REF], Color Filter Array (CFA) pattern and interpolation algorithms [START_REF] Swaminathan | Nonintrusive component forensics of visual sensors using output images[END_REF][START_REF] Cao | Accurate detection of demosaicing regularity for digital image forensics[END_REF], and JPEG compression [START_REF] Choi | Source camera identification by JPEG compression statistics for image forensics[END_REF], which are proposed for camera model/brand identification.

In general, the image origin identification problem can be formulated into two frameworks: supervised classification [START_REF] Mitchell | Machine Learning[END_REF] and hypothesis testing [START_REF] Lehmann | Testing Statistical Hypotheses, 3rd Edition[END_REF]. Compared with hypothesis testing framework, supervised classification framework is utilized by most of existing methods in the literature to identify camera brands/models. Based on above fingerprint, a forensic feature set is designed and employed in a machine learning algorithm, e.g. Support Vector Machines (SVM) [START_REF] Vapnik | The Nature of Statistical Learning Theory[END_REF]. Supervised classification framework involves three main drawbacks. To achieve high accuracy, this framework requires an expensive training stage that comprises many images with different characteristics (e.g. image content or camera settings) from various sources for representing a real-world situation, which might be unrealistic in practice. Another drawback is the choice of an appropriate forensic feature, which importantly affects the detection performance of the classifier. Besides, the analytic establishment of statistical performances still remains an open problem in supervised classification framework [START_REF] Scott | Performance measures for Neyman-Pearson classification[END_REF].

Even though those methods perform efficiently, they have been designed with a very limited exploitation of hypothesis testing theory and statistical image models. Therefore, their performance can not be analytically established and is only evaluated based on a large image database. Moreover, in the operational context, it is crucial to warrant a prescribed false alarm probability. Our previous work [START_REF] Thai | Camera model identification based on the heteroscedastic noise model[END_REF] proposed to design a statistical test within hypothesis test-ing framework based on the heteroscedastic noise model for camera model identification from RAW images. Recently, this work has been extended in [START_REF] Thai | Camera model identification based on the generalized noise model in natural images[END_REF] for camera model identification from JPEG images based on the generalized noise model [START_REF] Thai | Generalized signal-dependent noise model and parameter estimation for natural images[END_REF]. The proposed tests can guarantee a false alarm probability while ensuring a high detection performance on a large database.

Main Contributions of the Paper

This paper aims to design a statistical test for camera model identification from JPEG images, which is based on the same methodology proposed in our previous works [START_REF] Thai | Camera model identification based on the heteroscedastic noise model[END_REF][START_REF] Thai | Camera model identification based on the generalized noise model in natural images[END_REF]. The approach is based on fingerprint extracted in the Discrete Cosine Transform (DCT) domain. The main motivation behind fingerprint extraction in the DCT domain is that the statistics of DCT coefficients change with different sensor noises combining with various in-camera processing algorithms. Relying on an accurate model of DCT coefficients allows us to capture this statistical difference in natural images taken by different camera models. The main contributions are the following:

• This paper is based on the state-of-the-art statistical model of DCT coefficients [START_REF] Thai | Statistical model of natural images[END_REF][START_REF] Thai | Steganalysis of jsteg algorithm based on a novel statistical model of quantized DCT coefficients[END_REF][START_REF] Thai | Statistical model of quantized DCT coefficients: Application in the steganalysis of Jsteg algorithm[END_REF] for fingerprint extraction. The parameters (c, d) that characterize the linear relation between two parameters α and β specified in the proposed model of DCT coefficients are exploited as camera fingerprint for camera model identification.

• Stating the camera model identification problem in hypothesis testing framework, the paper studies the optimal detector given by the Likelihood Ratio Test (LRT) in an ideal context where all model parameters are known. This optimal detector serves as an upper-bound of any statistical test for the camera model identification problem.

• In the practical context, the model parameters are unknown. The paper proposes two Generalized Likelihood Ratio Tests (GLRTs) to deal with the difficulty of unknown parameters. The statistical performance of the GLRTs is analytically established. Moreover, the proposed GLRTs allow the guaranteeing of a prescribed false-alarm rate and the setting of decision threshold independently of the image content, which is crucial in an operational context. Numerical experiments also show that the loss of power of GLRTs compared with the LRT is negligible.

Organization of the Paper

The paper is organized as follows. Section 2 presents the camera fingerprint that is further exploited for camera model identification. Section 3 designs an algorithm for camera fingerprint extraction in the DCT domain. Section 4 states the camera model identification problem in the framework of hypothesis testing theory and studies the LRT assuming that all model parameters are known in advance. Section 5 designs two GLRTs to address the difficulty of unknown parameters. Section 6 presents numerical results of two proposed GLRTs on simulated and real JPEG images. Finally, Section 7 concludes the paper.

Design of Camera Fingerprint

To design a camera fingerprint, it is vital to deeply understand image characteristics during various stages of image processing pipeline and study image statistics. Image processing pipeline involves several steps from light capturing to image storage performed in a digital camera, see details about image processing pipeline in [START_REF] Ramanath | Color image processing pipeline[END_REF][START_REF] Thai | Statistical model of natural images[END_REF][START_REF] Thai | Statistical model of quantized DCT coefficients: Application in the steganalysis of Jsteg algorithm[END_REF]. After measuring light intensity at each pixel, RAW image that contains exactly information recorded by the image sensor goes through some typical post-acquisition processes, e.g. demosaicing, white-balancing and gamma correction, to render a full-color high-quality output image, referred to as TIFF image. Image compression can be also performed for ease of storage and transmission.

The study of noise statistics in a natural image from RAW format to TIFF format has been performed in our previous work [START_REF] Thai | Generalized signal-dependent noise model and parameter estimation for natural images[END_REF]. Firstly, the approach starts from the heteroscedastic noise model proposed in [START_REF] Thai | Camera model identification based on the heteroscedastic noise model[END_REF], which is established by modeling various noise sources during RAW image acquisition [START_REF] Healey | Radiometric CCD camera calibration and noise estimation[END_REF]. The heteroscedastic noise model characterizes more accurately the RAW image than the Additive White Gaussian Noise (AWGN) model widely used in image processing since it takes into account the contribution of Poisson noise in the acquisition process. The heteroscedastic noise model is given as

x i ∼ N µ x i , aµ x i + b , (1)
where x i denotes a RAW pixel. The index of color channel is omitted for simplicity. By convention, µ X and σ 2 X denote respectively expectation and variance of a random variable X. The parameters (a, b) was proposed in our previous work [START_REF] Thai | Camera model identification based on the heteroscedastic noise model[END_REF] as fingerprint for camera model identification from RAW images. Then, assuming the operations of demosaicing and white-balancing are linear, that approach [START_REF] Thai | Generalized signal-dependent noise model and parameter estimation for natural images[END_REF] takes into account the non-linear effect of gamma correction to develop a generalized noise model, given as

σ 2 z i = 1 γ 2 µ 2-2γ z i (ãµ γ z i + b), (2)
where z i denotes an output pixel, γ is the correction factor, and (ã, b) differ from the parameters (a, b) due to the operations of demosaicing and white balancing. It is also shown in [START_REF] Thai | Generalized signal-dependent noise model and parameter estimation for natural images[END_REF] that the generalized noise model is relevant to characterize JPEG images with moderate-to-high quality factors (QF ≥ 70). More details of the generalized noise model are given in [START_REF] Thai | Generalized signal-dependent noise model and parameter estimation for natural images[END_REF]. Similarly, the parameters (ã, b, γ) are proposed in [START_REF] Thai | Camera model identification based on the generalized noise model in natural images[END_REF] as fingerprint to identify camera models from JPEG images.

The next step in image processing pipeline is JPEG compression that involves transforming the TIFF image into the DCT domain. To capture statistics of DCT coefficients accurately, it is necessary to study the model of DCT coefficients. Based on the assumption that the pixels are identically distributed within 8 × 8 block, our previous work [START_REF] Thai | Statistical model of natural images[END_REF][START_REF] Thai | Steganalysis of jsteg algorithm based on a novel statistical model of quantized DCT coefficients[END_REF][START_REF] Thai | Statistical model of quantized DCT coefficients: Application in the steganalysis of Jsteg algorithm[END_REF] has recently proposed a novel model of DCT coefficients, given by

f I (u) = 2 π |u| β 2 α-1 2 β α Γ(α) K α-1 2 |u| 2 β , (3)
(α 0,1 , β -1 0,1 ) (α 8,8 , β -1 8,8 )
Figure 1: Estimated parameters (α, β) at frequency (0, 1) and (8, 8) of uniform images generated using ã

= 0.1, b = 2, γ = 2.2.
where f X denotes the probability density function (pdf) with respect to a random variable X, α is a positive shape parameter, β is a positive scale parameter, Γ(•) denotes the gamma function and K ν (x) denotes the modified Bessel function [26, chap. 5.5]. The proposed model of DCT coefficients outperforms the Laplacian, Generalized Gaussian, and Generalized Gamma model, see more details in [START_REF] Thai | Statistical model of quantized DCT coefficients: Application in the steganalysis of Jsteg algorithm[END_REF]. The parameters (α, β) can be estimated following the Maximum Likelihood (ML) approach as proposed in [START_REF] Thai | Statistical model of quantized DCT coefficients: Application in the steganalysis of Jsteg algorithm[END_REF].

Since the parameters (ã, b, γ) also contain information about camera model, after transforming into DCT domain, this information is expanded over different frequencies. Therefore, it is proposed to establish the relation between the parameters (ã, b, γ) and (α, β) to capture such information in the DCT domain. For the sake of simplification, this relation is given by

β -1 = c α + d, (4)
where the parameters (c, d) depend on (ã, b, γ) (see more details in Appendix A). This suggests that the parameters (c, d) can be also used for camera model identification. It can be said that while the relations (1) and ( 2) characterize the non-stationarity of noise in the spatial domain, the relation (4) characterizes this property in the DCT domain.

The relation (4) can capture the difference of noise statistics in natural images taken by different camera models.

Moreover, the linearity in relation ( 4) can facilitate the estimation of the parameters (c, d). It should be noted that in an image whose each 8 × 8 block is uniform, the same parameters (α, β) and (c, d) are shared among DCT coefficients at different frequencies. The relation ( 4) on such images is illustrated in Figure 1.

Extraction of Camera Fingerprint

An important requirement when using the parameters (c, d) as camera fingerprint is that they should be invariant to image content. Furthermore, to guarantee the above mathematical framework, it is necessary to work on homogeneous blocks. These considerations are addressed in this section.

Because of heterogeneity and noise non-stationarity in a natural image, the energy tends to be more located in lower frequencies. Consequently, DCT coefficients at different frequencies do not share the same parameters (α, β)

and (c, d). Therefore, the estimation of parameters (α, β) and (c, d) should be performed on each frequency separately.

In this paper it is proposed to arrange DCT coefficients into 64 vectors of coefficients according to the zig-zag order.

Let I k = (I k,1 , . . . , I k,N ) with k ∈ {1, . . . , 64}, be the vector of length N that contains coefficients at the frequency k.

The coefficients (I 1,i , . . . , I 64,i ) are in the same block. Analogously, let denote the parameters (α k , β k ) and (c k , d k ) with respect to the AC coefficients I k .

Estimation of Parameters

(c k , d k )
The proposed algorithm for estimation of parameters (c k , d k ) consists of three fundamental steps: image denoising, homogeneous block detection, and Least-Squares (LS) estimation [START_REF] Rao | Linear models : Least Squares and Alternatives, 2nd Edition[END_REF]. Image denoising step aims to attenuate the impact of image content. The detection of homogeneous blocks is performed subsequently to provide appropriate sample data for parameter estimation. The LS approach is applied straightforwardly as the relation ( 4) is linear.

Let Z be a two-dimensional matrix representing a natural image. To remove the image content, a denoising filter D is employed so that the residual image W is given by

W = Z -D(Z). (5)
If Z is a color image, the denoising filter D is performed on each color component, then three residual components are combined into one residual image using the usual conversion from RGB to grayscale

W = 0.2989W r + 0.587W g + 0.114W b , (6)
where W r , W g , and W b are respectively residuals of red, green, and blue component. The residual image W is then transformed into the DCT domain

I = DCT(W), ( 7
)
where I is the image of DCT coefficients of the residual image W.

For homogeneous block detection, this paper proposes to calculate the standard deviation of each block and compare it with a threshold λ. The median of absolute deviations (MAD), which is considered as a robust estimator of standard deviation [START_REF] Rousseeuw | Alternatives to the median absolute deviation[END_REF], is utilized to calculate the standard deviation of each block. The standard deviation of block i is calculated in the DCT domain as

ŝi = 1.4826 • MAD I 2,i , . . . , I 64,i . (8)
The DC coefficient I 1,i is excluded in the calculation ŝi . The block i is selected if the standard deviation ŝi is smaller than the threshold λ. The number of selected homogeneous blocks is denoted as N b .

Suppose that L couples ( αk,l , βk,l ), l ∈ {1, . . . , L} are available, the LS estimates of the parameters (c k , d k ) are given by with

          ĉk dk           = (H T k H k ) -1 H T k υ k (9)
H k =                   αk,1 1 . . . . . . αk,L 1                   and υ k =                   β-1 k,1 . . . β-1 k,L                  
, where H T k and H -1 k denote respectively the transpose and inverse of the matrix H k . The LS estimates (ĉ k , dk ) are unbiased and asymptotically equivalent to ML estimates in large samples [START_REF] Rao | Linear models : Least Squares and Alternatives, 2nd Edition[END_REF].

As showed above, the LS approach requires several couples (α k , β k ) for estimation of parameters (c k , d k ). One can collect L images and estimate a couple (α k,l , β k,l ) from all homogeneous blocks of each image following the ML approach [START_REF] Thai | Statistical model of quantized DCT coefficients: Application in the steganalysis of Jsteg algorithm[END_REF]. However, from the practical point of view, it is necessary to estimate the parameters (c k , d k ) from a single image. This is accomplished by extracting randomly a subset of n b homogeneous blocks from N b blocks, then performing the ML estimation of parameters (α k,l , β k,l ) on this subset.

Property of Parameters

(c k , d k )
When the image content is removed perfectly, the parameters (c k , d k ) remain identical for images with different image contents. However, in practice, due to the fact that the perfect denoising filter D is difficult to obtain, the DCT coefficients at low frequencies are still affected by image content. Meanwhile, the coefficients at high frequencies contain mostly noises because of the energy compaction property of DCT operation. Thus they are more relevant to exploit for camera model identification. Figure 2 shows the the linear relation (4) at frequency [START_REF] Li | Color-decoupled photo response non-uniformity for digital image forensics[END_REF][START_REF] Li | Color-decoupled photo response non-uniformity for digital image forensics[END_REF] of natural JPEG images taken by Canon Ixus 70 and Nikon D200 camera models. It should be noted that each point (α, β) in Figure 2 corresponds to one image. Figure 2 Moreover, for camera model identification problem, it is necessary to verify the discriminability of parameters (c k , d k ) for different camera models. The Figure 3 shows the parameters (c k , d k ) estimated from JPEG images at frequency [START_REF] Li | Color-decoupled photo response non-uniformity for digital image forensics[END_REF][START_REF] Li | Color-decoupled photo response non-uniformity for digital image forensics[END_REF] for different camera models. This figure clearly shows their discriminability between different camera models. Therefore, the parameters (c k , d k ) are proposed to be exploited as camera fingerprint to identify camera models in this paper.

Optimal Detector for Camera Model Identification Problem

Hypothesis Testing Formulation

Let analyze two camera models S 0 and S 1 . Each camera model S j , j ∈ {0, 1}, is characterized by the parameters (c k, j , d k, j ), k ∈ {1, . . . , K}, where K is the number of usable frequencies for camera model identification. For obvious reasons, it is assumed that (c k,0 , d k,0 ) (c k,1 , d k,1 ). In a binary hypothesis testing, the inspected image Z is either acquired by camera model S 0 or camera model S 1 . The goal of the test is to decide between two hypotheses defined by ∀k ∈ {1, . . . , K}, ∀i ∈ {1, . . . ,

N b }            H 0 = I k,i ∼ P θ k,0 , β -1 k,0 = c k,0 α k + d k,0 H 1 = I k,i ∼ P θ k,1 , β -1 k,1 = c k,1 α k + d k,1 , (10)
where P θ k, j , θ k, j = (α k , c k, j , d k, j ), denotes the probability distribution of DCT coefficients I k,i under hypothesis H j . As previously explained, this paper focuses on designing a test that allows the guaranteeing of a prescribed false-alarm probability. Hence, let

K α 0 = δ : sup θ 0 P H 0 δ(Z) = H 1 ≤ α 0
be the class of tests whose the false alarm probability is upper-bounded by the prescribed rate α 0 . Here θ 0 = (θ 1,0 , . . . , θ K,0 ) is the vector containing all parameters, P H j E stands for the probability of event E under hypothesis H j , j ∈ {0, 1}, and the supremum over θ has to be understood as whatever model parameters might be. Among all the tests in the class K α 0 , it is aimed at finding a test δ which maximizes the power function, defined by the correct detection probability:

β δ = P H 1 δ(Z) = H 1 .
The problem [START_REF] Choi | Source camera identification using footprints from lens aberration[END_REF] highlights three fundamental difficulties of the camera model identification. First, even when all model parameters (α k , c k, j , d k, j ) are known, the most powerful test, namely the LRT, has never been studied in the literature. The second difficulty concerns unknown parameters α k in practice. Finally, the camera parameters (c k, j , d k, j ) are also unknown, thus the hypothesis H j becomes composite.

Suppose that the camera model S 0 is available, thus forensic analysts can have access to its characteristics, or its fingerprints, i.e. its camera parameters (c k,0 , d k,0 ) can be known. Therefore, they can make a decision by checking whether the image under investigation Z contains the fingerprint (c k,0 , d k,0 ). It is proposed to solve the problem [START_REF] Choi | Source camera identification using footprints from lens aberration[END_REF] when the alternative hypothesis H 1 is composite, i.e. the camera parameters (c k,1 , d k,1 ) are unknown. It can be noted that a test that maximizes the correct detection probability whatever (c k,1 , d k,1 ) might scarcely exist [START_REF] Lehmann | Testing Statistical Hypotheses, 3rd Edition[END_REF]. The main goal of this paper is to study the LRT and to design the GLRTs to address the second and third difficulties.

Likelihood Ratio Test for Two Simple Hypotheses

When all model parameters are known, in virtue of the Neyman-Pearson lemma [15, theorem 3.2.1], the most powerful test δ * solving the problem [START_REF] Choi | Source camera identification using footprints from lens aberration[END_REF] is the LRT given by the following decision rule

δ * (Z) =                      H 0 if Λ(Z) = K k=1 N b i=1 Λ(I k,i ) < τ * H 1 if Λ(Z) = K k=1 N b i=1 Λ(I k,i ) ≥ τ * (11)
where the decision threshold τ * is the solution of the equation

P H 0 Λ(Z) ≥ τ * = α 0 ( 12
)
to ensure that the LRT is in the class K α 0 and the LR Λ(I k,i ) is defined as

Λ(I k,i ) = log P θ k,1 I k,i P θ k,0 I k,i , (13)
assuming that the DCT coefficients are statistically independent. From (3), it can be noted that the expression of the LR Λ(I k,i ) is difficult to exploit for subsequent stages, e.g. the design of the GLRT and analytic establishment of its statistical performance. Therefore it is proposed to simplify the LR Λ(I k,i ) to facilitate the study in the manner that it does not cause any loss of optimality.

Using the Laplace's approximation [START_REF] Tierney | Accurate approximations for posterior moments and marginal densities[END_REF][START_REF] Butler | Laplace approximations for hypergeometric functions with matrix argument[END_REF] (see more details in Appendix B), the function f I (u) can be approximated as

f I (u) ≈ |u| α-1 (2β) α 2 Γ(α) exp         -|u| 2 β         . (14)
Consequently, the LR Λ(I k,i ) can be simplified as

Λ(I k,i ) = log |I k,i | α k -1 (2β k,1 ) α k 2 Γ(α k ) exp -|I k,i | 2 β k,1 |I k,i | α k -1 (2β k,0 ) α k 2 Γ(α k ) exp -|I k,i | 2 β k,0 = α k 2 log β -1 k,1 β -1 k,0 - √ 2|I k,i | β -1 k,1 -β -1 k,0 . (15)
It should be noted that other polynomial expansions for the modified Bessel function K ν (x) are provided in [START_REF] Ryzhik | Tables of Integrals, Series, and Products[END_REF], so a polynomial approximation of f I (u) can be derived. However, those approximations are not considered in this paper.

The main advantage of the Laplace's approximation ( 14) is to provide an approximation of the form of exponential family function, which allows us to simplify the expression of the LR Λ(I k,i ). The approximating function ( 14) is used only for simplification of the LR. The estimation of parameters (α k , β k ) is always based on the exact function [START_REF] Friedman | The trustworthy digital camera: restoring credibility to the photographic image[END_REF].

In order to analytically establish the statistical performance of the LRT, it is necessary to characterize the statistical distribution of the LR Λ(Z) under each hypothesis H j . To this end, it is proposed to rely on the Lindeberg Central Limit Theorem (CLT) [15, theorem 11.2.5] that requires to calculate the expectation and variance of Λ(I k,i ).

Proposition 1. Under hypothesis H j , the first two moments of the LR Λ(I k,i ) are given by

m k, j E H j Λ(I k,i ) = α k 2 log β -1 k,1 β -1 k,0 - 2 √ π β 1 2 k, j Γ(α k + 1 2 ) Γ(α k ) β -1 k,1 -β -1 k,0 (16)
v k, j Var H j Λ(I k,i ) = 2 β -1 k,1 -β -1 k,0 2 × α k β k, j - 2β k, j π Γ 2 (α k + 1 2 ) Γ 2 (α k ) . ( 17
)
where E H j [•] and Var H j [•] respectively denote the mathematical expectation and variance under hypothesis H j .

Proof. of Proposition 1 is given in Appendix C.

In virtue of Lindeberg CLT, the statistical distribution of the LR Λ(Z) under hypothesis H j is derived as

Λ(Z) d → N m j , v j , (18)
where the notation d → denotes the convergence in distribution and the expectation m j and variance v j are given by

m j = K k=1 N b i=1 E H j Λ(I k,i ) = K k=1 N b m k, j (19)
v j = K k=1 N b i=1 Var H j Λ(I k,i ) = K k=1 N b v k, j . (20)
Since a natural image is heterogeneous, it is proposed to normalize the LR Λ(Z) in order to set the decision threshold independently of the camera parameters. The normalized LR is defined by Λ (Z) = Λ(Z)-m 0 √ v 0 . Accordingly, the corresponding LRT δ is rewritten as follows

δ (Z) =            H 0 if Λ (Z) < τ H 1 if Λ (Z) ≥ τ (21)
where the decision threshold τ is the solution of the equation P H 0 Λ (Z) ≥ τ = α 0 . The decision threshold τ and the power β δ are given in following theorem.

Theorem 1. In an ideal context where all the model parameters (α k , c k, j , d k, j ) are exactly known, the decision threshold and the power function of the LRT δ are given by

τ = Φ -1 (1 -α 0 ) (22)
β δ = 1 -Φ m 0 -m 1 + τ √ v 0 √ v 1 , (23)
where Φ(•) and Φ -1 (•) denotes respectively the cumulative distribution function of the standard Gaussian random variable and its inverse.

The test power β δ serves as an upper-bound of any statistical test for the camera model identification problem.

The test δ allows to warrant a prescribed false alarm rate and maximizes the detection probability. Since its statistical performance is analytically established, it can provide an analytically predictable result for any false alarm probability α 0 .

Practical Context: Generalized Likelihood Ratio Test

The scenario studied in the LRT may not be realistic because the parameters (α k , c k,1 , d k,1 ) are unknown in practice.

This section designs two GLRTs to deal with unknown parameters. It is proposed to replace unknown parameters by their ML estimates in the LR Λ(I k,i ) (15).

Generalized Likelihood Ratio Test with Unknown Parameters α k

In this subsection it is assumed that the camera parameters (c k, j , d k, j ) are known and we only deal with unknown nuisance parameters α k . By replacing unknown parameter α k by its ML estimate αk in the LR Λ(I k,i ) (15) (see more details about ML estimation of parameters (α k , β k ) in [START_REF] Thai | Statistical model of quantized DCT coefficients: Application in the steganalysis of Jsteg algorithm[END_REF]), the Generalized Likelihood Ratio (GLR) Λ 1 (I k,i ) can be given by

Λ 1 (I k,i ) = αk 2 log c k,1 αk + d k,1 c k,0 αk + d k,0 - √ 2|I k,i | c k,1 αk + d k,1 -c k,0 αk + d k,0 . (24)
The ML estimate αk is asymptotically consistent [START_REF] Lehmann | Testing Statistical Hypotheses, 3rd Edition[END_REF], i.e. it asymptotically converges in probability to its true value: αk p → α k . Therefore, from the Slutsky's theorem [15, theorem 11.2.11], the statistical distribution of the GLR

Λ 1 (Z) = K k=1 N b i=1 Λ 1 (I k,i
) under each hypothesis H j can be approximated as

Λ 1 (Z) d → N m j , v j , (25)
where the expectation m j and variance v j are given in [START_REF] Thai | Camera model identification based on the generalized noise model in natural images[END_REF] and [START_REF] Thai | Generalized signal-dependent noise model and parameter estimation for natural images[END_REF], respectively.

Similarly, the normalized GLR Λ 1 (Z) is defined by

Λ 1 (Z) = Λ 1 (Z)-m 0 √ v 0
. However, the expectation m 0 and variance v 0 can not be defined in practice since the parameters α k are unknown. Therefore, this paper proposes to replace α k by αk in [START_REF] Thai | Camera model identification based on the generalized noise model in natural images[END_REF] and ( 20) to obtain the estimates of m 0 and v 0 , denoted m0 and v0 . The normalized GLR Λ 1 (Z) can be given

in practice as Λ 1 (Z) = Λ 1 (Z)-m0 √ v0
. Since the estimates m0 and v0 are consistent, it also follows from Slutsky's theorem

that              Λ 1 (Z) d → N(0, 1) under H 0 , Λ 1 (Z) d → N m 1 -m 0 √ v 0 , v 1 v 0 under H 1 . (26)
Finally, the GLRT δ 1 based on the normalized GLR Λ 1 (Z) is given by

δ 1 (Z) =            H 0 if Λ 1 (Z) < τ 1 H 1 if Λ 1 (Z) ≥ τ 1 ( 27
)
where the decision threshold τ 1 is the solution of the equation 26), the decision threshold and the power of the GLRT δ 1 can be accordingly defined as in the Theorem 1.

P H 0 Λ 1 (Z) ≥ τ 1 = α 0 . From (

Generalized Likelihood Ratio Test with Unknown Parameters

(α k , c k,1 , d k,1 )
Before designing the GLRT, the LS estimation of camera parameters (c k,1 , d k,1 ) is performed on the inspected image Z; see Section 3. The LS estimates (ĉ k,1 , dk,1 ) are asymptotically equivalent to ML estimates in large samples [START_REF] Rao | Linear models : Least Squares and Alternatives, 2nd Edition[END_REF]. Moreover, they are unbiaised and follow the asymptotic bivariate Gaussian distribution

          ĉk,1 dk,1           ∼ N                     c k,1 d k,1           ,           σ 2 c k,1 σ c k,1 d k,1 σ c k,1 d k,1 σ 2 d k,1                     , ( 28
)
where

σ 2 c k,1 , σ 2 d k,1 , σ c k,1 d k,1
denote the variance of ĉk,1 , variance of dk,1 , and covariance between ĉk,1 and dk,1 , respectively (see also discussions in Section 6.1 for this covariance matrix). The parameters (c k,1 , d k,1 ) would characterize an unknown camera model. It is require to take into account the variability of LS estimates (ĉ k,1 , dk,1 ) in the analytic establishment of performance of the GLRT.

By replacing unknown parameters (α

k , c k,1 , d k,1 ) by ( αk , ĉk,1 , dk,1 ) in the LR Λ(I k,i ) (15), the GLR Λ 2 (I k,i ) is given by Λ 2 (I k,i ) = αk 2 log ĉk,1 αk + dk,1 c k,0 αk + d k,0 - √ 2|I k,i | ĉk,1 αk + dk,1 -c k,0 αk + d k,0 . (29)
Proposition 2. Under hypothesis H j , from the Delta method [15, theorem 11.2.14], the first two moments of the GLR Λ 2 (I k,i ) can be approximated as

E H j Λ 2 (I k,i ) = m k, j (30)
Var

H j Λ 2 (I k,i ) = v k, j + β 2 k,1 α k (α k + 2) 4 α 2 k σ 2 c k,1 + σ 2 d k,1 + 2α k σ c k,1 d k,1 . (31)
Proof. of Proposition 2 is given in Appendix D.

For brevity, let denote ṽk, j = Var H j Λ 2 (I k,i ) . It can be noted that the second term in [START_REF] Mihçak | Spatially adaptive statistical modeling of wavelet image coefficients and its application to denoising[END_REF] aims to take into account the variability of LS estimates (ĉ k,1 , dk,1 ). In virtue of Lindeberg CLT, the GLR

Λ 2 (Z) = K k=1 N b i=1 Λ 2 (I k,i ) follows the Gaussian distribution under each hypothesis H j Λ 2 (Z) d → N m j , ṽ j , ( 32
)
where the expectation m j is given in [START_REF] Thai | Camera model identification based on the generalized noise model in natural images[END_REF] and the variance ṽ j is defined as

ṽ j = K k=1 N b ṽk, j . (33)
Finally, the GLRT δ 2 based on the normalized GLR Λ 2 (Z) = Λ 2 (Z)-m0 √ v0 is written as

δ 2 (Z) =            H 0 if Λ 2 (Z) < τ 2 H 1 if Λ 2 (Z) ≥ τ 2 ( 34
)
where the decision threshold τ 2 is the solution of the equation P H 0 Λ 2 (Z) ≥ τ 2 = α 0 and m0 and v0 are estimates of m 0 and ṽ0 by replacing unknown parameters (α k , c k,1 , d k,1 ) by ( αk , ĉk,1 , dk,1 ) in ( 19) and [START_REF] Lam | A mathematical analysis of the DCT coefficient distributions for images[END_REF], respectively. From the Slutsky's theorem [15, theorem 11.2.11], the decision threshold and the power of the GLRT δ 2 are given in the following theorem.

Theorem 2. When the image Z is tested against the known camera model S 0 characterized by the parameters (c k,0 , d k,0 ), the decision threshold and the power of the GLRT δ 2 are given by

τ 2 = Φ -1 (1 -α 0 ) ( 35
)
β δ 2 = 1 -Φ       m 0 -m 1 + τ 2 √ ṽ0 √ ṽ1       . ( 36
)
The statistical performance of the proposed GLRTs δ 1 and δ 2 is analytically provided. Moreover, they allow us to warrant a prescribed false alarm rate and set the decision threshold independently of camera parameters (see [START_REF] Thai | Steganalysis of jsteg algorithm based on a novel statistical model of quantized DCT coefficients[END_REF] and ( 35)). It is worth noting that the GLRT dealing with unknown parameters α k while the camera parameters (c k, j , d k, j ) are known can be interpreted as a closed hypothesis testing since the decision is made only between two known camera models S 0 and S 1 . Meanwhile, the GLRT dealing with unknown camera parameters (c k,1 , d k,1 ) becomes an open hypothesis testing telling whether the given image is acquired by camera model S 0 or not. The given image is allowed to be acquired by an unknown camera model. Therefore, two proposed tests can be straightforwardly applied, depending on the requirements of the operational context.

Numerical Experiments

In this paper, the wavelet-based denoising filter proposed in [START_REF] Mihçak | Spatially adaptive statistical modeling of wavelet image coefficients and its application to denoising[END_REF][START_REF] Lukas | Digital camera identification from sensor pattern noise[END_REF] is employed to suppress image content because of its relative accuracy and computational efficiency. Besides, the selection of homogeneous blocks requires an appropriate threshold λ. This threshold should be fixed independently of image content. The threshold λ is set at λ = 0.5.

Detection Performance on Simulated Database

The implementation of the GLRT δ 2 requires to know the covariance matrix of LS estimates (ĉ k,1 , dk,1 ). However, the ML estimates ( αk , βk ) are solved numerically, which causes a difficulty of defining their statistical properties. Thus it seems impossible to establish the covariance matrix of (ĉ k , dk ) analytically. To overcome this difficulty, it is proposed to estimate the parameters (c k , d k ) on each image from 50 images taken by the camera model S 0 since this camera model is assumed to be available. Then the empirical covariance matrix can be calculated from previous couples (ĉ k , dk ). Speaking rigorously, this is the covariance matrix characterizing the variability of the camera parameters (c k,0 , d k,0 ). By doing so, it is expected that the parameters (ĉ k,1 , dk,1 ) fall into the neighborhood of (c k,0 , d k,0 ), namely that the inspected image Z is acquired by the camera model S 0 . This leads us to exploit this covariance matrix in the implementation of the GLRT δ 2 . This step is also performed in the test with real images.

The detection performance of proposed tests is first theoretically studied on simulated database. Suppose that the camera models S 0 and S 1 are characterized by the parameters (c 0 , d 0 ) = (11.8, -3.5) and (c 1 , d 1 ) = (13.5, -4.5), respectively. These parameters correspond to frequency [START_REF] Li | Color-decoupled photo response non-uniformity for digital image forensics[END_REF][START_REF] Li | Color-decoupled photo response non-uniformity for digital image forensics[END_REF] of JPEG images taken by Canon Ixus 70 and Nikon D200 camera models in the Dresden image database [START_REF] Gloe | The Dresden image database for benchmarking digital image forensics[END_REF], respectively (see Figure 2). They are used to generate randomly 5000 vectors of 1024 and 4096 coefficients under H 0 and H 1 . Because this paper proposes to simplify the LR Λ(I k,i ) to facilitate the study, it is desirable to compare the detection performance of the LRT based on the approximating LR with the one based on the exact LR. The expectation and variance of the exact LR are calculated numerically. Moreover, it is necessary to compare the detection performance of the proposed GLRTs with the LRT since the GLRTs utilize ML estimates of unknown parameters, which may cause a loss of power. Figure 4 and Figure 5 show the detection performance of all proposed tests for 1024 and 4096 coefficients, respectively. For clarity, only regions of interest are illustrated in the figures. It is worth noting that the loss of power between the theoretical LRT and approximating LRT is negligible. Besides, a small loss of power is revealed between the GLRTs and LRT due to the estimation of unknown parameters. Nevertheless this loss of power decreases when the number of coefficients increases. It can be also noted that the loss of power between two GLRTs δ 1 and δ 2 is negligible, i.e. the variability of estimates (ĉ k,1 , dk,1 ) are well taken into account in the GLRT δ 2 . The power function of all proposed tests is perfect (e.g. β δ = 1) from 2 14 coefficients for any false alarm rate α 0 . Moreover, it is desirable to study the detection performance of the proposed tests on simulated images that follow the image processing pipeline as described in Section 2. To this end, suppose the camera models S 0 and S 1 are characterized by the parameters (ã 0 , b0 , γ 0 ) = (0.1, 2, 2.2) and (ã 1 , b1 , γ 1 ) = (0.2, 2, 2.2). These parameters are used together with the reference image lena to generate randomly 5000 images under H 0 and H 1 . The simulated images are then compressed with quality factor of 90 and 75. The detection performance of the proposed GLRTs for 1024 coefficients at frequency [START_REF] Li | Color-decoupled photo response non-uniformity for digital image forensics[END_REF][START_REF] Li | Color-decoupled photo response non-uniformity for digital image forensics[END_REF] extracted randomly from those simulated images is shown in Figure 6. As expected, a small loss of power is revealed with the decline of quality factor. It is important to remind that the proposed GLRTs are designed in the framework of hypothesis testing theory where the reference camera parameters (c k,0 , d k,0 ) under hypothesis H 0 are assumed to be known in advance. Therefore, those parameters need to be defined accurately in practice. To this end, the parameters (c k , d k ) are estimated on 50 images of the camera model S 0 and the reference parameter (c k,0 , d k,0 ) is calculated as the average of previous estimates (ĉ k , dk ). Evidently, using more images will get a better estimate but it is also less realistic. The number of 50 is a good trade-off.

τ α 0 -3 -2 -1 0 1 2 3 10 -3 10 -2 10 -1 1 Theo. FAP: α 0 = 1 -Φ(τ) δ 1 : Empirical FAP δ 2 : Empirical FAP
To highlight the relevance of the proposed GLRTs, two Canon Ixus 70 and Nikon D200 camera models of the Dresden image database [START_REF] Gloe | The Dresden image database for benchmarking digital image forensics[END_REF] are chosen to conduct experiments. The Canon Ixus 70 and Nikon D200 cameras are respectively set at H 0 and H 1 . All available JPEG images of each camera model are used in this experiment. The reference camera parameters are estimated as discussed above. The Figure 7 shows the detection performance of the GLRTs δ 1 and δ 2 for 1024 and 4096 coefficients extracted randomly at frequency [START_REF] Li | Color-decoupled photo response non-uniformity for digital image forensics[END_REF][START_REF] Li | Color-decoupled photo response non-uniformity for digital image forensics[END_REF] of natural JPEG images taken by Canon Ixus 70 and Nikon D200 camera models. We can note a similar behavior to the detection performance on simulated database. Besides, there is a small loss of power between the two GLRTs because different estimates (ĉ k,1 , dk,1 ) used in the design of the GLRT δ 2 are still influenced by image content. Nevertheless, this loss of power also decreases when the number of coefficients increases. Besides, Figure 8 illustrates detection performance of the GLRT δ 2 for 4096 coefficients randomly extracted at different frequencies. It can be noted that the detection performance decreases with the reverse zig-zag order.

Meanwhile, the Figure 9 shows the comparison between the theoretical and empirical false alarm probability, which are plotted as a function of decision threshold τ. The two proposed GLRTs δ 1 and δ 2 show an ability of guaranteeing a prescribed false alarm rate, even though there is a slight difference in some cases (typically α 0 ≤ 10 due to the influence of image content and the inaccuracy of the CLT for modeling tails.

Detection Performance on a Large Image Database

Experiments are then conducted on a large database to verify the efficiency of the proposed approach. The public Dresden image database [START_REF] Gloe | The Dresden image database for benchmarking digital image forensics[END_REF] is chosen in our experiments. Technical specifications of the cameras are shown in Table Table 3: Detection performance of SVM-based detector [START_REF] Gloe | The Dresden image database for benchmarking digital image forensics[END_REF] incorrect detection in some groups of camera models, such as (C,N1,Pe), (F,Ro), and (K,N3), is important. This may be justified due to a similarity in JPEG compression scheme used in the camera. To deal with this scenario, a second testing round involves performing the GLRT δ 1 on the camera models of conflict. The detection performance of the test δ 1 is shown in Table 5. The images are almost correctly classified. It should be noted that the GLRT δ 1 aims to give a decision rule between two different known camera models, thus the experiment on the diagonal of Table 5 is not performed. This paper also presents the SVM-based detector for comparison, which has been already performed in [START_REF] Gloe | The Dresden image database for benchmarking digital image forensics[END_REF] on the Dresden database using 46 different features to capture characteristics of different camera components of a digital camera. The detection performance of this SVM-based detector is shown in Table 3. The proposed detector δ 2 is slightly equivalent to the SVM-based detector in terms of average correct detection performance but the misclassification of the former is more severe. The PRNU-based detector [START_REF] Fridrich | Digital image forensics[END_REF] is also performed in this experiment.

This PRNU-based detector is only conducted on one device per model. Its detection performance is shown in Table 4. Overall, the two proposed detectors provide an equivalent detection performance compared with the other ones in the literature, but the latter can not allow to guarantee a prescribed false alarm probability like the proposed detectors.

Remark 1. The present paper proposes to exploit the state-of-the-ar model of DCT coefficients provided in [START_REF] Thai | Statistical model of quantized DCT coefficients: Application in the steganalysis of Jsteg algorithm[END_REF]. This model is not only more accurate than prior-art models in the literature, but also is the only one that is mathematically justified based on a statistical analysis of images' properties according to the image processing pipeline, as provided in [START_REF] Thai | Statistical model of quantized DCT coefficients: Application in the steganalysis of Jsteg algorithm[END_REF]. Therefore, relying on the proposed model allows us to capture accurately statistics of DCT coefficients as well as to analyze camera fingerprint that can be exploited for camera model identification. The function g(t) reaches the minimum at t * = |u| β 2 and its second derivative is defined by g (t) = u 2 t 3 . Consequently, the function f I (u) can be approximated as

f I (u) ≈ |u| α-1 (2β) α 2 Γ(α) exp         -|u| 2 β         . (B.8)
It can be noted that this approximating model is a special case of the GΓ model [START_REF] Chang | Image probability distribution based on generalized gamma function[END_REF] when γ = 1 (the variable γ is given in [34, Eq. ( 6)]).

Appendix C. Expectation and Variance of the LR Λ(I k,i ) under Hypothesis H j It can be noted from [START_REF] Lehmann | Testing Statistical Hypotheses, 3rd Edition[END_REF] that it is necessary to calculate the expectation and variance of the random variable |I|.

Given a constant variance σ 2 b , the random variable I is normally distributed with zero-mean and variance σ 2 b . Thus, the random variable |I| follows the half-Normal distribution [START_REF] Johnson | Continuous univariate distributions, 2nd Edition[END_REF]. Therefore, we obtain

E I|σ 2 b |I| | σ 2 b = 2 π σ b . (C.1)
Based on the law of total expectation, the mathematical expectation of |I| is given by

E I |I| = E σ 2 b E I|σ 2 b |I| | σ 2 b = 2 π E σ 2 b [σ b ] = 2 π β 1 2 Γ(α + 1 2 ) Γ(α) . (C.2)
Besides, the variance of |I| is given by

Var I |I| = E I |I| 2 -E 2 I |I| = αβ - 2β π Γ 2 (α + 1 2 ) Γ 2 (α) . (C.3)
Consequently, the expectation and variance of the LR Λ(I k,i ) under hypothesis H j can be defined by

E H j Λ(I k,i ) = α k 2 log β -1 k,1 β -1 k,0 - 2 √ π β 1 2 k, j Γ(α k + 1 2 ) Γ(α k ) β -1 k,1 -β -1 k,0 (C.4) Var H j Λ(I k,i ) = 2 β -1 k,1 -β -1 k,0 2 × α k β k, j - 2β k, j π Γ 2 (α k + 1 2 ) Γ 2 (α k ) . (C.5)
Appendix D. Asymptotic Expectation and Variance of the GLR Λ 2 (I k,i ) under Hypothesis H j It is assumed that the variance of the estimate αk is negligible when the number of coefficients N b is very large.

Thus it is proposed to treat the estimate αk as a constant α k . Besides, since the estimates (ĉ k,1 , dk,1 ) are consistent, the asymptotic mathematical expectation of the GLR Λ 2 (I k,i ) under hypothesis H j does not change, i.e. E H j Λ 2 (I k,i ) = m k, j .

Meanwhile, the variance of the GLR Λ 2 (I k,i ) needs to take into account the variability of the estimates (ĉ k,1 , dk,1 ).

Based on the definitions of mathematical expectation and variance, we have

E
+ 2 β -1 k,1 -β -1 k,0 2 α k β k, j - 2β k, j π Γ 2 (α k + 1 2 ) Γ 2 (α k ) . (D.5)

Figure 2 :Figure 3 :

 23 Figure 2: Estimated parameters (α, β) at frequency (8, 8) of natural JPEG images issued from Canon Ixus 70 and Nikon D200 camera models.

  involves the JPEG images with different imaged scenes, different camera settings, different devices per model, and different environmental conditions. This indicates that the parameters (c k , d k ) remain similar under those conditions.

δ 2 Figure 4 :

 24 Figure 4: Detection performance of proposed tests on simulated vectors with 1024 coefficients.

Figure 5 :

 5 Figure 5: Detection performance of proposed tests on simulated vectors with 4096 coefficients.

Figure 6 :

 6 Figure 6: Detection performance of proposed tests for 1024 coefficients at frequency (8, 8) extracted randomly from simulated images with different quality factors.

Figure 7 :

 7 Figure 7: Detection performance of proposed tests for different number of coefficients at frequency (8, 8) of natural JPEG images taken by Ixus 70 and Nikon D200 camera models.

Figure 8 :

 8 Figure 8: Detection performance of the GLRT δ 2 for 4096 coefficients at different frequencies of natural JPEG images taken by Canon Ixus 70 and Nikon D200 camera models.

Figure 9 :

 9 Figure 9: Comparison between the theoretical false alarm probability (FAP) and the empirical FAP, plotted as a function of decision threshold τ.

Table 1 :

 1 Camera Model Used in Experiments 6.2. Detection Performance on Two Canon Ixus 70 and Nikon D200 Camera Models

	Camera Model		No. devices	Resolution	No. images
	Canon Ixus 70	C	3	3072 × 2304	350
	Fujifilm FinePix J50	F	3	3264 × 2448	350
	Kodak M1063	K	5	3664 × 2748	550
	Nikon Coolpix S710	N1	5	4352 × 3264	550
	Nikon D200	N2	2	3872 × 2592	250
	Nikon D70	N3	2	3008 × 2000	250
	Pentax Optio A40	Pe	4	4000 × 3000	450
	Praktica DCZ 5.9	Pr	5	2560 × 1920	550
	Ricoh Capilo GX100	Ri	5	3648 × 2736	550
	Rollei RCP-7325XS	Ro	3	3072 × 2304	350
	Sony DSC-H50	S	2	3456 × 2592	250
		11	39		4450

Table 2 :

 2 Detection performance of proposed detector δ 2 (the symbol * represents values smaller than 2%)

  Firstly, the GLRT δ 2 is used to verify whether a given image is acquired by the camera model of interest. The decision threshold τ 2 is given by the Theorem 2 corresponding to the false alarm rate α 0 = 10 -5 . If the normalized GLR Λ 2 (Z) is smaller than the decision threshold τ 2 , the hypothesis H 0 is accepted, i.e. the given image is taken from the camera model of interest. On the contrary, the hypothesis H 1 is accepted. It is proposed to use the last 21 high frequencies for the test. The detection performance of the test δ 2 is shown in Table2. In this table, each camera model is considered as hypothesis H 0 (row) and all images (column) are tested against H 0 . The values in the table indicate the percentage of images that are detected taken by the camera model H 0 . The table in this paper is not used in the same way as in the classification in which the sum for each class yields 100%. The inspected image is brought into the binary testing of the known camera model H 0 against the others, thus the sum of a class may not yield 100%. It could lead to a scenario that an image is declared taken by at least two camera models. It can be noted from Table2that the

			Cn2	N1	Pe	F	Ro	K	N3	N2	Pr	Ri	S	Avg.
		Cn2	98.26	*	*	*	*	*	*	2.61	*	*	*
		N1	*	100	3.97	*	*	*	*	*	*	*	*
		Pe	*	*	90.86	*	*	*	*	*	*	*	*
		F	*	*	*	99.94	*	*	*	*	*	*	*
		Ro	*	*	*	*	99.79	*	*	*	*	*	*
	H 0	K	*	*	*	*	*	99.19	*	*	*	*	*
		N3	*	*	*	*	*	*	84.57	*	*	*	*
		N2	*	*	*	*	*	*	*	93.73	*	*	*
		Pr	*	*	*	*	*	*	*	*	98.42	*	*
		Ri	*	*	*	*	*	*	*	*	*	100	*
		S	*	*	*	*	*	*	*	*	*	*	98.23
														96.63

1, see more details in

[START_REF] Gloe | The Dresden image database for benchmarking digital image forensics[END_REF]

. The database covers different devices per camera model, different imaged scenes, different camera settings and different environmental conditions. All images are acquired with the highest available JPEG quality setting and maximum available resolution. For each camera model, the set contains 50 images per camera model for estimation of reference camera parameters and 100 images per device for testing, which are randomly selected from the Dresden database.

  ) expg(t * ) . (B.5) As described in[START_REF] Thai | Statistical model of quantized DCT coefficients: Application in the steganalysis of Jsteg algorithm[END_REF], the DCT coefficient model f I (u) is rewritten as follows

	Cn2	Cn2 100	N1 * * I ≈ exp -g(t * ) Pe F Ro K N3 * * * * exp -g (t * ) 2	N2 (t -t * ) 2 dt. Pr	Ri	S	Avg.	(B.3)
	This integral takes the form of Gaussian integral. We derive
			I ≈		2π |g (t * )|	exp -g(t * ) .	(B.4)
	H 0
	A generalization has been made in [30] with an arbitrary function h(t)
		I =	h(t) exp -g(t) dt ≈ h(t f I (u) = 2π |g (t * )| ∞ 1 √ 0 2πβ α Γ(α) h(t) exp -g(t * ) dt	(B.6)
	where
			g(t) =	t β	+	u 2 2t	and h(t) = t α-3 2 .	(B.7)

*

  H j ĉk,1 αk + dk,1 = β -1 Var H j ĉk,1 αk + dk,1 = α 2 k σ 2 c k,1 + σ 2 d k,1 + 2α k σ c k,1 d k,1 . (D.2)Subsequently, from the Delta method [15, theorem 11.2.14], we derive thatVar H j log ĉk,1 αk + dk,1 c k,0 αk + d k,0 = Var H j log ĉk,1 αk + dk,1 = Var H j ĉk,1 αk + dk,1 E 2 H j ĉk,1 αk + dk,1 = β 2 k,1 α 2 k σ 2 c k,1 + σ 2 d k,1 + 2α k σ c k,1 d k,1 , (D.3)andVar H j ĉk,1 αk + dk,1c k,0 αk + d k,0 = Var H j ĉk,1 αk + dk,1 = Var H j ĉk,1 αk + dk,1 4E H j ĉk,1 αk + dk,1 + σ 2 d k,1 + 2α k σ c k,1 d k,1 . (D.4) Var H j Λ 2 (I k,i ) = α 2 k 4 Var H j log ĉk,1 αk + dk,1 c k,0 αk + d k,0 + 2Var H j ĉk,1 αk + dk,1c k,0 αk + d k,0 E H j |I| 2 ĉk,1 αk + dk,1c k,0 αk + d k,0 Var H j |I| 2 + σ 2 d k,1 + 2α k σ c k,1 d k,1

	= c k,1 + 2E 2 β k,1 4 α 2 k σ 2 H j = β 2 k,1 α k (α k + 2) 4	k,1 c k,1 k σ 2 α 2	(D.1)

The term source means an individual camera instance, a camera model, or a camera brand.

Acknowledgements

The authors would like to thank the Digital Data Embedding Laboratory, University of Binghamton, especially J.

Fridrich, for the use of the High Performance Computing server, allowing a considerable reduction in execution time of our algorithm on the large image database.

100 96.99 Table 4: Detection performance of PRNU-based detector [START_REF] Fridrich | Digital image forensics[END_REF]

Table 5: Detection performance of proposed detector δ 1

Conclusion

The goal of this paper is to design a statistical test for camera model identification from JPEG images within hypothesis testing framework. The approach is based on the state-of-the-art model of DCT coefficients to capture their statistical difference, which jointly results from different sensor noises and in-camera processing algorithms.

The parameters (c, d) characterizing the simplistic linear relation between α and β -1 , which are two parameters of the DCT coefficient model, are proposed to be exploited as camera fingerprint for camera model identification. Based on the parametric model of DCT coefficients, this paper studies the most powerful LRT and proposes two GLRTs that can be straightforwardly applied in practice. The strength of the proposed approach is that statistical performance of the tests can be analytically established as well as they can warrant a prescribed false alarm rate while ensuring a high detection performance.

Appendix A. Relation between the Parameters ( ã, b, γ) and (α p,q , β p,q )

Given an arbitrary image Z, the DCT operation is performed in each 8 × 8 block of Z as follows

where z m,n denotes a pixel within a 8 × 8 block of Z, 0 ≤ m ≤ 7, 0 ≤ n ≤ 7 and I p,q denotes the two-dimensional DCT coefficient and

For the sake of simplification, it is assumed that the pixels are independent and identically distributed within each In fact, the proposed model of DCT coefficients (3) is based on the doubly stochastic model given in [START_REF] Lam | A mathematical analysis of the DCT coefficient distributions for images[END_REF] that takes into account the statistics of DCT coefficients in a block whose variance is constant and the variability of block variance in a natural image. Given a constant block variance σ 2 b , the AC coefficient I may be approximately distributed as zero-mean Gaussian in virtue of the CLT. Besides, the block variance σ 2 b can be approximately modeled by the Gamma distribution G(α, β). As the pdf of I is symmetric, the odd moments vanish. Based on the law of total expectation, the short calculation shows that

On the other hand, it follows from [START_REF] Thai | Generalized signal-dependent noise model and parameter estimation for natural images[END_REF] that

Consequently, one derives that

S p,q + 1 α p,q β p,q + 3β p,q = 4S p,q µ 2 z m,n .

(A.11)

Resolving this system of equations, the relation between α p,q and β -1 p,q is given as

(S p,q + 1)α p,q + 3 4 b 2 γ S p,q × γ 2 α p,q b (S p,q + 1)α p,q + 3 4S p,q

This relation is too complicated to exploit. Therefore, it is proposed to employ the polynomial expansion and only keep the first two terms β -1 p,q = c p,q α p,q + d p,q , (A. [START_REF] Choi | Source camera identification by JPEG compression statistics for image forensics[END_REF] where the parameters (c p,q , d p,q ) depend on the parameters (ã, b, γ). Numerical experiments show that this simplified equation sufficiently characterizes the relation between the parameters (ã, b, γ) and (α p,q , β p,q ) (see Figure 1 and Figure 2).

Appendix B. Laplace's Approximation of DCT Coefficient Model

Let us briefly describe the idea behind the Laplace's approximation [START_REF] Tierney | Accurate approximations for posterior moments and marginal densities[END_REF]. The Laplace's method aims to provide an approximation for integrals of the form