IEEE TRANSACTIONS, TIP-03364-2007-FINAL 1
Practical Poissonian-Gaussian noise modeling and
Þtting for single-image raw-data
Alessandro Foi, Mejdi Trimeche, Vladimir Katkovnik, and Karen Egiazarian, Senior member, IEEE
Abstract— We present a simple and usable noise model for the
raw-data of digital imaging sensors. This signal-dependent noise
model, which gives the pointwise standard-deviation of the noise
as a function of the expectation of the pixel raw-data output, is
composed of a Poissonian part, modeling the photon sensing, and
Gaussian part, for the remaining stationary disturbancies in the
output data. We further explicitly take into account the clipping
of the data (over- and under-exposure), faithfully reproducing
the nonlinear response of the sensor. We propose an algorithm
for the fully automatic estimation of the model parameters given
a single noisy image. Experiments with synthetic images as well
as with real raw-data from various sensors prove the practical
applicability of the method and the accuracy of the proposed
model.
Index Terms— clipping, digital imaging sensors, noise estima-
tion, noise modeling, overexposure, Poisson noise, raw-data.
I. INTRODUCTION
Progress in hardware design and manufacturing has intro-
duced digital imaging sensors having a dramatically increased
resolution. This is mainly achieved by an increase of the
pixel density. Despite the electrical and thermal characteristics
of the sensors have noticeably improved in the last decade
[18], [15], with the size of each pixel becoming smaller and
smaller the sensor output signal’s susceptibility to photon noise
has become greater and greater. As of now, this source of
noise appears as the most signiÞcant contributor of the overall
noise in a digital imaging sensor [1]. This makes the noise
component of the raw-data output of the sensor markedly
signal-dependent, thus far from the conventional additive white
Gaussian noise modeling so widely used in image processing.
Further, with the intention of making full use of the rather
limited dynamic range of digital sensors, pictures are usually
taken with some areas purposely overexposed or clipped,
i.e. accumulating charge beyond the full-well capacity of
the individual pixels. These pixels obviously present highly
nonlinear noise characteristics, which are completely different
than those of normally exposed pixels.
The raw-data which comes from sensor always under-
goes various processing stages (e.g., denoising, demosaicking,
deblurring, compression) before the Þnal “cooked” image
reaches the user. In order to process the data and/or to attenuate
A. Foi, V. Katkovnik, and K. Egiazarian are with the Department of Signal
Processing, Tampere University of Technology, P.O. Box 553, 33101 Tampere,
Finland. E-mail: Þrstname.lastname@tut.Þ
M. Trimeche is with the Multimedia Technologies Laboratory, Nokia
Research Center, Tampere, Finland, E-mail: Þrstname.lastname@nokia.com
This work was supported by the Finnish Funding Agency for Technology
and Innovation (Tekes), AVIPA/AVIPA2 projects, and by the Academy of
Finland, project No. 213462 (Finnish Centre of Excellence program 2006 -
2011).
the noise in the most efÞcient and effective way, it is vital that
a proper modeling of the noise is considered during the various
stages of digital image processing. However, the technical
datasheets of the de vices usually provide vague and inadequate
Þgures for the noise that are of a global nature (i.e., “av erage”
values which are meant to be valid for the whole sensor)
[17]. Consequently, raw-data Þltering algorithms either assume
independent stationary noise models or, if a signal-dependent
model is assumed, the correct parameters for the noise are
often not speciÞed. Such rough noise estimates are inadequate
for the high-quality image processing Þlters which are rapidly
becoming an integral part of the imaging chain.
Two are the contributions in this paper. First, we present
a simple noise model which can accurately be used for the
raw-data. Based on the above considerations, it is a signal-
dependent noise model based on a Poissonian part, modeling
the photon sensing, and Gaussian part, for the remaining sta-
tionary disturbances in the output data. We explicitly take into
account the problem of clipping (over- and under-exposure),
faithfully reproducing the nonlinear response of the sensor.
Only two parameters are sufÞcient to fully describe the model.
These parameters are explained in relation to the sensor’s
hardware characteristics (quantum efÞciency, pedestal, gain).
As a second and most important contribution, we propose
an algorithm for the fully automatic estimation of the model
parameters given a single noisy image.
The paper is organized as follows. In Section II we present
the model in its basic form, which ignores the clipping. The
parameter estimation algorithm is then presented in Section
III. The general model with clipping requires more involved
mathematics, and it is given in Section IV, followed by
the modiÞed estimation algorithm in Sections V and VI.
Throughout these sections, we demonstrate the accurac y of
the algorithm with synthetic test images, for which the exact
noise parameters are known. Experiments with real raw-data
are presented in Section VII; these experiments prove the
practical applicability of the method and conÞrm that the
raw-data noise can indeed be accurately modeled as a clipped
Poissonian-Gaussian process. Further comments and details on
the algorithm and its implementation are giv en in Section VIII.
II. P
OISSONIAN-GAUSSIAN MODELING
Let us consider the generic signal-dependent noise observa-
tion model of the form
z (x)=y ( x)+σ (y (x)) ξ (x) (1)
where x ∈ X is the pixel position in the domain X, z : X →
R is the observed (recorded) signal, y : X → R is the original
2 IEEE TRANSACTIONS, TIP-03364-2007-FINAL
(unknown) signal, ξ : X → R is zero-mean independent
random noise with standard deviation equal to 1,andσ :
R → R
+
is a function of y that gives the standard deviation
of the overall noise component. Throughout the paper, we
denote the expected value (or mathematical expectation) of
a random variable as E {·}, its variance as var {·},andits
standard deviation as std {·} =
p
var {·}; when any of these
operators is applied to a sequence (resp. matrix) of random
variables, its output is deÞned as the sequence (resp. matrix)
of the operator’s outputs for the individual random variables.
The symbol σ is used exclusively to denote this function of
the model (1). From E {ξ (x)} =0follows that E {z (x)} =
y (x), i.e. the original signal can be deÞned as the expected
value of the noisy observations. Consequently, we have that
std {z (x)} = σ (E {z (x)}), i.e. the standard deviation of the
noise is a function, namely σ, of the expectation of the noisy
signal.
In our modeling, we assume that the noise term is com-
posed of two mutually independent parts, a Poissonian signal-
dependent component η
p
and a Gaussian signal-independent
component η
g
:
σ (y (x)) ξ (x)=η
p
(y (x)) + η
g
(x) .(2)
In terms of distributions, these two components are character-
ized as follows,
χ
¡
y (x)+η
p
(y (x))
¢
∼ P (χy (x)) ,η
g
(x) ∼ N (0,b) ,
where χ>0 and b ≥ 0 are real scalar parameters and P and
N denote the Poisson and normal (i.e., Gaussian) distributions.
From the elementary properties of the Poisson distribution, we
obtain the following equation for the mean and variance
E
©
χ
¡
y (x)+η
p
(y (x))
¢ª
=var
©
χ
¡
y (x)+η
p
(y (x))
¢ª
=
= χy (x) .
Since E
©
χ
¡
y (x)+η
p
(y (x))
¢ª
= χy (x)+χE
©
η
p
(y (x))
ª
and χ
2
var
©
η
p
(y (x))
ª
= χy (x),itfollowsthat
E
©
η
p
(y (x))
ª
=0 and var
©
η
p
(y (x))
ª
= y (x) /χ.
Thus, the Poissonian η
p
has varying variance that depends on
the value of y (x), var
©
η
p
(y (x))
ª
= ay (x),wherea = χ
−1
.
The Gaussian component η
g
has instead constant variance
equal to b.
Consequently, the overall variance of z in (1) has the afÞne
form
σ
2
(y (x)) = ay (x)+b,(3)
which giv es the standard deviation σ as the square root
σ (y (x)) =
p
ay (x)+b, (4)
and, in particular, σ (0) =
√
b and σ (1) =
√
a + b.
Some examples of standard-deviation functions σ for dif-
ferent combinations of the constants a and b are sho wn, as an
illustration, in Figure 1 (solid lines).
Figure 2 presents a simple piecewise smooth image which
is degraded by Poissonian and Gaussian noise with parameters
χ = 100 (a =0.01)andb =0.04
2
. As illustrated in Figure
1, these parameters imply that the noise standard-deviation in
the brightest parts of the image is more than twice as large as
in the darker ones.
A. Raw-data modeling
The Poissonian-Gaussian model (1-2) is naturally suited
for the raw-data of digital imaging sensors. The Poissonian
component η
p
models the signal-dependent part of the errors,
which is essentially due to the photon-counting process, while
the Gaussian η
g
accounts for the signal-independent errors
such as electric and thermal noise. We brießymentionhow
the above model parameters relate to elementary aspects of
the digital sensor’s hardware.
1) Quantum efÞciency: The parameter χ of η
p
is related to
the quantum efÞciency of the sensor: the larger the number
of photons necessary to produce a response of the sensor
(generation of an electron), the smaller the χ.
2) Pedestal parameter: In digital imaging sensors, the
collected charge is always added to some base “pedestal” level
p
0
∈ R
+
. This constitutes an offset-from-zero of the output
data and it can be rewritten as a shift in the argument of the
signal-dependent part of the noise:
z (x)=y (x)+σ (y (x) − p
0
) ξ (x)=
= y (x)+η
p
(y (x) − p
0
)+η
g
(x) .
3) Analog gain: We model the analog gain as an ampli-
Þcation of the collected charge. Let us denote the variables
before ampliÞcation by the circle superscript û,
ûz (x)=ûy (x)+ûη
p
(ûy (x) − p
0
)+ûη
g
(x) .
We formalize the ampliÞcation Θ of ûz as the multiplication of
the noise-free signal, of the Poissonian noise, and of a part of
the Gaussian noise, by a scaling constant θ>1,
z(x)=Θ (ûz(x)) = θ
¡
ûy (x)+ûη
p
(ûy(x)−p
0
)+ûη
0
g
(x)
¢
+ûη
00
g
(x).
Here, the Gaussian noise term ûη
g
has been split in two
components ûη
0
g
and ûη
00
g
, ûη
0
g
+ûη
00
g
= η
g
,whereûη
00
g
represents the
portion of the noise that is introduced after the ampliÞcation
and thus not affected by the factor θ. The expectation and
variance for z are
E {z (x)} = y (x)=θûy (x) , var {z (x)} =
= θ
2
χ
−1
(ûy (x) − p
0
)+θ
2
var
©
ûη
0
g
(x)
ª
+var
©
ûη
00
g
(x)
ª
.
Hence, we come again to a model of the form (3)-(4) with
a = χ
−1
θ, b = θ
2
var
©
ûη
0
g
(x)
ª
+var
©
ûη
00
g
(x)
ª
− θ
2
χ
−1
p
0
.
Note that now this b can be negative, provided a large pedestal
p
0
and a small variance of ûη
g
. This does not mean that there is
a “negative” variance. Indeed, because of the pedestal, y ≥ θp
0
and therefore ay + b ≥ 0.
In digital cameras, the analog gain (i.e., θ) is usually
controlled by the choice of the ISO sensitivity setting. This
can be done manually by the user, or automatically by the
camera (“auto mode”). Large ISO numbers (e.g., 800 or 1600)
correspond to large θ, and thus worse signal-to-noise ratio
(SNR). Lower values (e.g., ISO 50) yield a better SNR but at
the same time produce darker images, unless these are taken
with a longer exposure time (which corresponds to having
larger values of ûy before the multiplication by θ).
Figure 3 shows few examples of the standard-deviation
functions σ which can typically be found for the raw data.
Two of these examples have b<0, which corresponds to a
pedestal p
0
> 0.
FOI ET AL., PRACTICAL POISSONIAN-GAUSSIAN NOISE MODELING AND FITTING FOR SINGLE-IMAGE RAW-DATA 3
Fig. 1. Some examples of the standard-deviation functions σ (solid lines) from the model (1) for different combinations of the constants a and b of Equation
(4): (left) a =0.02
2
, 0.06
2
, 0.10
2
, b =0.04
2
and (right) a =0.4
2
, b =0.02
2
, 0.06
2
, 0.10
2
. The dashed lines show the corresponding functions ˜σ of the
clipped observation model (30), as functions of the clipped ˜y (see Section IV). The small black triangles indicate the points (˜y, ˜σ (˜y)) which correspond to
y =0and y =1.
Fig. 2. A piecewise smooth test image of size 512 × 512: original y and
observation z degraded by Poissonian and Gaussian noise with parameters
χ =100(a =0.01)andb =0.04
2
.
B. Heteroskedastic normal approximation
Throughout the following sections, we need to derive a few
results and relations which depend not only on the mean and
variance, but also on the particular distribution of the processed
samples. For the sake of simpliÞcation, we exploit the usual
normal approximation of the Poisson distribution, which gives
P (λ) ≈ N (λ, λ) . (5)
The accuracy of this approximation increases with the para-
meter λ and in practice, for large enough
1
λ, a Poissonian
process can be treated as a special heteroskedastic Gaussian
one. We thus obtain the following normal approximations of
the errors
σ (y (x)) ξ (x)=
p
ay (x)+bξ (x) ' η
h
(y (x)) , (6)
where η
h
(x) ∼ N (0,ay(x)+b).
III. T
HE ALGORITHM
Our goal is to estimate the function σ : R → R
+
of the
observation model (1) from a noisy image z . The proposed
algorithm is divided in two main stages: local estimation
of multiple expectation/standard-deviation pairs and global
parametric model Þtting to these local estimates. An initial
1
How large λ is enough really depends on the considered application
and desired accuracy. The fact that the Poisson distribution is discrete is
a secondary aspect, because quantization of the digital data makes anyway
discrete ev en errors due to continuous distributions.
For the considered standard-deviation estimation problem, we found exper-
imentally that already with λ =10(corresponding to χ =20for the middle
intensity y =0.5) there is virtually no difference between the estimation
accuracy of a truly Poissonian variable and that of its Gaussian approximation.
Fig. 3. Some examples of the standard-deviation functions σ (solid lines)
which are often found for the raw data. In these three examples the parameters
(a, b) from Equation (4) are
!
1.5 · 10
−3
, 10
−4
"
,
!
6 · 10
−4
, −5 · 10
−6
"
,and
!
10
−4
, −8 · 10
−6
"
. The dashed lines show the corresponding functions ˜σ of
the clipped observation model (30), as functions of the clipped ˜y (see Section
IV). The small black triangles indicate the points (˜y, ˜σ (˜y)) which correspond
to y =0and y =1.
preprocessing stage, in which the data is transformed to the
wavelet domain and then segmented into non-overlapping level
sets where the data is smooth, precedes the estimation.
A. Wavelet domain analysis
Similar to [4], we facilitate the noise analysis by considering
wavelet detail coefÞcients z
wdet
deÞned as the downsampled
convolution
z
wdet
=↓
2
(z ~ ψ) ,
where ψ is a 2-D wavelet function with zero mean and unity
5
2
-norm,
P
ψ =0, kψk
2
=1,and↓
2
denotes the decimation
operator that discards ev e ry second row and every second
column. Analogously, we deÞne the normalized approximation
coefÞcients as
z
wapp
=↓
2
(z ~ ϕ) ,
where ϕ is the corresponding 2-D wavelet scaling function,
which we specially normalize so that
P
ϕ =1.
For noisy images, the detail coefÞcients z
wdet
contain mostly
noise and, due to the normalizations of the convolution kernels,
we have
std
©
z
wdet
ª
= ↓
2
(std {z ~ ψ})=↓
2
µ
q
var {z} ~ ψ
2
¶
' (7)
' ↓
2
(std {z}kψk
2
)=↓
2
(std {z})=
= ↓
2
(σ (y)) = σ (↓
2
y)=σ (↓
2
(y
P
ϕ)) '
' σ (↓
2
(y ~ ϕ)) = σ (E {z
wapp
}) ,(8)
4 IEEE TRANSACTIONS, TIP-03364-2007-FINAL
Fig. 4. From left to right: wa velet approximation and detail coefÞcients z
wapp
and z
wdet
, restricted on the set of smoothness X
smo
,andtwolevel-setsS
i
(13) computed for ∆
i
= ∆ =1/300. The scale of this Þgure is half that of Fig. 2.
with the approximate equalities ' becoming accurate at points
in regions where y (and hence std {z}) is uniform, as we can
assume that the distribution of z does not change over the
small support of the wavelets. Thus, in particular, at a point
x in such uniform regions, we can assume that
z
wdet
(x) ∼ N (0,σ(E {z
wapp
(x)})) ,(9)
and, because of decimation and orthogonality properties of
wavelet functions, that the noise degrading z
wdet
,aswellas
the noise degrading z
wapp
, are independent ones.
Note that, always, kϕk
2
6=1. Therefore, when considering
std {z
wapp
}, the above equations can be repeated, replacing ψ
with ϕ, only provided that the factor kϕk
2
is kept. Thus, we
come to
std {z
wapp
}'kϕk
2
σ (z
wapp
) .
In our implementation, we use separable kernels ψ = ψ
1
~
ψ
T
1
and ϕ = ϕ
1
~ ϕ
T
1
where ψ
1
and ϕ
1
are 1-D Daubechies
wavelet and scaling functions
ψ
1
=[0.035 0.085 − 0.135 − 0.460 0.807 − 0.333] , (10)
ϕ
1
=[0.025 − 0.060 − 0.095 0.325 0.571 0.235] .
B. Segmentation
Like in our previous work [6], we segment the data into
level sets, in each of which the image can be reasonably
assumedtobeuniformlyclosetoacertainvalue.Having
nothing but a noisy image a t our disposal, we shall employ
spatial smoothing (as opposed to temporal smoothing, used in
[6]) in order to attenuate the noise and an edge-detector in
order to stay clear from edges when analyzing the data, thus
enabling the conditions (7)-(8).
There e xist a myriad of different methods which can be
used for smoothing or for edge detection. However, for our
purposes, the follo wing simple and non-adaptive methods
proved adequate for all considered experimental cases.
1) Smoothed approximation: From z
wapp
, we compute a
smoothed (low-pass) image z
smo
,
z
smo
= z
wapp
~ 7, (11)
where 7 is positive smoothing kernel, 7 ≥ 0 and k7k
1
=
1. The smoothing action of the k ernel should be especially
strong, so to effectively suppress most of the noise. In our
implementation, we use a uniform 7 × 7 kernel for 7.
In the corresponding regions where y itself is smooth, z
smo
is approximately equal to E {z
wapp
}, and thus to ↓
2
y.Thisis
a reasonable assumption provided that the support of 7 does
not intersect edges during the calculation of the convolution
(11).
2) Edges and set of smoothness: To detect edges, we use
the conventional approach where some smoothed derivatives
of the image are thresholded against an estimate of the local
standard deviation. Exploiting the fact that the mean of the
absolute deviations of N (0, 1) is equal to
p
2/π [7], we can
deÞne a rough estimate of the local standard-de viations of z
wdet
as the map
s =
r
π
2
¯
¯
z
wdet
¯
¯
~ 7.
We deÞne the set of smoothness X
smo
as
X
smo
= {x ∈↓
2
X :
: |∇(Λ (z
wapp
)) (x)| + |Λ (z
wapp
)(x)| <τ· s (x)}, (12)
Λ (z
wapp
)=∇
2
medÞlt (z
wapp
) ,
where ∇ and ∇
2
are, respectively, gradient and Laplacian
operators, medÞlt denotes a 3 × 3 median Þlter, ↓
2
X is the
decimated domain of the wavelet coefÞcients z
wapp
,andτ>0
is positive threshold constant. We realize both the Laplacian
operator ∇
2
and the gradient operator ∇ as convolutions
against 9 × 9 kernels. Thresholding the sum of the moduli of
the Laplacian and of its gradient is a heuristic way to obtain
“thickened” edges.
In Figure 4, we show the wavelet approximation and detail
coefÞcients z
wapp
and z
wdet
, restricted on the set of smoothness
X
smo
(whose complement thus appears as white in the Þgure),
calculated for the test image z of Figure 2. Note that some of
the weakest edges have not been detected as such.
3) Level sets (segments): In the set of smoothness X
smo
,
we can assume that edges of the image did not interfere with
the smoothing (11), hence, that the conditions (7)-(8) hold and
that, for x ∈ X
smo
,
z
smo
(x)=E {z
wapp
(x)} = E {(↓
2
z)(x)} =(↓
2
y)(x) ,
std
©
z
wdet
(x)
ª
=std{(↓
2
z)(x)} =(↓
2
(σ (y)))(x).
We identify in the smoothness set X
smo
acollectionof
N non-overlapping level sets (segments) S
i
⊂ X
smo
, i =
1,...,N of the smoothed image z
smo
. Each level set, char-
acterized by its centre v alue u
i
and allowed deviation ∆
i
> 0,
is deÞned as
S
i
= {x ∈ X
smo
: z
smo
(x) ∈ [u
i
− ∆
i
/2,u
i
+ ∆
i
/2)} . (13)
FOI ET AL., PRACTICAL POISSONIAN-GAUSSIAN NOISE MODELING AND FITTING FOR SINGLE-IMAGE RAW-DATA 5
By non-overlapping we mean that S
i
∩ S
j
= ∅ if
i 6= j. In practice, assuming a signal in the range [0, 1],
one can take Þxed ∆
i
≡ ∆ and equispaced u
i
∈
©
∆j, j =1,...,
¯
N =
§
∆
−1
¨ª
,wheretheb·c brackets indi-
cate the rounding to the nearest larger or equal integer. Further,
we require that the le vel sets are non-trivial, in the sense that
each set S
i
mustcontainatleasttwosamples
2
; thus, N ≤
¯
N
and
S
N
i=1
S
i
⊆ X
smo
⊆↓
2
X. Figure 4 shows two of the level
sets computed for the example in Figure 2 for ∆ =1/300.
Observe that these sets are meager and quite fragmented.
C. Local estimation of expectation/standard-deviation pairs
For each level set S
i
,wedeÞne the (unknown) variable
y
i
=
1
n
i
n
i
X
j=1
E {z
wapp
(x
j
)} , {x
j
}
n
i
j=1
= S
i
. (14)
Note that y
i
and u
i
might not coincide. The level set S
i
is used
as a one domain for the computation of a pair of estimates
(ˆy
i
, ˆσ
i
),whereˆy
i
is an estimate of y
i
and ˆσ
i
is an estimate
of σ (y
i
). In what follows, although we shall refer explicitly
to y
i
, this variable is always used implicitly and, in the Þnal
estimation of the function y 7→ σ (y),themanyy
i
, i,...,N,
remain “hidden” variables which are modeled as unknown.
Similarly, the smoothed data z
smo
and the values u
i
and ∆
i
used for the construction of S
i
do not appear in the following
estimation, where only z
wdet
, z
wapp
and S
i
areusedinorder
to compute the estimates ˆy
i
and ˆσ
i
.
1) Estimation of y
i
: We estimate y
i
as the sample mean of
the approximation coefÞcients z
wapp
on S
i
ˆy
i
=
1
n
i
n
i
X
j=1
z
wapp
(x
j
) , {x
j
}
n
i
j=1
= S
i
. (15)
2) Estimation of σ (y
i
): The estimate ˆσ
i
is calculated as the
unbiased sample standard-deviation of the detail coefÞcients
z
wdet
on S
i
ˆσ
i
=
1
κ
n
i
s
P
n
i
j=1
¡
z
wdet
(x
j
) − ¯z
wdet
i
¢
2
n
i
− 1
, (16)
where ¯z
wdet
i
=
1
n
i
P
n
i
j=1
z
wdet
(x
j
) and the factor κ
−1
n
i
is
deÞned [7]
κ
n
=
r
2
n − 1
Γ
¡
n
2
¢
Γ
¡
n−1
2
¢
=1−
1
4n
−
7
32n
2
+ O
µ
1
n
3
¶
. (17)
This factor, which comes from the mean of the chi-distribution
with n − 1 degrees of freedom, makes the estimate unbiased
for normally and identically independently distributed (i.i.d.)
z
wdet
(x
j
).
3) Unbiasedness: Clearly from the deÞnition (14), ˆy
i
is an
unbiased estimator of y
i
.
The unbiasedness of ˆσ
i
as an estimator of σ (y
i
) is a
more complex issue. As observed above, ˆσ
i
is an unbiased
estimator of σ (y
i
) provided that z
wdet
is normally i.i.d. on
the level set S
i
. However, we cannot claim, in general,
that z
wdet
is identically distributed on S
i
.Weremarkthat
2
The smoothness threshold τ (12) can be automatically increased in the
rare event of N<2, i.e. when there are not enough non-trivial level sets for
the estimation. Note that X
smo
is monotonically enlarging to ↓
2
X with τ,
X
smo
%
τ →∞
↓
2
X.
the assumed validity of (7)-(8) concerns individual points.
It does not mean that std
©
z
wdet
ª
is constant over S
i
.As
a matter of fact, especially for large ∆
i
, E {z
wapp
(x)} is
not constant for x ∈ S
i
, which implies that the standard
deviations of the wavelet detail coefÞcients (8) are not constant
over S
i
. Lacking any particular hypothesis on the image y,
it is nevertheless reasonable to assume that {E {z
wapp
(x)},
x ∈ S
i
} has a symmetric (discrete) distribution centred at y
i
(with diameter bounded by ∆
i
). Because of (3) and (9), we
have that {var
©
z
wdet
(x)
ª
, x ∈ S
i
} has also a symmetrical
distribution, which is centred at var {y
i
}. This makes κ
2
n
i
ˆσ
2
i
an unbiased
3
estimator of var {y
i
} and, since κ
n
→
n→∞
1,
ˆσ
i
is an asymptotically unbiased estimator of σ (y
i
).This
asymptotic unbiasedness is relevant in the practice, since a
large ∆
i
corresponds to large n
i
.
We further note that, despite the segmentation and removal
of edges, the presence of sharp image features, singularities, or
e ven texture in the segment S
i
is not completely ruled out. This
can be effectively compensated by means of non-linear robust
estimators of the standard de viation, such as the well-known
median of absolute deviations (MAD) [12]. F or the sake of
expository simplicity, in the current and in the next section
we restrict ourself to the basic estimator (16) and postpone
considerations on robust estimation of the standard-deviation
to Section VI.
4) Variance of the estimates: The variance of the estimates
ˆy
i
and ˆσ
i
depends directly on the variances of the samples
used for the estimation, which are degraded by independent
noise. With arguments similar to Section III-C.3, the variances
of the estimates can be expressed as
var {ˆy
i
} = σ
2
(y
i
) c
i
, var {ˆσ
i
} = σ
2
(y
i
) d
i
, (18)
c
i
=
kϕk
2
2
n
i
,d
i
=
1 − κ
2
n
i
κ
2
n
i
=
1
2n
i
+
5
8n
2
i
+ O
µ
1
n
3
i
¶
, (19)
where these expressions coincide with those for the perfect
case when var
©
z
wdet
ª
and var {z
wapp
} are constant on S
i
[7].
5) Distribution of the estimates: The estimates ˆy
i
and ˆσ
i
are distributed, respectively, following a normal distribution
and a scaled non-central chi-distribution, which can also
be approximated, very accurately for large n
i
,asanormal
distribution [7]. Thus, in what follows, we treat both ˆy
i
and
ˆσ
i
as normally distributed random variables and, in particular,
as
ˆy
i
∼ N
¡
y
i
,σ
2
(y
i
) c
i
¢
, ˆσ
i
∼ N
¡
σ (y
i
) ,σ
2
(y
i
) d
i
¢
, (20)
where c
i
and d
i
are deÞned as in (19).
D. Maximum-likelihood Þtting of a global parametric model
The maximum-likelihood (ML) approach is used to Þta
global parametric model of the function σ on the estimates
{ˆy
i
, ˆσ
i
}
N
i=1
. Depending on the parameters a and b,wehave
σ
2
(y )=ay + b. For reasons of numerical consistency (note
that formally this σ
2
(y) may be zero or negative), for the
Þtting we deÞne a simple regularized variance-function σ
2
reg
as
σ
2
reg
(y)=max
¡
ε
2
reg
,σ
2
(y)
¢
(21)
3
This can be proved easily since, for x ∈ X
smo
,wecantreatE
#
z
wdet
(x)
$
as zero.
