This paper was submitted to IEEE Transactions on Information Forensics &
Security the July 28, 2022.

# An Adaptive Method for Camera Attribution under Complex Radial Distortion
Corrections

Andrea Montibeller, and Fernando Pérez-González A. Montibeller is with the
Department of Information Engineering and Computer Science, University of
Trento, Trento 38123, Italy (e-mail: andrea.montibeller@unitn.it).F. Pérez-
González is with the atlanTTic Research Center, Department of Signal Theory
and Communications, University of Vigo, Vigo 36310, Spain (e-mail:
fperez@gts.uvigo.es).Manuscript received xxx xx, xxxx; revised xxx xx, xxxx.

###### 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.

###### Index Terms:

Image forensics, source attribution, PRNU, photo response non-uniformity,
radial correction, distortion correction, PCE, adaptive processing.

##  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

$\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$}}||},$ (1)

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]:

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

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

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

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]:

${\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}.$ (4)

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

$\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$}))},$ (5)

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

$\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}},$ (6)

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

$\displaystyle T_{\alpha}:\mathbb{R}^{2}$ $\displaystyle\rightarrow$
$\displaystyle\mathbb{R}^{2}$ $\displaystyle(x,y)$ $\displaystyle\mapsto$
$\displaystyle(x^{\prime},y^{\prime})$ (7)

where

$x^{\prime}=x_{p}+(x-x_{p})(1+\alpha r^{2});\\\ $ (8)
$y^{\prime}=y_{p}+(y-y_{p})(1+\alpha r^{2}),\\\ $ (9)

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

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

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

$r^{\prime}=T_{\bm{\alpha}}(r)=r\left(1+\sum_{i=1}^{n}\alpha_{i}r^{2i}\right),\\\
$ (11)

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

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

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

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

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

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

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

$\text{PCE}_{\mathsf{max}}^{*}\doteq\max_{\bm{\alpha}\in{\mathcal{A}}}\text{PCE}_{\mathsf{max}}(\bm{\alpha}).$
(15)

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

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

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.

Figure 1: Values of $\alpha$ maximizing $\text{PCE}_{\mathsf{inv}}(\alpha)$ vs
inner radius of the annulus. Values are linearly interpolated. Canon 1200D
camera with EF-S 10-18mm lens, corrected with Adobe Lightroom. Focal length:
10mm. Shutter speed: 1/100 sec. Aperture: f7.1. ISO 800. The PRNU is estimated
with 20 natural images all taken with those settings.

##  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.

Figure 2: $\text{PCE}_{\mathsf{max}}$ as a function of $\alpha$ for a
Panasonic DMC-ZS7 camera.

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^{*}$.

Figure 3: $\text{PCE}_{\mathsf{inv}}(\alpha)$ for: $\alpha=-0.01$ and
$\alpha=0.05$.

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:

${\mathcal{R}}_{k}\doteq\\{(u,v)\in{\mathbb{R}}^{2}:r_{k}^{2}\leq
u^{2}+v^{2}<(r_{k}+\Delta_{k})^{2}\\}.$ (17)

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$.

Figure 4: Annular partition used in the proposed method.

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

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

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

$\displaystyle T_{\alpha_{k},c}(r)$ $\displaystyle\doteq$ $\displaystyle
r(1+\alpha_{k}\cdot r^{2});\ \ $ $\displaystyle T_{\alpha_{k},l}(r)$
$\displaystyle\doteq$ $\displaystyle r(1+\alpha_{k}),\ \
r\in{\mathcal{R}}_{k},$ (19)

whereas the corresponding inverse mappings are

$\displaystyle T_{\alpha_{k},c}^{-1}(r^{\prime})$ $\displaystyle\approx$
$\displaystyle r^{\prime}(1-\alpha_{k}\cdot r^{\prime
2}+3\alpha_{k}^{2}r^{\prime 4}),\ \ r^{\prime}\in
T_{\alpha_{k},c}({\mathcal{R}}_{k});\ \ $ $\displaystyle
T_{\alpha_{k},l}^{-1}(r^{\prime})$ $\displaystyle=$
$\displaystyle\frac{r^{\prime}}{1+\alpha_{k}},\ \ r^{\prime}\in
T_{\alpha_{k},l}({\mathcal{R}}_{k}),$ (20)

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.,

${\mathcal{P}}_{k}\doteq\left(D_{2}\cdot{\mathcal{R}}_{k}\right)\cap\mathcal{I},\
\ k=1,\cdots,L,$ (21)

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).

Figure 5: Illustration of the application of transforms $T^{-1}_{\alpha_{k}}$
and $T_{\alpha_{k}}$, and related domains.

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

${\mathcal{Q}}_{k,\mathsf{inv}}(\alpha_{k})=D_{2}\cdot
T_{\alpha_{k}}\left(\left[\left(D_{2}\cdot
T^{-1}_{\alpha_{k}}({\mathcal{P}_{k}}/D_{2})\right)\cap{\mathcal{B}}\right]/D_{2}\right).$
(22)

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

${\mathcal{Q}}_{k,\mathsf{dir}}(\alpha_{k})=D_{2}\cdot
T^{-1}_{\alpha_{k}}\left(\left[\left(D_{2}\cdot
T_{\alpha_{k}}({\mathcal{P}_{k}}/D_{2})\right)\cap{\mathcal{B}}\right]/D_{2}\right).$
(23)

###  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,

$\displaystyle\mathsf{\Phi}_{k,\mathsf{inv}}(\alpha_{k})$
$\displaystyle\doteq$
$\displaystyle\sum_{(i,j)\in{\mathcal{Q}}_{k,\mathsf{inv}}(\alpha_{k})}\hat{K}^{\prime}_{i,j}\cdot\left[T^{-1}_{\alpha_{k}}(\tilde{\mbox{$\bf
W$}})\right]_{i,j},$ (24)
$\displaystyle\mathsf{E}_{k,\mathsf{inv}}(\alpha_{k})$ $\displaystyle\doteq$
$\displaystyle\sum_{(i,j)\in{\mathcal{Q}}_{k,\mathsf{inv}}(\alpha_{k})}\left[T^{-1}_{\alpha_{k}}(\tilde{\mbox{$\bf
W$}})\right]_{i,j}^{2},$ (25)

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

$\alpha^{*}_{k}=\arg\max_{\alpha_{k}\in{\mathcal{A}}_{k}}\varphi_{k,\mathsf{inv}}(\alpha_{k}),$
(26)

where

$\varphi_{k,\mathsf{inv}}(\alpha_{k})\doteq\frac{\mathsf{\Phi}_{k,\mathsf{inv}}(\alpha_{k})}{\mathsf{E}_{k,\mathsf{inv}}(\alpha_{k})}.$
(27)

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

$k_{0}=\arg\max_{k=1,\cdots,L}\varphi_{k,\mathsf{inv}}(0).$ (28)

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

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

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

$e_{k}\doteq\alpha_{k}^{*}-\hat{\alpha}_{k},$ (30)

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,

$\mbox{$\bf u$}_{k}=\mbox{$\bf u$}_{k-1}+\mu e_{k}\bm{\beta}_{k-1},\ \
k=k_{0}+1,\cdots,L,$ (31)

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:

$\mbox{$\bf u$}_{k}=\mbox{$\bf u$}_{k+1}+\mu e_{k}\bm{\beta}_{k+1},\ \
k=k_{0}-1,\cdots,1,$ (32)

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:

${\mathcal{A}}_{k}=\\{\hat{\alpha}_{k}+\lambda_{k}\cdot
n:n\in{\mathbb{Z}}\cap[-(A_{k}-1)/2,(A_{k}-1)/2]\\}.$ (33)

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

$\text{ssq}(\langle\hat{\mbox{$\bf
K$}}^{\prime},T^{-1}_{\bm{\alpha}^{*}}(\tilde{\mbox{$\bf
W$}})\rangle)=\sum_{k=1}^{L}\text{ssq}(\mathsf{\Phi}_{k,\mathsf{inv}}(\alpha^{*}_{k})),$
(34)

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

$\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}\approx\kappa\cdot\hat{\sigma}_{\hat{K}^{\prime}}^{2}\cdot\sum_{k=1}^{L}\mathsf{E}_{k,\mathsf{inv}}(\alpha^{*}_{k}),$
(35)

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.

Figure 6: Sample variance of the transformed residual for different annuli and
different values of $\alpha$. Camera and parameters are the same as in Fig. 1.

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,

$\text{CPCE}_{n,\mathsf{inv}}(\hat{\mbox{$\bf K$}}^{\prime},\mbox{$\bf
W$})\doteq\frac{\sum_{k=\xi(1)}^{\xi(n)}\text{ssq}\left(\mathsf{\Phi}_{k,\mathsf{inv}}(\alpha^{*}_{k})\right)}{\hat{\sigma}_{\hat{K}^{\prime}}^{2}\left(\sum\limits_{k=\xi(1)}^{\xi(n)}\mathsf{E}_{k,\mathsf{inv}}(\alpha^{*}_{k})+\sum\limits_{k=\xi(n+1)}^{\xi(L)}\mathsf{E}_{k,\mathsf{inv}}(\alpha^{*}_{\mathsf{f}})\right)}.$
(36)

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,

$\lambda_{k}=\begin{cases}0.1&\text{if }|\alpha_{k}-\alpha_{k-1}|>0.1,\\\
0.01&\text{if }0.01<|\alpha_{k}-\alpha_{k-1}|\leq 0.1,\\\ 0.001&\text{if
}|\alpha_{k}-\alpha_{k-1}|\leq 0.001.\end{cases}$ (37)

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).

| $\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 I: Thresholds required to achieve FPR=0.05 and FPR=0.01 and
corresponding TPRs for different methods and variants.

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.

|  $\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

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. Figure 7: ROCs obtained with the variants of our method, [8]
(with and without DS) and [9] on the test dataset. For a better visualization,
the zoomed-in box corresponding to low FPRs uses a log scale on the x-axis.

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:

$[T^{-1}_{\alpha_{k}^{(n)}}(\tilde{\mbox{$\bf
W$}})]_{i,j}=\gamma_{i,j}^{(n)}[T^{-1}_{\alpha^{(n)}_{k}}(T_{\alpha^{\dagger}_{k}}(\hat{\mbox{$\bf
K$}}^{\prime}))]_{i,j}+N_{i,j}^{(n)}$ (38)

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

$\sum_{(i,j)\in{\mathcal{Q}}_{k,\mathsf{inv}}(\alpha_{k}^{(n)})}\left(\gamma_{i,j}^{(n)}\right)^{2}[T^{-1}_{\alpha^{(n)}_{k}}(T_{\alpha^{\dagger}_{k}}(\hat{\mbox{$\bf
K$}}^{\prime}))]_{i,j}^{2}\ll\mkern-36.0mu\sum_{(i,j)\in{\mathcal{Q}}_{k,\mathsf{inv}}(\alpha_{k}^{(n)})}\mkern-36.0muN_{i,j}^{2}$
(39)

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

$\displaystyle\mathsf{\Phi}_{k,\mathsf{inv}}(\alpha^{(n)}_{k})$
$\displaystyle\sim$ $\displaystyle{\mathcal{N}}(0,(\sigma^{(n)})^{2}),\
n=1,\cdots,A_{k},n\neq l$ (40)
$\displaystyle\mathsf{\Phi}_{k,\mathsf{inv}}(\alpha^{(l)}_{k})$
$\displaystyle\sim$ $\displaystyle{\mathcal{N}}(\mu^{(l)},(\sigma^{(l)})^{2})$
(41)

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

$\displaystyle
f(\mathsf{\Phi}_{k,\mathsf{inv}}(\alpha^{(1)}_{k}),\cdots,\mathsf{\Phi}_{k,\mathsf{inv}}(\alpha^{(A_{k})}_{k})|{\mathcal{E}}_{l})=$
$\displaystyle
f_{\mathcal{N}}(\mathsf{\Phi}_{k,\mathsf{inv}}(\alpha^{(l)}_{k});\mu^{(l)},\sigma^{(l)})\cdot\prod_{\stackrel{{\scriptstyle
n=1}}{{n\neq
l}}}^{A_{k}}f_{\mathcal{N}}(\mathsf{\Phi}_{k,\mathsf{inv}}(\alpha^{(l)}_{k});0,\sigma^{(n)})$
(42)

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

$\displaystyle
L(\mathsf{\Phi}_{k,\mathsf{inv}}(\alpha^{(1)}_{k}),\cdots,\mathsf{\Phi}_{k,\mathsf{inv}}(\alpha^{(A_{k})}_{k})|{\mathcal{E}}_{l})$
$\displaystyle=\frac{f_{\mathcal{N}}(\mathsf{\Phi}_{k,\mathsf{inv}}(\alpha^{(l)}_{k});\mu^{(l)},\sigma^{(l)})}{f_{\mathcal{N}}(\mathsf{\Phi}_{k,\mathsf{inv}}(\alpha^{(l)}_{k});0,\sigma^{(l)})}$
(43)

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

$l^{*}=\arg\max_{l=1,\cdots,A_{k}}\mathsf{\psi}^{(l)}$ (44)

where

$\mathsf{\psi}^{(l)}\doteq\frac{\mu^{(l)}\cdot\mathsf{\Phi}_{k,\mathsf{inv}}(\alpha^{(l)}_{k})}{(\sigma^{(l)})^{2}}-\frac{1}{2}\frac{(\mu^{(l)})^{2}}{(\sigma^{(l)})^{2}}$
(45)

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

$\mathsf{\psi}^{\prime(l)}\doteq{\mathsf{\Phi}_{k,\mathsf{inv}}(\alpha^{(l)}_{k})}/{(\sigma^{(l)})^{2}}$
(46)

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

$\omega_{n,l}\doteq\frac{(\mu^{(l)})^{2}}{(\sigma^{(l)})^{2}}+\frac{(\mu^{(n)})^{2}}{(\sigma^{(n)})^{2}}$
(47)

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

$\psi^{\prime(l)}-\psi^{\prime(n)}\sim{\mathcal{N}}\left(\frac{\mu^{(l)}}{(\sigma^{(l)})^{2}},\left[\frac{1}{(\sigma^{(l)})^{2}}+\frac{1}{(\sigma^{(n)})^{2}}\right]\right)$
(48)

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

$\text{SNR}_{\psi^{\prime}}=\frac{(\mu^{(l)})^{2}/(\sigma^{(l)})^{2}}{\frac{(\sigma^{(l)})^{2}}{(\sigma^{(n)})^{2}}+1}$
(49)

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

$\frac{\text{SNR}_{\psi}}{\text{SNR}_{\psi^{\prime}}}=\frac{1+\left(\frac{\mu^{(n)}}{\mu^{(l)}}\right)^{2}\cdot\left(\frac{\sigma^{(l)}}{\sigma^{(n)}}\right)^{2}}{4}\cdot\left(\left(\frac{\sigma^{(l)}}{\sigma^{(n)}}\right)^{2}+1\right)$
(50)

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

$\frac{\text{SNR}_{\psi}}{\text{SNR}_{\psi^{\prime}}}=\frac{(1+\beta_{n,l})^{2}}{4\beta_{n,l}}=1+\frac{(1-\beta_{n,l})^{2}}{4\beta_{n,l}}$
(51)

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.

## References

  * [1] J. Lukas, J. Fridrich, and M. Goljan, “Digital camera identification from sensor pattern noise,” IEEE Transactions on Information Forensics and Security, vol. 1, no. 2, pp. 205–214, 2006.
  * [2] K. Rosenfeld and H. T. Sencar, “A study of the robustness of prnu-based camera identification,” in Media Forensics and Security, vol. 7254, p. 72540M, International Society for Optics and Photonics, 2009.
  * [3] S. Taspinar, M. Mohanty, and N. Memon, “Camera fingerprint extraction via spatial domain averaged frames,” IEEE Transactions on Information Forensics and Security, vol. 15, pp. 3270–3282, 2020.
  * [4] P. Korus and J. Huang, “”Multi-Scale Analysis Strategies in PRNU-Based Tampering Localization”,” IEEE Transactions on Information Forensics and Security, vol. 12, no. 4, pp. 809–824, 2017.
  * [5] M. Goljan, “Digital camera identification from images–estimating false acceptance probability,” in International workshop on digital watermarking, pp. 454–468, Springer, 2008.
  * [6] S. Taspinar, M. Mohanty, and N. Memon, “Source camera attribution using stabilized video,” in 2016 IEEE International Workshop on Information Forensics and Security (WIFS), pp. 1–6, IEEE, 2016.
  * [7] D. Morshedi, M. Hosseini, and M. Goljan, “Camera identification from hdr images,” in Proceedings of the ACM Workshop on Information Hiding and Multimedia Security, pp. 69–76, 2019.
  * [8] M. Goljan and J. Fridrich, “Sensor-fingerprint based identification of images corrected for lens distortion,” in Media Watermarking, Security, and Forensics 2012, vol. 8303, p. 83030H, International Society for Optics and Photonics, 2012.
  * [9] M. Goljan and J. Fridrich, “Estimation of lens distortion correction from single images,” in Media Watermarking, Security, and Forensics 2014, vol. 9028, p. 90280N, International Society for Optics and Photonics, 2014.
  * [10] M. Chen, J. Fridrich, M. Goljan, and J. Lukás, “Determining image origin and integrity using sensor noise,” IEEE Transactions on information forensics and security, vol. 3, no. 1, pp. 74–90, 2008.
  * [11] M. K. Mihcak, I. Kozintsev, K. Ramchandran, and P. Moulin, “Low-complexity image denoising based on statistical modeling of wavelet coefficients,” IEEE Signal Processing Letters (SPL), vol. 6, pp. 300–303, 1999.
  * [12] M. Goljan, J. Fridrich, and T. Filler, “Large scale test of sensor fingerprint camera identification,” Proceedings of SPIE - The International Society for Optical Engineering, February 2009.
  * [13] X. Kang, Y. Li, Z. Qu, and J. Huang, “Enhancing source camera identification performance with a camera reference phase sensor pattern noise,” IEEE Transactions on Information Forensics and Security, vol. 7, no. 2, pp. 393–402, 2012.
  * [14] W. Hugemann, “Correcting lens distortions in digital photographs,” Ingenieurbüro Morawski+ Hugemann: Leverkusen, Germany, vol. 20, 2010.
  * [15] H. Li and R. Hartley, “A non-iterative method for correcting lens distortion from nine point correspondences,” OMNIVIS 2005, vol. 2, p. 7, 2005.
  * [16] J. Perš and S. Kovacic, “Nonparametric, model-based radial lens distortion correction using tilted camera assumption,” in Proceedings of the Computer Vision Winter Workshop, vol. 1, pp. pp–286, 2002.
  * [17] P. Drap and J. Lefèvre, “An exact formula for calculating inverse radial lens distortions,” Sensors, vol. 16, no. 6, p. 807, 2016.
  * [18] M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables.  New York: Dover, ninth dover printing, tenth gpo printing ed., 1964.
  * [19] A. Montibeller and F. Pérez-González, “Technical report additional material for “an adaptive method for camera attribution under complex radial distortion corrections”,” 2022.
http://dx.doi.org/10.13140/RG.2.2.20038.96323.

  * [20] B. Widrow and S. Stearns, Adaptive Signal Processing.  Englewood Cliffs, NJ: Prentice-Hall, 1985.
  * [21] F. Pérez-González, M. Masciopinto, I. González-Iglesias, and P. Comesaña, “Fast sequential forensic detection of camera fingerprint,” in 2016 IEEE International Conference on Image Processing (ICIP), pp. 3902–3906, IEEE, 2016.
  * [22] M. Goljan, “Blind detection of image rotation and angle estimation,” Electronic Imaging, vol. 2018, no. 7, pp. 158–1, 2018.
  * [23] S. Mandelli, P. Bestagini, L. Verdoliva, and S. Tubaro, “Facing device attribution problem for stabilized video sequences,” IEEE Transactions on Information Forensics and Security, vol. 15, pp. 14–27, 2019.
  * [24] M. Iuliani, M. Fontani, and A. Piva, “A leak in prnu based source identification—questioning fingerprint uniqueness,” IEEE Access, vol. 9, pp. 52455–52463, 2021.