An Adaptive Method for Camera Attribution under Complex Radial Distortion Corrections

TL;DR

This paper introduces an adaptive method that effectively handles complex radial distortion corrections in images, improving camera attribution accuracy and efficiency over existing approaches.

Contribution

The proposed algorithm divides images into concentric annuli to better handle sophisticated radial distortions, and introduces the CPCE statistic for efficient early stopping.

Findings

Improves camera attribution accuracy on radially distorted images.

Reduces computational cost compared to previous methods.

Effective on complex out-camera correction software like Adobe Lightroom.

Abstract

Radial correction distortion, applied by in-camera or out-camera software/firmware alters the supporting grid of the image so as to hamper PRNU-based camera attribution. Existing solutions to deal with this problem try to invert/estimate the correction using radial transformations parameterized with few variables in order to restrain the computational load; however, with ever more prevalent complex distortion corrections their performance is unsatisfactory. In this paper we propose an adaptive algorithm that by dividing the image into concentric annuli is able to deal with sophisticated corrections like those applied out-camera by third party software like Adobe Lightroom, Photoshop, Gimp and PT-Lens. We also introduce a statistic called cumulative peak of correlation energy (CPCE) that allows for an efficient early stopping strategy. Experiments on a large dataset of in-camera and out-camera radially corrected images show that our solution improves the state of the art in terms of both accuracy and computational cost.

Tables

Table 1. Table I : Thresholds required to achieve FPR=0.05 and FPR=0.01 and corresponding TPRs for different methods and variants.

	${\tau }_{0.05}$	${\tau }_{0.01}$	TPR0.05	TPR0.01
$\overrightarrow{\text{ID}}$, Lin	98.86	112.71	0.96	0.94
$\overrightarrow{\text{DI}}$, Lin	90.01	105.63	0.97	0.95
$\overrightarrow{\text{ID}}$, Cub	73.48	90.28	0.99	0.98
$\overrightarrow{\text{DI}}$, Cub	71.13	84.75	0.99	0.98
Inv, Lin	97.66	111.72	0.93	0.91
Dir, Lin	90.01	105.63	0.95	0.92
Inv, Cub	73.48	90.29	0.98	0.98
Dir, Cub	71.13	84.57	0.99	0.98
2W, Cub	71.12	84.34	0.99	0.99
[8]	4.81	8.36	0.86	0.82
[8] no DS	7.65	10.54	0.94	0.93
[9]	2.83	5.50	0.78	0.73

Table 2. Table II : TPRs of the different variants of our method and [ 8 ] (with and without DS), and the average time consumed to declare a match for specific devices/software. In bold we highlight the fastest and most accurate variants of our method.

	$\overrightarrow{\text{ID}}$, Lin		$\overrightarrow{\text{DI}}$, Lin		$\overrightarrow{\text{ID}}$, Cub		$\overrightarrow{\text{DI}}$, Cub		Inv, Cub		Dir, Cub		2W, Cub		[8] no DS		[8]
	TPR	time	TPR	time	TPR	time	TPR	time	TPR	time	TPR	time	TPR	time	TPR	time	TPR	time
GIMP [$3456\times 5184$]	0.98	220.1	0.98	303.6	0.98	193.9	0.99	204.8	0.98	182	0.98	196.2	0.99	326.4	0.97	962.3	0.97	85.3
LIGHTROOM [$3456\times 5184$]	0.93	237.8	0.95	330.3	0.94	288.4	0.97	309.1	0.94	274.6	0.97	295.1	0.98	467.4	0.6	930.8	0.44	87.9
LIGHTROOM* [$3456\times 5184$]	0.95	137.5	0.92	222.6	0.99	158.8	0.98	126.5	0.99	127.7	0.98	129.7	0.98	220.5	0.51	930.9	0.4	91.6
PHOTOSHOP [$3456\times 5184$]	0.97	151.2	0.98	159.5	0.97	165.4	0.97	176.2	0.97	170.7	0.97	167.3	0.97	293	0.96	974.6	0.91	107.2
PT LENS [$3456\times 5184$]	0.97	181.3	0.99	247.3	0.98	248.9	0.99	253.6	0.98	256.3	1	238.2	1	412.4	0.91	918.3	0.8	100.3
S9100 [$3000\times 4000$]	0.99	114.2	0.98	117.3	0.99	72.9	0.99	99.8	0.99	92.9	1	92.5	1	160.3	1	598.3	0.98	81.7
SX210 [$3240\times 4320$]	0.96	256.8	0.97	262.5	1	127.5	1	139.1	1	121.8	1	126.8	1	210.3	1	723.4	0.99	81.3
SX230 [$1584\times 2816$]	0.87	33.8	0.9	29	1	23	1	24.2	0.98	41	0.98	41.8	0.98	70.5	0.98	229.0	0.83	25.1
SX40 [$2664\times 4000$]	1	105.1	0.96	145	1	93.4	1	103.4	1	88.6	1	93.6	1	154	1	526.3	0.98	56.3
ZS7 [$1920\times 2560$]	0.76	42.7	0.81	35	0.98	27.7	0.98	29.8	0.93	50.5	0.99	46.9	0.99	79.6	0.98	255.2	0.79	27.8

Equations

$$\rho(\mbox{$\bf X$},\mbox{$\bf Y$})=\frac{\langle\mbox{$\bf X$}-\bar{\mbox{$\bf X$}},\mbox{$\bf Y$}-\bar{\mbox{$\bf Y$}}\rangle}{||\mbox{$\bf X$}-\bar{\mbox{$\bf X$}}||\cdot||\mbox{$\bf Y$}-\bar{\mbox{$\bf Y$}}||},$$

$$\mbox{$\bf I$}=\mbox{$\bf I$}_{0}+\mbox{$\bf I$}_{0}\circ\mbox{$\bf K$}+\mathbf{\Theta},$$

$$\mbox{$\bf W$}^{(l)}=\mbox{$\bf I$}^{(l)}-F(\mbox{$\bf I$}^{(l)}),\ \ l=1,\cdots,L.$$

$${\hat{\mbox{$\bf K$}}}=\left(\sum_{l=1}^{L}\mbox{$\bf I$}^{(l)}\circ\mbox{$\bf W$}^{(l)}\right)\circ\left(\sum_{l=1}^{L}\mbox{$\bf I$}^{(l)}\circ\mbox{$\bf I$}^{(l)}\right)^{\circ-1}.$$

$$\text{PCE}(\hat{\mbox{$\bf K$}}^{\prime},\mbox{$\bf W$})=\frac{\text{sgn}(\rho(\hat{\mbox{$\bf K$}}^{\prime},\mbox{$\bf W$}))\cdot\rho^{2}(\hat{\mbox{$\bf K$}}^{\prime},\mbox{$\bf W$})}{\frac{1}{|{\mathcal{I}}\backslash{\mathcal{S}}|}\sum_{\mbox{$\bf s$}\in{\mathcal{I}}\backslash{\mathcal{S}}}\rho^{2}(\hat{\mbox{$\bf K$}}^{\prime},C(\mbox{$\bf W$},\mbox{$\bf s$}))},$$

$$\text{PCE}(\hat{\mbox{$\bf K$}}^{\prime},\mbox{$\bf W$})=\frac{\text{ssq}(\langle\hat{\mbox{$\bf K$}}^{\prime},\tilde{\mbox{$\bf W$}}\rangle)}{\frac{1}{|{\mathcal{I}}\backslash{\mathcal{S}}|}\sum_{\mbox{$\bf s$}\in{\mathcal{I}}\backslash{\mathcal{S}}}\langle\hat{\mbox{$\bf K$}}^{\prime},C(\tilde{\mbox{$\bf W$}},\mbox{$\bf s$})\rangle^{2}},$$

$$T_{\alpha }:{\mathbb{R}}^2$$

$$(x,y)$$

$$x^{\prime }=x_p+(x-x_p)(1+\alpha r^2);$$

$$y^{\prime }=y_p+(y-y_p)(1+\alpha r^2),$$

$$r^{\prime }=T_{\alpha }(r)=r(1+\alpha r^2).$$

$$r^{\prime }=T_{\mathbf{\alpha }}(r)=r\left(1+\sum _{i=1}^n{\alpha }_i{r}^{2i}\right),$$

$$\text{PCE}_{\mathsf{dir}}({\bm{\alpha}})\doteq\text{PCE}(T_{\bm{\alpha}}(\hat{\mbox{$\bf K$}}^{\prime}),\mbox{$\bf W$}),$$

$$\text{PCE}_{\mathsf{inv}}({\bm{\alpha}})\doteq\text{PCE}(\hat{\mbox{$\bf K$}}^{\prime},T_{\bm{\alpha}}^{-1}(\mbox{$\bf W$})).$$

$${\text{PCE}}_{\mathsf{max}}(\mathbf{\alpha })\doteq \max \{{\text{PCE}}_{\mathsf{dir}}(\mathbf{\alpha }),{\text{PCE}}_{\mathsf{inv}}(\mathbf{\alpha })\}.$$

$${\text{PCE}}_{\mathsf{max}}^{\ast }\doteq \underset{\mathbf{\alpha }\in \mathcal{A}}{\max }{\text{PCE}}_{\mathsf{max}}(\mathbf{\alpha }).$$

$$r=T_{\alpha }^{-1}(r^{\prime })=r^{\prime }(1-\alpha {r^{\prime }}^2+3{\alpha }^2{r^{\prime }}^4+O({r^{\prime }}^6)).$$

$${\mathcal{R}}_k\doteq \{(u,v)\in {\mathbb{R}}^2:r_k^2\le u^2+v^2<(r_k+{\Delta }_k{)}^2\}.$$

$$r^{\prime }=T_{\alpha (r)}(r)=r(1+\alpha (r)\cdot r^2).$$

$$T_{{\alpha }_k,c}(r)$$

$$T_{{\alpha }_k,l}(r)$$

$$T_{{\alpha }_k,c}^{-1}(r^{\prime })$$

$$T_{{\alpha }_k,l}^{-1}(r^{\prime })$$

$${\mathcal{P}}_k\doteq \left(D_2\cdot {\mathcal{R}}_k\right)\cap \mathcal{I},k=1,\cdots ,L,$$

$${\mathcal{Q}}_{k,\mathsf{inv}}({\alpha }_k)=D_2\cdot T_{{\alpha }_k}\left(\left[\left(D_2\cdot T_{{\alpha }_k}^{-1}({\mathcal{P}}_k/D_2)\right)\cap \mathcal{B}\right]/D_2\right).$$

$${\mathcal{Q}}_{k,\mathsf{dir}}({\alpha }_k)=D_2\cdot T_{{\alpha }_k}^{-1}\left(\left[\left(D_2\cdot T_{{\alpha }_k}({\mathcal{P}}_k/D_2)\right)\cap \mathcal{B}\right]/D_2\right).$$

$${\mathsf{\Phi }}_{k,\mathsf{inv}}({\alpha }_k)$$

$${\mathsf{E}}_{k,\mathsf{inv}}({\alpha }_k)$$

$${\alpha }_k^{\ast }=\arg \underset{{\alpha }_k\in {\mathcal{A}}_k}{\max }{\phi }_{k,\mathsf{inv}}({\alpha }_k),$$

$${\phi }_{k,\mathsf{inv}}({\alpha }_k)\doteq \frac{{\mathsf{\Phi }}_{k,\mathsf{inv}}({\alpha }_k)}{{\mathsf{E}}_{k,\mathsf{inv}}({\alpha }_k)}.$$

$$k_0=\arg \underset{k=1,\cdots ,L}{\max }{\phi }_{k,\mathsf{inv}}(0).$$

$$\hat{\alpha}_{k}=\mbox{$\bf u$}_{k-1}^{T}\bm{\beta}_{k-1},$$

Full text

I Introduction

During the past years, camera fingerprints based on the Photo Response Non-Uniformity (PRNU) have gained broad popularity in forensic applications thanks to their ability to identify the device that captured a certain image. The PRNU is a multiplicative spatial pattern that owes its uniqueness to manufacturing imperfections that cause sensor elements to have minute area and substrate material differences that make them capture different amounts of energy even under a perfectly uniform flat field [1]. Applications of the PRNU in multimedia forensics go beyond camera identification from images [2] or videos [3], as they have also been used in detecting inconsistencies that reflect image manipulations [4].

Unfortunately, the fact that the PRNU can be accurately modeled as a white random process explains its sensitivity to geometric transformations that alter the image coordinates. Unless those transformations are reverted, standard detection statistics will perform poorly as they are roughly based on cross-correlations that yield very small values under grid misalignment. In the literature several methods have been proposed to deal with those spatial transformations, including digital zoom [5], video stabilization [6], high dynamic range (HDR) processing [7], and radial distortion corrections [8, 9]. It is in the context of the latter that we have developed the methodology presented in this work.

Radial distortion correction aims at digitally removing the distortion introduced by the camera lens. This kind of processing is becoming more pervasive as devices increase their computing capabilities; in-camera correction is now common in compact models, tablets and smartphones. On the other hand, out-camera corrections can be performed with powerful software like Adobe Lightroom, which are able to invert distortions almost perfectly by matching the model of the lens mounted on the camera. This is not done by applying conventional radial distortion models such as barrel or pincushion but by making use of complex models (i.e., with a large number of parameters). As a consequence, existing methods [8, 9] relying on models with at most two parameters will only partially succeed in dealing with camera attribution under these complex out-camera processing. Increasing the number of model parameters often constitutes an undesirable path because reverting the distortion corrections entails a grid search whose computational load grows exponentially with the number of unknowns.

In this work we propose a novel approach to PRNU-based camera attribution under radial corrections that is able to deal with complex models without significantly increasing the computational burden. The main idea is to divide the image under test and the PRNU into a series of concentric annuli that are thin enough to be locally describable with a simple (i.e., linear or cubic) distortion model which allows for an equally simple inverse transformation. The annuli are traversed sequentially by keeping track of the cumulative peak-to-correlation energy ratio (CPCE), which is a statistic introduced in this work and used to decide whether the radially corrected test image contains the reference PRNU. In fact, the sequential nature of the procedure makes it possible to implement an early stopping strategy to declare a match without having to process all the annuli and thus saving computational time. Another key feature of our method is adaptivity: instead of carrying out a wide-interval search for the distortion parameters describing each annulus, an adaptive Least-Mean-Squares-like predictor updates the parameters of the previously processed annulus in order to narrow down the current parameter search. This leads to a large computational efficiency without giving up flexibility. In order to steer the search we propose and justify mathematically a new objective function.

Different variants of our method are evaluated in terms of accuracy and speed. We compare our method with the state of the art in [8] and [9] on a large dataset composed of images taken with: 1) compact devices and radially corrected in-camera, and 2) a reflex camera and radially corrected out-camera using different software tools. Our results show considerable performance improvements, especially on low-resolution images and in presence of complex radial corrections.

The rest of the paper is organized as follows: Sect. II provides the mathematical background and formulates the addressed problem. Sect. III discusses the relevant state of the art. Sect. IV is devoted to discussing the proposed method which in Sect. V is validated and compared with [8] and [9]. Finallt, Sect.VI presents our conclusions.

II Problem Formulation and Modeling

II-A Notation

In this paper we will consider gray-scale images (the extension to color images being straightforward). Bi-dimensional signals will be denoted with boldface. For every such signal, a domain ${\mathcal{S}}\subset{\mathbb{Z}}^{2}$ will be specified; for instance, a signal $\bf X$ with domain ${\mathcal{S}}_{X}$ is a collection of values $X_{i,j}\in{\mathbb{R}}$ defined for all locations $(i,j)\in{\mathcal{S}}_{X}$. For the case of images of size $M\times N$, the original domain is ${\mathcal{I}}=\{1,\cdots,M\}\times\{1,\cdots,N\}\subset{\mathbb{Z}}^{2}$; however, we will often find ourselves working with domains that are subsets of ${\mathcal{I}}$. We will denote by $D_{2}$ half of the diagonal of domain ${\mathcal{I}}$ measured in pixels. Notice that the set ${\mathcal{I}}$ can be expressed as ${\mathcal{I}}={\mathcal{B}}\cap{\mathbb{Z}}^{2}$, with ${\mathcal{B}}\subset{\mathbb{R}}^{2}$ denoting the image bounding box.

The inner product of two signals $\bf X$ and $\bf Y$ with respective domains ${\mathcal{S}}_{X}$ and ${\mathcal{S}}_{Y}$ can be defined by extending the Frobenius product of matrices as \langle\mbox{\bf X},\mbox{\bf Y}\rangle\doteq\sum_{(i,j)\in{\mathcal{S}}}X_{i,j}Y_{i,j}, where ${\mathcal{S}}={\mathcal{S}}_{X}\cap{\mathcal{S}}_{Y}$ is assumed to be non-empty. The Frobenius norm of $\bf X$ with domain ${\mathcal{S}}_{X}$ induced by this inner product is ||\mbox{\bf X}||\doteq\langle\mbox{\bf X},\mbox{\bf X}\rangle=\sum_{(i,j)\in{\mathcal{S}}_{X}}X_{i,j}^{2}. The product of signals $\bf X$ and $\bf Y$, denoted by \mbox{\bf X}\circ\mbox{\bf Y}, is the element-wise product, i.e., (\mbox{\bf X}\circ\mbox{\bf Y})_{i,j}=X_{i,j}\cdot Y_{i,j} and is defined for all $(i,j)\in{\mathcal{S}}_{X}\cap{\mathcal{S}}_{Y}$. The multiplicative inverse of $\bf X$ is denoted by \mbox{\bf X}^{\circ-1} and is such that (\mbox{\bf X}^{\circ-1})_{i,j}=X_{i,j}^{-1}. For a signal $\bf X$ with domain ${\mathcal{S}}_{X}$, we denote by \bar{\mbox{\bf X}} a constant signal with the same support as $\bf X$ and whose value is the sample mean $\sum_{(i,j)\in{\mathcal{S}}_{X}}X_{i,j}/|{\mathcal{S}}_{X}|$, where $|{\mathcal{S}}_{X}|$ denotes the cardinality of ${\mathcal{S}}_{X}$. The normalized cross-correlation (NCC) between $\bf X$ and $\bf Y$ is defined as

[TABLE]

with the inner product and norms defined as above.

Given a signal $\bf X$ with rectangular domain ${\mathcal{I}}$ and a vector \mbox{\bf s}=(s_{1},s_{2})\in{\mathbb{Z}}^{2}, we denote by C(\mbox{\bf X},\mbox{\bf s}) the cyclic shift of $\bf X$ by vector $\bf s$, so that the $(i,j)$th component of C(\mbox{\bf X},\mbox{\bf s}) is $X_{(i+s_{1})\text{mod}M,(j+s_{2})\text{mod}N}$. Note that the domain of C(\mbox{\bf X},\mbox{\bf s}) is also ${\mathcal{I}}$. Finally, the all-zeros image is denoted by $\bf 0$.

II-B PRNU estimation

As previously indicated, the PRNU is a multiplicative noise-like signal that serves as a sensor fingerprint [1],[5]. Because the PRNU is a very weak signal, it is necessary to separate it from both the true image and other noise components. If \mbox{\bf I}_{0} denotes the image in absence of noise, and $\bf K$ is the PRNU, it is possible to derive the following simplified model [10]:

[TABLE]

where $\mathbf{\Theta}$ is uncorrelated with both \mbox{\bf I}_{0} and $\bf K$, and summarizes noise components of different nature, and all signals are defined over ${\mathcal{I}}$. The fingerprint $\bf K$ of a camera can be extracted from $L$ images \mbox{\bf I}^{(l)}, $l=1,\cdots,L$, taken with the camera under analysis. Let \mbox{\bf W}^{(l)} denote the noise residual obtained by applying a generic denoising filter $F(\cdot)$ to the $i$th image \mbox{\bf I}^{(l)}, as

[TABLE]

In all our reported experiments, we have used Mihcak’s wavelet-based denoiser [11] for it yields an excellent trade-off between performance and complexity. Then, the PRNU can be estimated as follows [10]:

[TABLE]

The estimate so obtained is customarily post-processed to remove some systematic artifacts that are present in most cameras. Here, we will follow [10] and apply a mean-removal operation by columns and rows, and a Wiener filter in the DFT aimed at removing periodic spatial artifacts. For color images, the fingerprints are estimated separately for the RGB channels and then linearly combined into gray-scale as in [12]. Given an image under investigation $\bf I$ and its corresponding residual \mbox{\bf W}\doteq\mbox{\bf I}-F(\mbox{\bf I}) a binary hypothesis test can be formulated to decide whether $\bf I$ contains a certain PRNU \mbox{\bf K}^{\prime} for which an estimate \hat{\mbox{\bf K}}^{\prime} is available. We will denote the null hypothesis of this test (i.e., $\bf I$ does not contain {\mbox{\bf K}}^{\prime}) by $H_{0}$ and the alternative (i.e., $\bf I$ contains {\mbox{\bf K}}^{\prime}) by $H_{1}$.

The most popular decision statistic for the test is the Peak-to-Correlation Energy ratio (PCE) which computes the peak cross correlation between the test image residual $\bf W$ and the estimated PRNU \hat{\mbox{\bf K}}^{\prime} from the candidate camera, and normalizes it by an estimate of the correlation noise under $H_{0}$ [12].

For non-cropped images the PCE simplifies to

[TABLE]

where, following the improvement proposed in [13], we have included the sign of the NCC to exclude negative values that would be never expected under $H_{1}$. In (5) $\mathcal{S}$ is a cyclic exclusion neighborhood of $(0,0)$ of small size (e.g., $11\times 11$ pixels) to avoid contamination from cross-correlation peaks when estimating the cross-correlation noise when $H_{1}$ holds. Noticing that for every $\bf s$, ||C(\mbox{\bf W},\mbox{\bf s})-{\overline{C(\mbox{\bf W},\mbox{\bf s})}}||=||\mbox{\bf W}-\bar{\mbox{\bf W}}||, and letting \tilde{\mbox{\bf W}}\doteq\mbox{\bf W}-\bar{\mbox{\bf W}}, (5) can be alternatively written as

[TABLE]

where we have assumed that the mean of \hat{\mbox{\bf K}}^{\prime} is zero due to the zero-meaning operation discussed above, and the signed-squared function $\text{ssq}(\cdot)$ is such that $\text{ssq}(x)\doteq\text{sgn}(x)\cdot x^{2}$.

II-C Lens Distortion Models

To describe radially symmetric barrel/pincushion distortions we adopt the same models presented in [8] and explained in [14, 15, 16]. If we denote the coordinates before and after the radial distortion by $(x,y)$ and $(x^{\prime},y^{\prime})$, the invertible geometrical mapping $T_{\alpha}$ is given by

[TABLE]

where

[TABLE]

[TABLE]

and $(x_{p},y_{p})$ is the optical center of the image and $r^{2}\doteq[(x-x_{p})^{2}+(y-y_{p})^{2}]/D_{2}^{2}$ is the normalized squared radial distance from point $(x,y)$ to the optical center. This normalization by $D_{2}^{2}$ is for convenience, so that $r=1$ corresponds to half of the image diagonal [8]. Parameter $\alpha\in\mathbb{R}$ in (8-9) models the type of radial distortion: $\alpha>0$ for pincushion distortion, and $\alpha<0$ for barrel distortion. Alternatively, given $(x_{p},y_{p})$ and assuming that $T_{\alpha}(x_{p},y_{p})=(x_{p},y_{p})$, the transformation can be written in normalized polar coordinates. Since the phase is preserved under $T_{\alpha}(\cdot)$, with a slight abuse of notation we will drop the phase component and sometimes write the radial transformation as $T_{\alpha}:\mathbb{R}^{+}\cup\{0\}\rightarrow\mathbb{R}^{+}\cup\{0\}$ such that

[TABLE]

More complex radial corrections [17], [9] can be expressed through an $n$th order model:

[TABLE]

where $\bm{\alpha}\doteq[\alpha_{1},\cdots,\alpha_{n}]^{T}$ is a real parameter vector.

Again, with some abuse of notation, and following [8], given a signal $\bf X$ with domain ${\mathcal{S}}_{X}$, the mapping \mbox{\bf Y}=T_{\bm{\alpha}}(\mbox{\bf X}) is produced as follows. Let \mbox{\bf X}^{\prime} be the signal with domain ${\mathcal{S}}_{X^{\prime}}=T_{\bm{\alpha}}({\mathcal{S}}_{X})$ such that, for every $(u,v)\in{\mathcal{S}}_{X}$, and with $(u^{\prime},v^{\prime})=T_{\bm{\alpha}}(u,v)$, $X^{\prime}_{u^{\prime},v^{\prime}}=X_{u,v}$. Then, given an ouput domain ${\mathcal{S}}_{Y}$, the signal \mbox{\bf Y}=T_{\bm{\alpha}}(\mbox{\bf X}) is obtained by interpolating the signal \mbox{\bf X}^{\prime} defined on ${\mathcal{S}}_{X}^{\prime}$ at the points in ${\mathcal{S}}_{Y}$. Of course, precautions must be taken when specifying ${\mathcal{S}}_{Y}$ so that the interpolation is computable at all points in ${\mathcal{S}}_{Y}$. This aspect will be made clearer in Sect. IV, when we present our method.

II-D Direct and inverse approaches to PCE computation

When the image under analysis has been subjected to a radial distortion correction, the statistic \text{PCE}(\hat{\mbox{\bf K}}^{\prime},\mbox{\bf W}) is expected to perform poorly under $H_{1}$ in the hypothesis test, because the grids supporting \hat{\mbox{\bf K}}^{\prime} and $\bf W$ will not coincide (recall that the PRNU has a very narrow spatial autocorrelation function).

The approach explored in [8] is to take into account the distortion correction when computing the PCE. If the parameter vector $\bm{\alpha}$ of the radial mapping is known, there are essentially two possibilities, which we will term direct and inverse. In the direct approach, the candidate PRNU \hat{\mbox{\bf K}}^{\prime} is transformed in order for its grid to match that of $\bf W$. Then, the test statistic becomes

[TABLE]

where the domain $\mathcal{I}_{T}$ of T_{\bm{\alpha}}(\hat{\mbox{\bf K}}^{\prime}) is the largest rectangular subset of ${\mathcal{I}}$ for which the interpolation is computable (see discussion at the end of Sect. II-C) and, accordingly, ${\mathcal{I}}_{T}$ replaces ${\mathcal{I}}$ in the denominator of (5).

In the inverse approach $\bf W$ is mapped back to the original domain, so that its grid coincides with that of \hat{\mbox{\bf K}}^{\prime}. Then, the test statistic in this case is

[TABLE]

where, as above, the domain $\mathcal{I}_{T}$ of T_{\bm{\alpha}}^{-1}(\mbox{\bf W}) is the largest rectangular subset of ${\mathcal{I}}$ for which the interpolation is computable and ${\mathcal{I}}_{T}$ replaces ${\mathcal{I}}$ in the denominator of (6).

Since one is interested in finding the best possible match, [8] suggests using the following statistic

[TABLE]

When the parameter vector $\bm{\alpha}$ is not known, which is often the case in practice, it must be estimated. In [8] this is done by maximizing the test statistic in (12-13), which makes sense from a maximum likelihood point of view. Let $\mathcal{A}\subset{\mathbb{R}}^{n}$ be the set of feasible vectors $\bm{\alpha}$; then, the statistic used in the hypothesis test is

[TABLE]

For the case of scalar $\bm{\alpha}$ in (10) the inverse radial correction $T_{\bm{\alpha}}^{-1}(\mathbf{W)}$ needed in (13) can be approximated via the Lagrange Inversion Theorem [18, 3.6.6.] which yields

[TABLE]

Using the approach described above, the radial correction can be approximately inverted in many practical cases by finding the optimal value of $\alpha$ [8]. However, when more complex radial corrections as in (11) have been applied, a single parameter $\alpha$ may be not sufficient. To illustrate this fact, we consider the example of an image of size $3456\times 5184$ taken with a Canon 1200D camera using a Canon EF-S 10-18mm as lens and radially corrected with Adobe Lightroom (with settings for the mounted lens, using the strongest correction). We partitioned the image into non-overlapping annuli of width 64 pixels and found for each annulus—through exhaustive search—the value of $\alpha$ that maximizes $\text{PCE}_{\mathsf{inv}}(\alpha)$ in (13), where inversion is done via (16). The result is plotted in Fig. 1 as a function of the inner radius of the annulus. As it is quite apparent, there is a dependence of $\alpha$ with $r$ that indicates that one parameter alone is not sufficient to describe the radial transformation and that a more intricate relationship—even if parametric—must be sought.

III State of the Art

The PCE is very sensitive to the correct alignment of the locations corresponding to the estimated PRNU and the residual; this means that unless a value of $\bm{\alpha}$ very close to the true one is used in the mappings in (12) or (13), the resulting PCE will be very small, and hypothesis $H_{1}$ is likely to be rejected when it is in force. To illustrate this phenomenon, in Fig. 2 we show the function $\text{PCE}_{\mathsf{max}}(\alpha)$ for an image taken with a Panasonic DMC-ZS7 camera, shutter speed: 1/400 s, aperture: f4.4, focal length: 19.5 mm, and ISO 100. The stepsize in $\alpha$ is $2\cdot 10^{-3}$. As we can observe, under $H_{1}$ the function is very spiky, with the consequence that a sufficiently dense grid must be used; otherwise, it is easy to miss the peak. In addition, this spikiness precludes the use of gradient-based algorithms, because they would only work in the very close vicinity of the peak.

Therefore, any search grid in the parameter space has to be fine enough to be able to locate the maximum. The method in [8] considers that the transformations (both the direct and the inverse) are parameterized by a scalar $\alpha$ and starts by selecting a search interval $[-A,A]$ which is progressively made finer so that at each iteration $k$, with $k=1,\cdots,k_{\mathsf{max}}$, a grid with $2^{k}+1$ points is generated. Note that at the $k+1$-th iteration only $2^{k}$ new points are produced. A threshold $\tau_{1}$ is set so that if, after all $k_{\mathsf{max}}$ iterations, no $\alpha$ exists in the grid such that $\text{PCE}_{\mathsf{max}}(\alpha)>\tau_{1}$, then the search is stopped and a mismatch is declared (i.e., $H_{0}$ is decided). At every iteration, $\text{PCE}_{\mathsf{max}}$ is maximized over all grid points; this requires computing it only for the new points. Let $\alpha^{\circ}$ denote the grid point for which the maximum is obtained; if at some iteration $\text{PCE}_{\mathsf{max}}(\alpha^{\circ})>\tau_{1}$, then the search stops and the algorithm proceeds to the second stage in order to refine the value of $\alpha^{\circ}$. However, in order to speed up the process, the maximization skips the exhaustive enumeration of all grid points provided that $k>4$ whenever $\alpha^{\dagger}$ is found such that $\text{PCE}_{\mathsf{max}}(\alpha^{\dagger})>\tau_{2}$ (with $\tau_{2}>\tau_{1}$). In this case, the algorithm proceeds to the second stage by searching around $\alpha^{\dagger}$. The second stage takes the value of $\alpha$ with which the first stage was exited and constructs an interval with its two neighboring points in the grid. If $k^{*}$ is the exit value of $k$ for the first stage, then this interval has width $A/2^{k^{*}-1}$. Next, a golden section search is performed until the width of the interval is approximately $1/(8D_{2})$, with $D_{2}$ the half-diagonal of the image. Let $\alpha^{*}$ be the value found with the golden section search; then, if $\text{PCE}_{\mathsf{max}}(\alpha^{*})>\tau_{3}$ hypothesis $H_{1}$ is accepted, else, $H_{0}$ is declared. The thresholds suggested in [8] are $\tau_{1}=15$ and $\tau_{2}=\tau_{3}=75$, and $k_{\mathsf{max}}=7$. To reduce the computational load [8] downsamples the signals by a factor of two in each dimension; since this has an impact on accuracy in some cases, in the experimental section, we will consider both the downsampled (DS) and non-downsampled versions.

The method in [9] takes a different approach to perform the inversion of radially-corrected barrel distortions by employing the so-called linear patterns that are present in the residuals and are due to artifacts of the capturing device. These patterns are typically removed towards source attribution, but when kept, they serve as pilot signals that may be used to infer the radial correction distortion. The feature that is used to steer the parameter estimation is the energy of the linear pattern, defined for a given residual $\bf W$ as E(\mbox{\bf W})\doteq||\mbox{\bf c}||^{2}+||\mbox{\bf r}||^{2}, where $\bf c$ and $\bf r$ are vectors containing respectively the column and row averages of $\bf W$. Then, considering the set of fourth-order transformations $T_{\bm{\alpha}}(r)=r(1+\alpha_{2}r^{2}+\alpha_{4}r^{4}),$ where $\bm{\alpha}=(\alpha_{2},\alpha_{4})$, the method in [9] seeks to maximize $E(T_{\bm{\alpha}}^{-1}(\mathbf{W}))$ with respect to $\bm{\alpha}$, with the rationale that when the correct inverse transformation is applied, the linear pattern is recovered; otherwise, the column and row averages will be expected to produce low values. The fact that the transformation is now parameterized by two variables $\alpha_{2}$ and $\alpha_{4}$ gives more flexibility in inverting the transformation, but potentially incurs a larger computational cost. To make the optimization more manageable, a first stage consists in fitting a second-degree polynomial on variable $\alpha_{2}$ to values of $E(T_{\bm{\alpha}}^{-1}(\mathbf{W}))$ sampled on a grid for $\alpha_{2}\in[\alpha_{\mathsf{min}},\alpha_{\mathsf{max}}],\alpha_{\mathsf{min}}>0$, and $\alpha_{4}=0$. The reason for this choice of $\alpha_{4}$ is that in practice the contribution of $\alpha_{4}$ to $T_{\bm{\alpha}}(r)$ is only significant for large $r$, that is, far from the image center. This first stage yields the value $\alpha_{2}^{(1)}$ of $\alpha_{2}$ that maximizes the difference from the energy of the linear pattern and its polynomial fit. The second stage employs a Nelder-Mead optimization (using the linear pattern energy as cost function) that is initialized with three points derived from $\alpha_{2}^{(1)}$. This produces the two optimal radial correction parameters $(\alpha_{2}^{*},\alpha_{4}^{*})$. Due to noise, the previous procedure will yield an optimum $\alpha_{2}\neq 0$ regardless of whether radial correction was applied. Then, the decision is confirmed only if the cost function evaluated in a neighborhood of $(\alpha_{2}^{*},\alpha_{4}^{*})$ corroborates the existence of a significant peak; otherwise, the image is deemed to be not radially corrected.

Even though, as we will see in Sect. V, the performance of the two methods outlined above is rather good, they have two main intrinsic limitations that we aim at overcoming with our work: 1) their corresponding first stages employ an exhaustive search on a fixed grid. This fact, together with the high sensitivity of the PCE with respect to changes in the parameter vector $\bm{\alpha}$ about the correct one that results in a very spiky objective function, advise the use of a relatively tight grid to minimize the risk of missing the optimum. Unfortunately, this tightness entails a significant computational cost. 2) Again, due to the computational cost of an exhaustive search, the transformations $T_{\bm{\alpha}}$ and $T^{-1}_{\bm{\alpha}}$ use a small number of parameters: one in [8], and two in [9]. Therefore, these parameterization are unable to capture more complex radial corrections, such as those employed by editing programs, a trend that is likely to increase, as the capabilities of out-of-camera processing improve.

IV Proposed Method

In order to motivate the method proposed in this paper, we will rely on an example generated with the popular photo editing software Adobe Lightroom that will give us the necessary clues. Images were taken with a Canon 1200D camera and then radially corrected with Lightroom. In Fig. 3 we superimpose two $\text{PCE}_{\mathsf{inv}}$ maps (corresponding to $\alpha=-0.01$ and $\alpha=0.05$) in which $\text{PCE}_{\mathsf{inv}}(\alpha)$ is computed using (13) and (16) for non-overlapping blocks of size $64\times 64$. For mere illustrative purposes, and in order to enhance the visibility, the (radially corrected) image under analysis (from which $\bf W$ is computed) is one of the 20 flat-field images used to extract \hat{\mbox{\bf K}}^{\prime}. As we can see, the region where the PCE is significant is an annulus, and the position of the annulus depends on $\alpha$. This shows that if $L(r)$ denotes the radial correction induced by the software and $L^{-1}(r)$ its inverse, then for a given $\alpha=\alpha_{0}$, $T_{\alpha_{0}}^{-1}(r)\approx L^{-1}(r)$ only in a small neighborhood of some $r=r^{*}$.

This experiment clearly indicates that for complex radial corrections, an approach like (16) will not work. However, the fact that the inversion works locally suggests breaking the problem into non-overlapping concentric annuli as shown in Fig. 4, and solving each separately.

IV-A Set partitioning and transform computation

Let ${\mathcal{R}}_{k}$, $k=1,\cdots,L$, be the $k$th annulus described by an inner radius $r_{k}$ (recall that radii are scaled by $D_{2}$ so that $r=1$ corresponds to half of the image diagonal) and a width $\Delta_{k}$ as follows:

[TABLE]

The inner radii are generated as $r_{k+1}=r_{k}+\Delta_{k}$, with $r_{1}=0$, and the inner radius of the last annulus $r_{L}$ is such that $r_{L}<1<r_{L}+\Delta_{L}$ (see Fig. 4). This definition implies that the first annulus degenerates into a disk and the image is fully covered by annuli. Except for this degenerate annulus, in this work we will assume that $\Delta_{k}=\Delta$ for all $k$.

The experiment shown in Fig. 1 (obtained applying a brute force search for each annulus) suggests that a good modeling of the radial correction can be obtained by allowing $\alpha$ to vary with $r$, so (10) in this case becomes

[TABLE]

The idea is that by allowing $\alpha$ to be a function of $r$, we achieve much more flexibility in modeling complex distortions. Moreover, as long as the annuli are thin enough, the zero-th order approximation $\alpha(r)\approx\alpha(r_{k}+\Delta_{k}/2)\doteq\alpha_{k}$ will be reasonably good for all $r\in{\mathcal{R}}_{k}$. This local approximation will allow us to use (16) for the inverse transform. However, since we are allowing $\alpha$ to vary with $r$, instead of a locally cubic dependence, as in (18), it also makes sense to consider a locally linear one, i.e., $r^{\prime}=T_{\alpha(r)}(r)=r(1+\alpha(r))$. Even though for the generic mappings we will keep using $T_{\alpha_{k}}(r)$ and $T^{-1}_{\alpha_{k}}(r)$ for the sake of generality, we specialize them by adding the sub-indices $c$ to denote cubic, and $l$ to denote linear. Therefore, on each annulus we write

[TABLE]

whereas the corresponding inverse mappings are

[TABLE]

Note that the ranges of the inverse transforms in (IV-A) may be different because the image of each annulus will differ under the locally cubic and locally linear mappings.

Given a collection of annuli ${\mathcal{R}}_{k}$, $k=1,\cdots,L$, one can see the mapping $T_{\bm{\alpha}}(r)$ in (12) as a sequence of transformations $T_{\alpha_{k}}(r)$, $k=1,\cdots,L$, that is parameterized by a vector $\bm{\alpha}=[\alpha_{1},\cdots,\alpha_{L}]^{T}$. Obviously, the maximization of the PCE with respect to $\bm{\alpha}\in{\mathcal{A}}\doteq{\mathcal{A}}_{1}\times\cdots{\mathcal{A}}_{L}$, with ${\mathcal{A}_{k}}$ the feasible set for $\alpha_{k}$, would suffer from a combinatorial explosion due to the $L$ dimensions involved, so we will be interested in finding efficient alternative ways for performing an approximate maximization.

A first step is to treat each annulus separately and find the optimal value of $\alpha_{k}$ constrained to the $k$th annulus. There are several possible approaches at this stage. One would be to find $\alpha_{k}$ that maximizes the PCE constrained to the $k$th annulus; unfortunately, since the total PCE is not the sum of those constrained PCEs, it is quite difficult to work individually with each annulus using such a criterion. Instead, we have opted for a maximum likelihood estimation approach that aims at finding the $\alpha_{k}$ that has the highest likelihood of producing the observed cross-correlations with the estimated PRNU. Once we describe how the optimal $\alpha_{k}$ can be found for each annulus in an adaptive way (Sect. IV-C), we proceed by explaining how the PCE can be computed and updated (Sect. IV-D).

In the following, we give a formal description of the annuli for the inverse approach (i.e., using $T^{-1}_{\alpha_{k}}$) and afterwards indicate how to adapt the discussion to the direct approach. Let ${\mathcal{P}}_{k}$ be the set of points of the image grid that are contained in the $k$th annulus, i.e.,

[TABLE]

where multiplication of ${\mathcal{R}}_{k}$ by $D_{2}$ (i.e., half the diagonal in pixels) is necessary to re-scale the annulus back to integer-valued coordinates (recall that $r=1$ corresponds to half the diagonal).

Given \tilde{\mbox{\bf W}}=\mbox{\bf W}-\bar{\mbox{\bf W}} and ${\mathcal{P}}_{k}$, computation of T_{\alpha_{k}}^{-1}(\tilde{\mbox{\bf W}}) proceeds as follows (see Fig. 5). First, the image of the set ${\mathcal{P}}_{k}$ under $T^{-1}_{\alpha_{k}}$ , i.e. $T^{-1}_{\alpha_{k}}({\mathcal{P}}_{k})$ is calculated and the transformed points lying outside the image boundaries ${\mathcal{B}}$ are discarded, as the subsequent interpolation would not be computable. For the remaining points, T_{\alpha_{k}}^{-1}(\tilde{\mbox{\bf W}}) is obtained by interpolation from \tilde{\mbox{\bf W}}. We let ${\mathcal{Q}}_{k,\mathsf{inv}}(\alpha_{k})$ be the set of points of ${\mathcal{P}}_{k}$ for which their image under $T^{-1}_{\alpha_{k}}$ exists (the sub-index $\mathsf{inv}$ stands for ‘inverse approach’). Formally, this set is

[TABLE]

Notice that if the set ${\mathcal{P}}_{k}$ transformed via $T_{\alpha_{k}}^{-1}$ does not get out of the image bounds ${\mathcal{B}}$, then ${\mathcal{Q}}_{k,\mathsf{inv}}(\alpha_{k})={\mathcal{P}_{k}}$; otherwise, ${\mathcal{Q}}_{k,\mathsf{inv}}(\alpha_{k})\subset{\mathcal{P}_{k}}$. As a consequence, ${\mathcal{Q}}_{k,\mathsf{inv}}(\alpha_{k})\cap{\mathcal{P}}_{k}={\mathcal{Q}}_{k,\mathsf{inv}}(\alpha_{k})$. Also notice that, as explicitly indicated, the set ${\mathcal{Q}}_{k,\mathsf{inv}}(\alpha_{k})$ and, in particular, its cardinality, varies with $\alpha_{k}$.

For the direct approach, the considerations are similar. Basically, we have to exchange the roles of $T_{\alpha_{k}}$ and $T_{\alpha_{k}}^{-1}$. Recalling that the sub-index $\mathsf{dir}$ stands for ‘direct approach’, the set ${\mathcal{Q}}_{k,\mathsf{dir}}(\alpha_{k})$ can be formally written as

[TABLE]

IV-B Optimization with respect to $\alpha_{k}$

Once the annuli have been characterized, in this section we address the problem of finding the optimal values of $\alpha_{k}$ that parameterize the transformations $T_{\alpha_{k}}^{-1}$ and $T_{\alpha_{k}}$ for the $k$th annulus.

For the sake of compactness, we will find it useful to denote the cross-correlation and the energy of the transformed residual computed over ${\mathcal{Q}}_{k,\mathsf{inv}}(\alpha_{k})$ as, respectively,

[TABLE]

by making implicit the use of the inverse transformation $T^{-1}_{\alpha_{k}}(\cdot)$, and \hat{\mbox{\bf K}}^{\prime} and \tilde{\mbox{\bf W}}. Similarly, we denote by $\mathsf{\Phi}_{k,\mathsf{dir}}(\alpha_{k})$ and $\mathsf{E}_{k,\mathsf{dir}}(\alpha_{k})$ the cross-correlation and energy for the direct mapping $T_{\alpha_{k}}(\cdot)$ computed over ${\mathcal{Q}}_{k,\mathsf{dir}}(\alpha_{k})$.

In Appendix A we derive an estimator of $\alpha_{k}$ on the $k$th annulus. This estimator is rooted in the principle of maximum likelihood applied to the output of a bank of cross-correlations. For the inverse approach, this becomes

[TABLE]

where

[TABLE]

For the direct approach, the optimization is carried out after replacing the subindex $\mathsf{inv}$ by $\mathsf{dir}$ in both (26) and (27). We notice the proposed objective function is different from the PCE (constrained to the $k$th annulus); besides the theoretical justification in Appendix A, in [19] we provide empirical evidence that optimization of our objective function renders better global performance.

IV-C Adaptive optimization

One key observation from Fig. 1 is that the sequence $\alpha_{k}^{*}$, $k=1,\cdots,L$, changes smoothly for sufficiently small $\Delta_{k}$. This hints at the possibility of reducing the computational complexity of the exhaustive search by using an adaptive predictor. In our case, we will show experimentally that a linear predictor $\bf u$ with length $U$ suffices to achieve excellent results. In the following, we explain this adaptive procedure. As above, we will give the details for the inverse approach, as the direct one is methodologically identical.

First, we need to select an initial index that we will denote by $k_{0}$. To this end, we look for the annulus that gives the best results under no transformations (i.e., when $\alpha_{k}=0$). Formally, this implies that

[TABLE]

Once this initial point is found, the optimal value of $\alpha_{k_{0}}$ is found by exhaustive search in a discrete set around $\alpha_{k_{0}}=0$. Let ${\mathcal{A}}_{k_{0}}$ be such a neighborhood, then following (26), $\alpha_{k_{0}}^{*}=\arg\max_{\alpha_{k_{0}}\in{\mathcal{A}}_{k_{0}}}\varphi_{k_{0},\mathsf{inv}}(\alpha_{k_{0}})$.

We will find it useful to define an auxiliary sequence $\{\beta_{k}\}$ that is initialized as $\beta_{k}=\alpha^{*}_{k_{0}}\cdot\delta_{k-k_{0}}$, where $\delta_{k}$ is Kronecker’s delta.111Although from a notational point of view, it would be more correct to define a sequence for every iteration of the algorithm, we allow replacing values in this sequence in order to avoid overcomplicating the notation. This sequence is used to store the regressor values. Since the starting point is $k=k_{0}$, there are two possible directions for the prediction: forward (i.e., $k>k_{0}$), and backward (i.e, $k<k_{0}$).222Degenerate cases arise when $k_{0}=L$ or $k_{0}=1$, for which the forward and backward predictions, respectively, are not needed. We will describe how the former is carried out, and then indicate the modifications needed for the latter. We define the forward regressor at index $k$ as $\bm{\beta}_{k}^{T}\doteq[\beta_{k-U+1},\cdots,\beta_{k-1},\beta_{k}]$, where $U$ is the length. Notice that as a consequence of initializing the auxiliary sequence, $\bm{\beta}_{k_{0}}^{T}=[0,\cdots,0,\alpha^{*}_{k_{0}}]$. We also need a vector of weights at index $k$ that will be denoted by \mbox{\bf u}_{k}; this vector of length $U$ is initialized as \mbox{\bf u}_{k_{0}}^{T}=[0,\cdots,0,1]. Then, for $k>k_{0}$ the output of the predictor at index $k$ will be computed as

[TABLE]

for $k=k_{0}+1,\cdots,L$. This predicted value is refined by exhaustive search in a discrete neighborhood of $\hat{\alpha}_{k}$. Let ${\mathcal{A}}_{k}$ denote such a neighborhood; then $\alpha_{k}^{*}$ is obtained as in (26). The details on how the neighborhood ${\mathcal{A}}_{k}$ is constructed are given below. Before that, we explain the updating procedure for \mbox{\bf u}_{k} and $\bm{\beta}_{k}$. To that end, we define the a posteriori error at index $k$ as

[TABLE]

for $k=k_{0}+1,\cdots,L$. This error is used to drive the adaptive algorithm. It is easy to show that the gradient vector of $|e_{k}|^{2}$ with respect to the weights vector \mbox{\bf u}_{k-1} is equal to $-2e_{k}\bm{\beta}_{k-1}$. Then, following the Least Mean Squares algorithm [20], we propose to update the weights by taking a step in the direction of the negative gradient, that is,

[TABLE]

where $\mu$ is the so-called step-size. The update of the sequence $\{\beta_{k}\}$ containing the regressor is done by making $\beta_{k}=\alpha_{k}^{*}$; the forward regressor vector $\bm{\beta}_{k}$ is updated accordingly. This iterative procedure is then repeated by going back to (29) and proceeding until the sequence $\alpha_{k_{0}+1}^{*},\alpha_{k_{0}+2}^{*},\cdots,\alpha_{L}^{*}$ is produced.

The backward prediction proceeds in a similar way, but now vector $\bm{\beta}_{k}$ is defined as $\bm{\beta}_{k}^{T}\doteq[\beta_{k},\beta_{k+1},\cdots,\beta_{k+U-1}]$; this means that at the backward initialization, vector $\bm{\beta}_{k_{0}}$ will take advantage of the availability of values of $\alpha_{k}^{*}$ that have been already computed, i.e., $\bm{\beta}_{k_{0}}=[\alpha^{*}_{k_{0}},\alpha^{*}_{k_{0}+1},\cdots,\alpha^{*}_{k_{0}+U-1}]^{T}$. The weights vector for the backward prediction \mbox{\bf u}_{k_{0}} is initialized as \mbox{\bf u}_{k_{0}}=[1,0,\cdots,0]. Now this weights vector is updated in the reverse direction:

[TABLE]

and again the sequence $\{\beta_{k}\}$ containing the regressor is updated by making $\beta_{k}=\alpha_{k}^{*}$; the backward regressor vector $\bm{\beta}_{k}$ is updated accordingly. The algorithm thus generates the sequence $\alpha_{k_{0}-1}^{*},\alpha_{k_{0}-2}^{*},\cdots,\alpha_{1}^{*}$.

After both forward and backward predictions are finished, the optimal vector is $\bm{\alpha}^{*}_{\mathsf{inv}}=[\alpha_{1}^{*},\cdots,\alpha_{L}^{*}]^{T}\in{\mathcal{A}}$, where once again we have added the subindex $\mathsf{inv}$ to stress the fact that we are dealing with the inverse approach. The same procedure applied to the direct approach will yield an optimal vector $\bm{\alpha}^{*}_{\mathsf{dir}}$. The pseudo-code for the proposed algorithm is provided in the technical report [19].333The code is available at https://github.com/AMontiB/AdaptivePRNUCameraAttribution

One critical point of the algorithm is the refining of $\hat{\alpha}_{k}$ that produces $\alpha_{k}^{*}$. While smarter strategies might be possible, here we perform an exhaustive search around $\hat{\alpha}_{k}$ in a discrete set ${\mathcal{A}}_{k}$. Of course, the cardinality of this set must be kept at a small value in order to limit the computational burden. On the other hand, the discrete points must be generated finely enough to output a value that is sufficiently close to the optimal. We thus employ two parameters to describe the set: $\lambda_{k}$ that controls the resolution, and $A_{k}$ that is an odd integer that determines the number of points. Then, given $\hat{\alpha}_{k}$ and these parameters, the search set is constructed as:

[TABLE]

Note that this construction guarantees that $|{\mathcal{A}}_{k}|=A_{k}$. The parameter $\lambda_{k}$ is selected to be commensurate with $|\alpha^{*}_{k}-\alpha^{*}_{k-1}|$ in the forward case (resp. $|\alpha^{*}_{k}-\alpha^{*}_{k+1}|$ in the backward case), so that the smaller the change in $\alpha_{k}^{*}$, the finer the grid. In Sect. IV-G we give more details about the rules that were employed to generate $\lambda_{k}$ for the experiments. Regarding the size of the set $A_{k}$, this is updated in the same loop as the predictor; for the forward predictor, the rule is as follows: if for index $k$ the maximum $\alpha_{k}^{*}$ is found at one of the extremes of the set ${\mathcal{A}}_{k}$ (i.e., $\alpha^{*}_{k}=\hat{\alpha}_{k}-\lambda_{k}\cdot(A_{k}-1)/2$ or $\alpha^{*}_{k}=\hat{\alpha}_{k}+\lambda_{k}\cdot(A_{k}-1)/2$) then the size of the set is increased at the following iteration, i.e., $A_{k+1}=A_{k}+2$. Otherwise, if $A_{k}$ is already small, i.e., $A_{k}=A_{\text{min}}$ for some minimum size $A_{\text{min}}$, then $A_{k+1}=A_{\text{min}}$; else (i.e, if the maximum in ${\mathcal{A}}_{k}$ is not found at either of the extremes, and the set is large enough), the size is decreased at the following iteration, i.e., $A_{k+1}=A_{k}-2$. This update is intended to find a compromise between the size of the set and the objective of capturing the optimal $\alpha_{k}$. For the backward prediction the reasoning is identical, but updating $A_{k-1}$ from $A_{k}$ (see [19] for an example of the evolution of $A_{k}$).

IV-D PCE computation for the optimal $\bm{\alpha}$

As a result of the adaptive algorithm presented in the previous section, it is possible to compute the PCEs that are required in the hypothesis test, that is, $\text{PCE}_{\mathsf{inv}}(\bm{\alpha}^{*}_{\mathsf{inv}})$ and $\text{PCE}_{\mathsf{dir}}(\bm{\alpha}^{*}_{\mathsf{dir}})$, see the definitions in (12) and (13). In both cases, the numerator and denominator of the PCE are already available, as they are required for the optimization. The only additional computations are simple sums to accumulate the results corresponding to the different annuli. To see how this is so for the inverse approach, notice first that the right hand side of (13) requires computing the difference T^{-1}_{\bm{\alpha}^{*}}(\mbox{\bf W})-\overline{T^{-1}_{\bm{\alpha}^{*}}(\mbox{\bf W})} (cf. the expression of the PCE in (6)), which can be simplified by noticing that: 1) It is reasonable to write \overline{T^{-1}_{\bm{\alpha}^{*}}(\mbox{\bf W})}\approx T^{-1}_{\bm{\alpha}^{*}}(\bar{\mbox{\bf W}}) because $T^{-1}_{\bm{\alpha}^{*}}$ is a geometrical transformation that will not substantially alter the mean value of the residual.444Strict equality does not hold because \bar{\mbox{\bf W}} and T^{-1}_{\bm{\alpha}^{*}}(\bar{\mbox{\bf W}}) do not have the same support. 2) Due to zero-meaning on the residual, it is possible to write \bar{\mbox{\bf W}}=\mbox{\bf 0}. With these considerations, we can write T^{-1}_{\bm{\alpha}^{*}}(\mbox{\bf W})-\overline{T^{-1}_{\bm{\alpha}^{*}}(\mbox{\bf W})}\approx T^{-1}_{\bm{\alpha}^{*}}(\tilde{\mbox{\bf W}}), which is simpler to compute.

With this approximation, the numerator of \text{PCE}(\hat{\mbox{\bf K}}^{\prime},T^{-1}_{\bm{\alpha}^{*}}(\mbox{\bf W})) can be expanded as follows

[TABLE]

Now we can easily identify each of the $L$ summands in (34) as the numerator of $\varphi_{k,\mathsf{inv}}(\alpha^{*}_{k})$ in (27) which can be stored during the adaptive optimization process.

The denominator of the PCE requires more attention. With the approximation above, this denominator is \frac{1}{|{\mathcal{I}_{T}}\backslash{\mathcal{S}}|}\sum_{\mbox{\bf s}\in{\mathcal{I}_{T}}\backslash{\mathcal{S}}}\langle\hat{\mbox{\bf K}}^{\prime},C(T^{-1}_{\bm{\alpha}^{*}}(\tilde{\mbox{\bf W}}),\mbox{\bf s})\rangle^{2} which is nothing but a sample estimate of the variance of the cross-correlation of \hat{\mbox{\bf K}}^{\prime} and T^{-1}_{\bm{\alpha}^{*}}(\tilde{\mbox{\bf W}}). In [19, Sect. VIII] we derive and discuss a simpler sample estimate that is more statistically efficient (i.e., has a lower variance). This fully justifies the approximation

[TABLE]

where \hat{\sigma}_{\hat{K}^{\prime}}^{2}\doteq||\hat{\mbox{\bf K}}^{\prime}||^{2}/|{\mathcal{I}}| (recall that $\hat{K}^{\prime}_{i,j}$ exists for all $(i,j)\in{\mathcal{I}}$), and $\kappa$ is a factor that takes into account the fact that the cardinalities of ${\mathcal{I}}$ and $\bigcup_{k=1}^{L}{\mathcal{Q}}_{k,\mathsf{inv}}(\alpha^{*}_{k})$ are different. (In practice, $\kappa$ will be close to 1, so it can be dropped.) Once again, the $L$ summands in (35) are already available as the denominator of (27).

IV-E Early stopping

The partition into annuli offers one remarkable byproduct: taking inspiration from [21], it is possible to stop processing annuli (and declare that $H_{1}$ holds) if a cumulative PCE exceeds a predefined threshold. Following the approximations in the previous subsection, one might be tempted to compute a cumulative PCE by using the numerators and denominators already produced during the optimization. In this way, the optimization would not need to be carried out for all annuli but instead it could be stopped as soon as the PCE computed so far exceeds the threshold. Unfortunately, this approach would be incorrect, because while a fraction with sums in the numerator can be expanded into a sum of fractions, this is not the case when there are sums in the denominator. Therefore, if we want to implement an early stopping mechanism, we need to seek ways to further approximate the denominator of the PCE without actually computing all the elements of $\bm{\alpha}^{*}$. To this end, we can ask ourselves how sensitive is the right hand side of (35) to changes in $\alpha_{k}^{*}$; after all, since each of the $L$ summands is an estimate of the variance of the transformed residual inside an annulus, one would expect not much variation for realistic values of $\alpha$. If this were the case, then one might approximate the right hand side of (35) (which corresponds to the optimal vector $\bm{\alpha}_{\mathsf{inv}}^{*}$) by computing it for any reasonable value of $\bm{\alpha}_{\mathsf{inv}}$ without involving any optimization.

In order to illustrate the feasibility of this approximation, we show in Fig. 6 the values of the sample variance of a transformed residual computed in each annulus, i.e., $\frac{\mathsf{E}_{k,\mathsf{inv}}(\alpha_{k})}{|{\mathcal{Q}}_{k,\mathsf{inv}}(\alpha_{k})|}$ as a function of $\alpha_{k}$ for several annuli (i.e., $k=18,22,33$) and for cubic inverse mappings, see (IV-A). Fig. 6 also shows the value of the variance estimated from the full-size transformed residual, i.e., \frac{1}{|{\mathcal{I}}_{T}|}||T^{-1}_{\alpha_{k}}(\tilde{\mbox{\bf W}})||^{2}. Bi-linear interpolation is used in all cases.

As we can see, the variance estimate is fairly constant for different values of $\alpha_{k}$, except in a neighborhood of zero. Moreover, this is similar to the variance estimate obtained from the whole transformed residual, so the latter can be used in place of the variance estimate for a specific annulus. The reason for the spike at $\alpha_{k}=0$ is that the interpolation that is needed for computing the inverse mapping when $\alpha_{k}\neq 0$ produces a reduction in the variance of the transformed residual. This reduction depends on the square magnitude of the interpolation filter at different sampling points. In general, the grids before and after the interpolation are not related through rational numbers, but for certain rings and values of $\alpha$, moiré patterns between the sampling grids may appear; this is why in Fig. 6 a ripple near zero is observed for the rings $k=18,22$. The energy reduction phenomenon has been reported in [22] in a different scenario but related to ours.

The invariance discussed in the previous paragraph suggests several ways of approximating the right hand side of (35); for instance, it is possible to pick any value of $\alpha$, say $\alpha_{\mathsf{f}}$, sufficiently far from $\alpha=0$ and for all the annuli use the same transformation $T_{\alpha_{\mathsf{f}}}^{-1}$ in place of $T_{\alpha^{*}_{k}}^{-1}$. We remark that the reason why the neighborhood of $\alpha=0$ should be excluded when selecting $\alpha_{\mathsf{f}}$ is the fact that inside such a neighborhood the denominator of the PCE is overestimated and, consequently, the PCE underestimated.

Another way of approximating the right hand side of (35) which offers a slightly better performance than the former is to use the values of $\alpha_{k}^{*}$ already available from the optimization to update the approximation. This comes at practically no cost because the corresponding term $\mathsf{E}_{k,\mathsf{inv}}(\alpha_{k}^{*})$ needs to be computed anyway during the optimization. For those annuli whose $\alpha_{k}^{*}$ is not available yet, the corresponding term is substituted by its approximation computed at $\alpha_{k}=\alpha_{\mathsf{f}}$.

We explain next how to compute the Cumulative PCE at the $n$th iteration which we will denote by \text{CPCE}_{n,\mathsf{inv}}(\hat{\mbox{\bf K}}^{\prime},\mbox{\bf W}). First, we need a mapping $\xi:\{1,\cdots,L\}\to\{1,\cdots,L\}$, from the natural order to the one induced by the proposed iterative procedure, i.e., $\xi(1)\mapsto k_{0},\xi(2)\mapsto k_{0}+1,\cdots,\xi(L-k_{0})\mapsto L,\xi(L-k_{0}+1)\mapsto k_{0}-1,\cdots,\xi(L)\mapsto 1$. Then,

[TABLE]

Thus, the early-stopping algorithm will declare a match and stop if for some $n=1,\cdots,L$, \text{CPCE}_{n,\mathsf{inv}}(\hat{\mbox{\bf K}}^{\prime},\mbox{\bf W})>\tau_{c} is satisfied. The value of $\tau_{c}$ is set experimentally to achieve the desired False Positive Rate (FPR).

Given the numerator and denominator of \text{CPCE}_{n,\mathsf{inv}}(\hat{\mbox{\bf K}}^{\prime},\mbox{\bf W}), and once $\alpha_{\xi(n+1)}^{*}$ is available, the numerator of \text{CPCE}_{n+1,\mathsf{inv}}(\hat{\mbox{\bf K}}^{\prime},\mbox{\bf W}) is updated by adding $\mathsf{\Phi}_{\xi(n+1),\mathsf{inv}}(\alpha^{*}_{\xi(n+1)})$, while the update of the denominator requires adding $\mathsf{E}_{\xi(n+1),\mathsf{inv}}(\alpha^{*}_{\xi(n+1)})$ and subtracting $\mathsf{E}_{\xi(n+1),\mathsf{inv}}(\alpha_{\mathsf{f}})$.

A similar definition follows for the Cumulative PCE in the direct approach \text{CPCE}_{n,\mathsf{dir}}(\hat{\mbox{\bf K}}^{\prime},\mbox{\bf W}) and the corresponding early stopping criterion.

IV-F Parameter inheritance

As we have discussed, the test decision statistic takes the maximum of the PCEs computed through the direct and the inverse approaches. This implies that it is necessary to compute the optimal vector $\bm{\alpha}^{*}$ for both approaches, so the computational complexity is roughly doubled. This also holds if the early stopping criterion introduced above is applied. In such a case, the iterations for both the direct and the inverse approaches are made in parallel, so that for every $k$ both \text{CPCE}_{n,\mathsf{dir}}(\hat{\mbox{\bf K}}^{\prime},\mbox{\bf W}) and \text{CPCE}_{n,\mathsf{inv}}(\hat{\mbox{\bf K}}^{\prime},\mbox{\bf W}) are checked against the threshold in order to stop as soon as possible.

There is one sub-optimal way to alleviate the computational burden due to keeping the two approaches. We term it parameter inheritance and basically consists in using for the direct approach the same vector $\bm{\alpha}^{*}$ that was computed for the inverse approach. Of course, the latter is not necessarily optimal for the direct approach, but the rationale is that inside each annulus ${\mathcal{R}}_{k}$ the direct and inverse transformations nearly correspond to each other for the same value of $\alpha_{k}$. Perfect correspondence does not exist because the inverse transformation is only an approximation and due to the fact that the search algorithm is prone to errors due to noise and insufficient resolution.

IV-G Parameter default values

In this section we provide the default values for the parameters of our algorithm and discuss some decisions regarding the initialization. These default values were used in the experiments reported in Sect. V. Specifically, the radius of the inner disk $r_{1}$ is such that $r_{1}\cdot D_{2}$ equals 250 pixels and the width of each annulus $\Delta_{k}$ is such that $\Delta_{k}\cdot D_{2}$ equals 64 pixels. Both values are chosen as a compromise between performance and computational cost. For the linear predictor we set $U=6$, $\mu=1$ and $A_{\mathsf{min}}=7$.

The initial search set $\mathcal{A}_{k_{0}}$ is given by $\mathcal{A}_{k_{0}}=\{-0.22,-0.21,\cdots,0.21,0.22\}$, which is the same range as used and justified in [8] to cover a variety of barrel and pincushion distortions. However, in our case we apply a coarser resolution for computational reasons and because the adaptive nature of our algorithm automatically adjusts to finer resolutions after a few iterations. We are aware that in [9] a wider range was preferred (even if just to invert pincushion distortions), so we carried out some experiments with images taken with the Canon 1200D camera and radially corrected with Adobe Lightroom using the lens distortion model of a different device (see Section V), since this combination produces some of the strongest and most variable radial corrections of our dataset. In these experiments, the search set was expanded to $\mathcal{A}_{k_{0}}=\{-0.50,-0,49,\cdots,0.49,0.50\}$. While it is true that this set allows in some cases to get closer to the proper $\alpha_{k_{0}}$, we found no significant differences in terms of performance with the previous initialization; as mentioned, this is due to our algorithm quickly finding the right range for $\alpha_{k}$ after few iterations. In contrast, the computational load of using the enlarged search set would be larger; for this reason, we recommend $\mathcal{A}_{k_{0}}=\{-0.22,-0.21,\cdots,0.21,0.22\}$. For an in-depth complementary discussion on the initial set, please see [19].

After the initial search, for the forward prediction ${\mathcal{A}}_{k_{0}+1}$ is given by (33) with $\lambda_{k_{0}+1}=0.001$ and $A_{k_{0}+1}=9$. For the following iterations,

[TABLE]

Identical considerations to the previous paragraph are made in regard to the backward prediction, where $k_{0}+1$ is replaced now by $k_{0}-1$ and in (37) $k-1$ is replaced by $k+1$.

V Experimental Results

In order to measure the performance of the methods presented in Sect. IV and compare them with the state of the art in [8] and [9], we built a test dataset composed of 3645 images, of which 2037 were taken with the following compact cameras and radially corrected “in-camera” (i.e., by the camera software): Canon SX230 HS (188 images), Panasonic ZS7 (170 images), Canon SX40 (57 images), Canon SX210 (82 images), and Nikon S9100 (1540 images). All these images were downloaded from Flickr, as done in [8] and [9]; for this reason, there is an uneven distribution of images per device. 1508 of the remaining images in the test dataset were taken with the Canon 1200D (a reflex camera not applying any type of in-camera post-processing) with the following Canon Zoom Lenses: 1) EF-S 10-18 mm 1:4-5.6 IS STM; 2) EF-S 18-55 mm 1:3.5-5.6; 3) EF 75-300 mm 1:4-5.6, all radially corrected “out-camera” with third-party editing software: Adobe Lightroom Classic CC 2017, Adobe Photoshop CC 2017, PT Lens v2.0 (Macbook) and Gimp 2.10.14. Specifically, 377 images were corrected with each of these tools. With Adobe Lightroom we applied the correction model specific to the lens used to take the picture, thanks to the database of radial correction models Lightroom is equipped with. For the other editing software we applied the strongest radial correction available, as those tools cannot be tuned to a specific lens model. The last 100 images in the test dataset were also taken with the Canon 1200D camera but corrected with Lightroom using models for other lenses (i.e., Nikon, Tamron, Apple, Huawei and DJI, with 20 images each), always applying the strongest radial correction. This latter subset will be labeled as “Lightroom*” in the following.

Images in the test dataset were JPEG compressed with a QFs in the range 90-98. For each device, the same QF is consistently used; see [19] for details. The reference PRNUs for carrying out the tests were estimated for each device using (4) with $L=20$ natural images (not used for testing) compressed with matching QFs to the test subset of that device. For the compact devices, since the in-camera corrections depend on the focal length, fixed specific values of the latter were sought in order to estimate the respective PRNUs; whenever enough images were available for a certain device and focal length, a different fingerprint was estimated and the results averaged for each device. In all cases, hypothesis $H_{1}$ was tested with images taken with focal lengths different from those used to estimate the fingerprints. We refer the reader to [19] for full details. When, under hypothesis $H_{0}$, the test images and the fingerprints have different sizes, we crop the central part of the largest to match its size to the smallest [7].

Next, we describe the identifiers used to refer to the different variants of our method in the figures and tables in this section. With “Dir” and “Inv” we indicate those cases where $\text{CPCE}_{n,\text{dir}}(\hat{\mathbf{K}^{\prime}},\mathbf{W})$ and $\text{CPCE}_{n,\text{inv}}(\hat{\mathbf{K}^{\prime}},\mathbf{W})$ are respectively used as the only test statistics. By “2W” we refer to the “two-way” case in which both the direct and the inverse approaches are used and $H_{1}$ is decided if either $\text{CPCE}_{n,\text{dir}}(\hat{\mathbf{K}^{\prime}},\mathbf{W})$ or $\text{CPCE}_{n,\text{inv}}(\hat{\mathbf{K}^{\prime}},\mathbf{W})$ are above the threshold for any $n\in\{1,\cdots,L\}$. To alleviate the computational load of the “two-way” parameter optimization, recall that in Sect.IV-F we proposed to inherit the parameters of one approach to the other. We will use the label $\overrightarrow{\text{DI}}$ to indicate inheritance of $\alpha_{n}^{*}$ from the direct approach to the inverse one; and $\overrightarrow{\text{ID}}$ vice versa. On the other hand, with the labels “Cub” and “Lin” we refer to the cubic and the linear radial correction models, respectively; see (IV-A). In all reported cases, the early stopping strategy from Sect. IV-E is imposed.

All the tests were run on a server with the following characteristics: 16 Cores, Processors 2xXeon E5-2667v3 3.2 GHz and RAM 192 GB; our implementation requires at most 5GB of RAM. In experimentally comparing the variants of our method with the algorithms proposed in [8] and [9], we noticed that [8] was tested on images of size $3000\times 4000$ that are, on average, larger than the in-camera corrected images in our dataset (refer to Table II for the image sizes in each subset). This explains the slightly worse performance measured here (with downsampling) compared to that reported in [8].

In Table I we provide the fixed thresholds $\tau_{c}$ (measured over the entire test dataset) that ensure False Positive Rates (FPR) of $0.05$ and $0.01$ together with the corresponding True Positive Rates (TPR) for the different variants of our method and those in [9] and [8] (with and without DS).

A breakdown by subset of the previous table is given in Table II which also shows the time consumed to declare a match (under $H_{1}$) by the different alternatives (with early stopping in those cases where it applies).555The time consumed when no early stopping is in force is given in the technical report [19]. For reasons of space, we have excluded the method in [9] which yields a modest performance, as well as the worst-performing variants of our method (cf. Table I); see [19] for fully comprehensive results. The Receiver Operating Characteristic (ROC) curves for the variants in Table II are plotted in Fig. 7, where we have also added for comparison the baseline (BL) obtained by using $\text{PCE}(\hat{\mathbf{K}^{\prime}},\mathbf{W})$ (i.e., with no transformations of either the PRNU or the residual) as test statistic.

From the results in Fig. 7 and Table I, it is possible to conclude that the best performing variants of our method correspond to the cubic correction model (“Cub”), with “$\overrightarrow{\text{DI}}$”, “Dir”, “Inv” and “2W” all achieving similar TPRs for the target FPRs. On the other hand, the average execution time of the “one way” variants, i.e. “Dir” and “Inv”, is lower because only one statistic has to be computed per iteration. The experiments, conducted on a dataset composed of a variety of radial corrections, show that our variants outperform [8] (both with or without DS) in terms of TPR. Moreover, thanks to our early stopping strategy, our fastest versions (i.e., “Dir, Cub” and “Inv, Cub”) achieve under $H_{1}$ execution times that are comparable to [8] with DS. We also note that the original solution proposed in [8] (i.e. with DS) achieves a limited performance both on low-resolution devices (i.e. SX230 and ZS7) and in presence of complex out-camera radial corrections as those applied by Adobe Lightroom. Adapting [8] to avoid DS results in a significant performance increase in those difficult cases, at the expense of a much more costly execution. Nevertheless, for some severe radial corrections like those in our “Lightroom*” subset, using the full resolution in [8] is still not sufficient. In contrast, our method is able to adapt to this high complexity and offers an excellent performance with an affordable execution time.

VI Conclusions

In this paper, we have proposed an adaptive method for PRNU-based camera attribution that is able to cope with complex radial distortion corrections, as those performed in-camera by most compact models and out-camera by image processing software. Existing approaches try to either “correct” the reference fingerprint or invert the correction by applying a further geometric transformation that, in order to avoid a combinatorial explosion, must use a reduced number of parameters. In turn, this limitation accounts for unsatisfactory performance when complex radial distortion corrections are in effect, an undesirable aspect in view of the trend of more elaborate transformations that are made possible by ever more powerful distortion correction firmware/software. Our approach is radically different: by applying a divide-and-conquer principle, embodied in the use of annuli, we are able to: 1) allow for complex distortion corrections, as locally the transformation undergone by each annulus is much simpler; 2) implement an early stopping strategy that offers large computational savings. The results presented in the paper clearly reveal that our algorithm (in most of its variants) outperforms the state of the art when accuracy and computational load are considered.

We believe that the adaptive approach proposed here could also be fruitful in other very challenging camera attribution scenarios with a number of latent parameters, such as in HDR images [7], in-camera-stabilized videos [23], and emerging in-camera processing [24].

Appendix A Derivation of the estimator of $\alpha_{k}^{*}$

In this Appendix we derive a plausible estimator of $\alpha_{k}^{*}$ under the inverse approach; the derivation would be identical for the direct approach and, hence, is skipped here. See the definition of $\alpha_{k}^{*}$ in (26). We introduce a super-index in $\alpha_{k}$ to enumerate the elements of the candidate set ${\mathcal{A}}_{k}$, i.e., $\{\alpha_{k}^{(n)}:n=1,\cdots,A_{k}\}={\mathcal{A}_{k}}$.

We assume that $H_{1}$ holds, i.e., $\bf I$ contains \mbox{\bf K}^{\prime}, and the following model for the residuals:

[TABLE]

for all $(i,j)\in{\mathcal{Q}}_{k,\mathsf{inv}}(\alpha_{k}^{(n)})$ and $n\in\{1,\cdots,A_{k}\}$. In (38) $\alpha^{\dagger}_{k}$ represents the true (locally for the $k$th annulus) value of $\alpha$. The multipliers $\gamma_{i,j}^{(n)}$ are non-negative and take into account both the multiplicative effect of the image $\bf I$ and the gain of the effective denoising filter (which also impacts on the estimate \hat{\mbox{\bf K}}^{\prime} of the true PRNU). We argue that these multipliers are very hard to estimate accurately; as a consequence, a full maximum likelihood decision will not be possible and some simplifications will be required. One such simplification is to consider that the cross-correlations between \hat{\mbox{\bf K}}^{\prime} and T^{-1}_{\alpha_{k}^{(n)}}(\tilde{\mbox{\bf W}}), for all $n=1,\cdots,A_{k}$, constitute a set of sufficient statistics for the estimation problem. Recall from (24) that these cross-correlations are denoted by $\mathsf{\Phi}_{k,\mathsf{inv}}(\alpha^{(n)}_{k})$.

We make the following hypotheses:

1. Spikiness: The $\alpha_{k}^{(n)}$ are sufficiently separated so that the $\mathsf{\Phi}_{k,\mathsf{inv}}(\alpha^{(n)}_{k})$ are mutually uncorrelated and $\mathbb{E}\{\mathsf{\Phi}_{k,\mathsf{inv}}(\alpha^{(n)}_{k})\}=0$, for all $n=1,\cdots,A_{k}$, except for $n=l$, where $l$ is such that $\alpha^{(l)}_{k}$ is the closest to the true value $\alpha^{\dagger}_{k}$ and the expectation is taken over the underlying distribution of \mbox{\bf K}^{\prime}. This hypothesis is reasonable in view of the spikiness of the PCE with $\alpha$ (see Fig. 2). We also assume that $\alpha^{(l)}_{k}$ is close enough to $\alpha^{\dagger}_{k}$ so that [T^{-1}_{\alpha^{(l)}_{k}}(T_{\alpha^{\dagger}_{k}}(\hat{\mbox{\bf K}}^{\prime}))]_{i,j}\approx\hat{K}^{\prime}_{i,j} for all $(i,j)\in{\mathcal{Q}}_{k,\mathsf{inv}}(\alpha_{k}^{(l)})$.

2. Uncorrelatedness: In (38), $N_{i,j}^{(n)}$ and \gamma_{i,j}^{(n)}[T^{-1}_{\alpha^{(l)}_{k}}(\hat{\mbox{\bf K}}^{\prime})]_{i,j} are zero-mean and mutually uncorrelated for all $(i,j)\in{\mathcal{Q}}_{k,\mathsf{inv}}(\alpha_{k}^{(n)})$ and $n\in\{1,\cdots,A_{k}\}$. For any $l,n\in\{1,\cdots,A_{k}\}$, $l\neq n$, the variables $\hat{K}^{\prime}_{i,j}\cdot N_{i,j}^{(n)}$ and $\hat{K}^{\prime}_{u,v}\cdot N_{u,v}^{(l)}$ are mutually uncorrelated for every $(i,j)\in{\mathcal{Q}}_{k,\mathsf{inv}}(\alpha_{k}^{(n)})$ and every $(u,v)\in{\mathcal{Q}}_{k,\mathsf{inv}}(\alpha_{k}^{(l)})$.

3. Weak PRNU: In (38), \left|\gamma_{i,j}^{(n)}[T^{-1}_{\alpha^{(n)}_{k}}(T_{\alpha^{\dagger}_{k}}(\hat{\mbox{\bf K}}^{\prime}))]_{i,j}\right|\ll|N_{i,j}| for a large number of pixels of each annulus; we write this more precisely as

[TABLE]

for all $n\in\{1,\cdots,A_{k}\}$.

As a consequence of the spikiness and uncorrelatedness assumptions above and the Central Limit Theorem (which is applicable if we assume that $|{\mathcal{Q}}_{k,\mathsf{inv}}(\alpha_{k}^{(n)})|$ is large for all $n\in\{1,\cdots,A_{k}\}$), the variables $\mathsf{\Phi}_{k,\mathsf{inv}}(\alpha^{(n)}_{k})$ will be well modeled by independent Gaussian distributions, so the cross-correlations will be

[TABLE]

where ${\mathcal{N}}(\mu,\sigma^{2})$ denotes a Gaussian with mean $\mu$ and variance $\sigma^{2}$, $\mu^{(l)}$ denotes the expected value of the cross-correlation for the value of $\alpha^{(n)}_{k}\in{\mathcal{A}}_{k}$ that is closest to $\alpha^{\dagger}_{k}$, and $(\sigma^{(n)})^{2}$, $n=1,\cdots,A_{k}$, denote the variances of the cross-correlations.

For all $n=1,\cdots,A_{k}$, the variances $(\sigma^{(n)})^{2}$ can be written as $\text{Var}\{\sum_{(i,j)\in{\mathcal{Q}}_{k,\mathsf{inv}}(\alpha_{k})}\hat{K}^{\prime}_{i,j}\cdot N_{i,j}\}\approx\hat{\sigma}_{\hat{K}^{\prime}}^{2}\sum_{(i,j)\in{\mathcal{Q}}_{k,\mathsf{inv}}(\alpha_{k})}N^{2}_{i,j}$. As a consequence of the weak PRNU assumption $\sum_{(i,j)\in{\mathcal{Q}}_{k,\mathsf{inv}}(\alpha_{k})}N^{2}_{i,j}\approx\mathsf{E}_{k,\mathsf{inv}}(\alpha_{k}^{(n)})$. Therefore, $(\sigma^{(n)})^{2}\approx\hat{\sigma}_{\hat{K}^{\prime}}^{2}\mathsf{E}_{k,\mathsf{inv}}(\alpha_{k}^{(n)})$.

Let $f_{\mathcal{N}}(Y;\mu,\sigma)$ denote the Gaussian pdf on random variable $Y\sim{\mathcal{N}}(\mu,\sigma^{2})$. Also, let $\mathcal{E}_{l}$ denote the event “$\alpha^{(l)}_{k}$, $l\in\{1,\cdots,A_{k}\}$ is the closest to the true value $\alpha^{\dagger}_{k}$”. Then the likelihood of jointly observing the cross-correlations $\mathsf{\Phi}_{k,\mathsf{inv}}(\alpha^{(n)}_{k})$ conditioned on ${\mathcal{E}}_{l}$ is

[TABLE]

The maximum likelihood estimator would be obtained by maximizing the likelihood in (A) with respect to $l$. The estimator will not change if we divide (A) by $\prod_{n=1}^{A_{k}}f_{\mathcal{N}}(\mathsf{\Phi}_{k,\mathsf{inv}}(\alpha^{(n)}_{k});0,\sigma^{(n)})$; this gives the following simpler likelihood function

[TABLE]

Taking the logarithm and simplifying, we find that the maximum likelihood estimator is equivalent to solving

[TABLE]

where

[TABLE]

and making $\alpha^{*}_{k}=\alpha_{k}^{(l^{*})}$.

Notice that, as discussed above, $(\sigma^{(l)})^{2}$ can be replaced by its estimator $\hat{\sigma}_{\hat{K}^{\prime}}^{2}\mathsf{E}_{k,\mathsf{inv}}(\alpha_{k}^{(l)})$ in (45). Unfortunately, producing a reliable estimator of $\mu^{(l)}$ is not feasible due to the unavailability of the gains $\gamma_{i,j}^{(l)}$. For this reason, we turn our attention to suboptimal estimators that can be practically implemented. If we assume that for all $l$ in $\{1,\cdots,A_{k}\}$ both $\mu^{(l)}$ and the ratio $\mu^{(l)}/\sigma^{(l)}$ do not vary significantly around their respective means, we can think of replacing $\mu^{(l)}$ and $\mu^{(l)}/\sigma^{(l)}$ in (45) by those means. This yields the simplified functional

[TABLE]

to be used in (44). After replacing $(\sigma^{(l)})^{2}$ in (46) by its estimator $\hat{\sigma}_{\hat{K}^{\prime}}^{2}\mathsf{E}_{k,\mathsf{inv}}(\alpha_{k}^{(l)})$, and dropping $\hat{\sigma}_{\hat{K}^{\prime}}^{2}$ because it is independent of $l$, we obtain the proposed (27).

It is interesting to evaluate the loss of performance that results when using (46) instead of (45). We do so by assuming w.l.o.g. that ${\mathcal{E}}_{l}$ holds and estimate the probabilities that a given $n\in\{1,\cdots,A_{k}\}$, $n\neq l$, produces a larger value than for $n=l$ in $\psi^{(n)}$ and $\psi^{\prime(n)}$. Then, we compare the two resulting probabilities in terms of the effective signal-to-noise ratios (SNR). Therefore, in this case, following (40), we have that for $n\neq l$, $\mathsf{\Phi}_{n,\mathsf{inv}}(\alpha^{(n)}_{k})\sim{\mathcal{N}}(0,(\sigma^{(n)})^{2})$, and $\mathsf{\Phi}_{k,\mathsf{inv}}(\alpha^{(l)}_{k})\sim{\mathcal{N}}(\mu^{(l)},(\sigma^{(l)})^{2})$. Thus, when ${\mathcal{E}}_{l}$ holds, $\psi^{(l)}\sim{\mathcal{N}}\left((\mu^{(l)})^{2}/(\sqrt{2}\sigma^{(l)})^{2},(\mu^{(l)})^{2}/(\sigma^{(l)})^{2}\right)$ and $\psi^{(n)}\sim{\mathcal{N}}\left(-(\mu^{(n)})^{2}/(\sqrt{2}\sigma^{(n)})^{2},(\mu^{(n)})^{2}/(\sigma^{(n)})^{2}\right)$, $n\neq l$. Since $\psi^{(l)}$ and $\psi^{(n)}$ are independent, the probability that $\psi^{(n)}\geq\psi^{(l)}$ when ${\mathcal{E}}_{l}$ holds is the probability that the random variable $\psi^{(l)}-\psi^{(n)}$ is less than zero. And since both variables are Gaussian, so is their difference. Therefore, $\psi^{(l)}-\psi^{(n)}\sim{\mathcal{N}}(\omega_{n,l}/2,\omega_{n,l})$, where

[TABLE]

If we define the effective SNR as the ratio between the squared mean and the variance of $\psi^{(l)}-\psi^{(n)}$, then we find that $\text{SNR}_{\psi}=\omega_{n,l}/4$, where the subindex $\psi$ indicates that we are using the estimator in (45).

For the simplified estimator in (46), a similar derivation leads to showing that

[TABLE]

for which the effective SNR, denoted as $\text{SNR}_{\psi^{\prime}}$ is now

[TABLE]

In order to compare the effective SNRs, we compute their ratio:

[TABLE]

To get a cleaner interpretation of this result, we can further assume that both $\mu^{(n)}$ and $(\sigma^{(n)})^{2}$ are proportional to the cardinality of the support set $|{\mathcal{Q}}_{k,\mathsf{inv}}(\alpha_{k}^{(n)})|$. This way, if we let $\beta_{n,l}\doteq|{\mathcal{Q}}_{k,\mathsf{inv}}(\alpha_{k}^{(n)})|/|{\mathcal{Q}}_{k,\mathsf{inv}}(\alpha_{k}^{(l)})|$, we can write that $\mu^{(n)}/\mu^{(l)}=\beta_{n,l}$ and $(\sigma^{(l)})^{2}/(\sigma^{(n)})^{2}=\beta_{n,l}^{-1}$. Then, substituting into (50) we find that

[TABLE]

which is clearly larger than one for all $\beta_{n,l}\geq 0$, $\beta_{n,l}\neq 1$. This confirms that, as expected, for any $\beta_{n,l}\neq 1$ there is a loss of effective SNR with respect to the optimal estimator. However, in practice this loss will be rather small: for instance, suppose that $|{\mathcal{Q}}_{k,\mathsf{inv}}(\alpha_{k}^{(n)})|$ is within the range of 20% larger and 20% smaller than $|{\mathcal{Q}}_{k,\mathsf{inv}}(\alpha_{k}^{(l)})|$, then the effective SNR for the suboptimal detector is at most 0.054 dB smaller than the corresponding to the optimal one.

This supports the plausibility of the proposed simplified detector.