A robust algorithm for sky background computation in CCD images

12 Pages·0.48 MB·English

by  F. Patat

Astronomy & Astrophysics manuscript no. preprint February 2, 2008
(DOI: will be inserted by hand later)
A robust algorithm for sky background computation in CCD
⋆
images
F. Patat
European Southern Observatory,K.Schwarzschild Str.2, 85748 - Garching - Germany
3
0
0 Received October 24, 2002; accepted Januar 24, 2003
2
n Abstract. In this paper we present a non-interactive algorithm to estimate a representative value for the sky
a backgroundonCCDimages.Themethodwehavedevisedusesthemodeasarobustestimatorofthebackground
J brightness in sub-windows distributed across theinput frame. The presenceof contaminating objects is detected

7 throughthestudyofthelocalintensitydistributionfunctionandtheperturbedareasarerejectedusingastatistical
2 criterion which was derived from numerical simulations. The technique has been extensively tested on a large
amount of images and it is suitable for fully automatic processing of large data volumes. The implementation

1 we discuss here has been optimized for the ESO-FORS1instrument, but it can be easily generalized to all CCD
v imagers with a sufficiently large field of view. The algorithm has been successfully used for the UBVRI ESO-
4
Paranal night sky brightness survey (Patat 2003).
3
5
Key words.techniques:photometric
1
0
3
0 1. Introduction changeablecollimators,whichgiveanimagingfieldofview
h/ of 6.′8 6.′8 (standard resolution, SR) and 3.′4 3.′4 (high
The study of the night sky brightness is fundamental to × ×
p resolution, HR) respectively.
- monitor the quality of very dark astronomical sites and
o sets the stage for the study of a whole series of effects FORS1 is offered during dark time both in Visitor
r Mode (VM) and Service Mode (SM). Imaging data ob-
t whichtakeplaceintheupperlayersofEarth’satmosphere
s tained during SMruns are bias andflat-fieldcorrectedby
(Leinert et al. 1998 and references therein).
a
the pipeline and undergo a series of quality checks before
: With the exception of a very few cases, all night sky
v they are finally distributed to the users. Due to the high
i brightness surveys are executed with photoelectric de-
X number of imaging frames produced by this instrument
vices coupled to small telescopes and usually span a lim-
(morethan4500fromApril2000to September 2001)and
r ited number of nights (Benn & Ellison 1998), distributed
a thevarietyofscientificcaseswhichdriveit,itisclearthat
across several years. For these reasons, the data are usu-
a complete and systematic study of the night sky bright-
ally rather scanty and suffer from the inclusion of bright
ness can be performed only by means of a robust and
stars (V 13; see for example Walker 1988).
≥ automatic procedure, capable of identifying and rejecting
Nowadays, with the availability of large telescopes
all the cases which are not suitable for sky background
equipped with CCD imagers and the possibility of re-
measurements (e.g. large galaxies, crowded stellar fields
ducing the data via dedicated pipelines, the paradigm is
and so on).
changing and different approaches become feasible.
In this work we present and discuss the algorithm we
In this spirit, during the beginning of year 2000 we
have specifically designed for this purpose, while the re-
have started a project to monitor the UBVRI night sky
sults of the night sky brightness survey are reported in
brightness at ESO-Paranal Observatory (Chile) as part
Patat (2003).
of the quality control (QC) procedures implemented for
the FOcal Reducer/low dispersion Spectrograph (here- The paper is organizedasfollows.In Sec.2 we discuss
afterFORS1).Thismulti-modeopticalinstrument,which the problems connected with the sky background mea-
is mounted at the Cassegrainfocus of ESO-Antu/Melipal surement in digital images, while Sec. 3 deals with the
8.2m telescopes (Szeifert 2002), has two remotely ex- technique we have adopted to compute the mode of the
image intensity distribution. The algorithm we have de-
Send offprint requests to: F. Patat; e-mail: [email protected] vised to identify the presence of contaminating objects in
⋆ Partially based on observations collected at ESO-Paranal. the field and the tests we have performed on real FORS1
2 FerdinandoPatat: Skybackground computation
data arepresentedin Sec.4and5 respectively.Finally,in
Sec. 6 we summarize our conclusions.
2. Problems in estimating the sky background
Widely available programs for object detection and pho-
tometry like DAOPHOT (Stetson 1987) and Sextractor
(Bertin & Arnouts 1996) use the mode of the image in-
tensity distribution to estimate the local background or
to construct the background map of a two-dimensional
frame.Infact,besidesitsrobustness,themodeisastatis-
tically powerful estimator, being the most probable value
of the background brightness in a given region of the im-
age.
The use of the mode as the maximum-likelihood es-
timator for the sky background implicitly assumes that
the image has some regions which are not seriously con-
taminated by any astronomical object. Of course distant
and faint galaxies (and stars) will always be present, but
we will consider them as part of the global background
throughoutthispaper.Therefore,whenwetalkaboutcon-
Fig.1. Illustration of the effect on the intensity distri-
taminatingobjects,wemeanobjectswhicharedetectable
bution when a Poissonian noise is added to a stellar field
in the image, well above the background noise. There is
withaconstantbackgroundI .Theverticalscaleinthe
also another assumption that one has to make, i.e. that sky
lower plots is arbitrary, while the intensity I is counted
it is really improbable that the images to be analysed are
from the input I value and normalised to σ = I .
filled by a single extended object with a roughly constant sky p sky
The solidthincurveinthe lowerleftplotdepicts the con-
surface brightness. In fact this is the only case where one p
tributiontothefinalintensitydistributionbythevalueof
wouldhaveanoverallbackgroundincreasewithouthaving
F(I) at I =I +σ .
anysecondaryeffectsontheintensitydistributionfunction sky p
andinthatcasetheprocesswouldleadtoanoverestimate
of the sky background. In the case of FORS1 frames, es-
pecially with the commonly used SR collimator field of 6.′
8 6.′8, this assumption seems reasonable.
×
Due to the field of view of FORS1 and the wide va-
Of course there is always the possibility that none of
riety of scientific projects which are carried out by this
the sub-windows is clean enough to allow for a reliable
multi-modeinstrument,oneexpectstodealwithverydif-
measurement. A few examples are galactic stellar fields
ferentastronomicalobjects,whichwouldperturbinadif-
with heavily saturated stars, outer parts of globular clus-
ferent way the local sky background. Possible examples
ters,nearbyinteractinggalaxiesandclose-bycomets,just
are comets, clusters of galaxies, outskirts of big spirals,
to cite a few real examples we encountered during this
large ellipticals,diffuse nebulosities, crowdedstellar fields
analysis.
andsoon.Inallthesecasestherestillmightexistpartsof
the image which are suitable for a sky background mea- Thedifficultiesonefacesinestimatingthebackground
surement. For this reason it is clear that the analysis has are easily understood from the following considerations.
to be performed using sub-windows distributed on a grid If the images to be analysed were noiseless and the sky
acrosstheinputimage.Thechoiceofthesub-windowsize background constant and equal to I , the solution of
sky
hastobedoneinsuchawaythatthisisneithertoosmall, theproblemwouldbetrivial.Infact,inthiscase,thecor-
becauseinthatcasethefractionofuncontaminatedpixels responding intensity distribution function F(I) would be
might become statistically insignificant, nor too big, oth- 0 for I <I , would then suddenly peak at I =I and
sky sky
erwisetheprobabilityofincludinglargediffuseobjectsbe- finally show an extended tail, whose shape depends on
comeslarge.Afterrunningsometestswehaveseenthata the number and intensity of contaminating objects. This
300 300px sub-window (whichcorrespondsto 1 arcmin2 is illustrated in the left partof Fig. 1, where we have pre-
×
for the SRcollimator)givessatisfactoryresults.Since the sented an artificial stellar field and its intensity distribu-
guide probe of FORS1 is sometimes vignetting the outer tion. As usual, Nature behaves in a more subtle way and
parts of the images, we have decided to use the central due to the photon statistics (and marginally to detector
1800 1800pxregionofthedetectoronly.Thisallowsone read-out)the imageswearegoingto dealwith arealways
×
toanalysetheimagesina6 6sub-windowsgridincluding affected by noise. From a mathematical point of view it
9 104 px each. × is very easy to predict the noise effect on the intensity
×
FerdinandoPatat: Skybackground computation 3
distribution becomes more and more skewed, with a tail
appearing at the highest intensities. Clearly effect A and
B will respectively lead to overestimates of the sky back-
ground I and the noise, which in the case of an un-
sky
contaminated field is expected to be σ = I , accord-
p sky
ing to the Poissonian statistics. One possible solution to
p
this problem is achieved by reconstructing the intensity
distribution one would have if no contamination effects
were present. A good example of such a solution is repre-
sentedbytheasymmetricclippingalgorithmdevelopedby
Ratnatunga & Newell (1984), which computes the back-
ground distribution via an iterative clipping of the right
wing of the perturbed intensity distribution.
Due to the high number of available FORS1 frames
and the purpose of this work, we can afford a different
approach. Instead of attempting to reconstruct the un-
derlyingskybackgroundintensitydistributionincontam-
inated regions,we rather try to identify these regionsand
exclude them from all further calculations. For this pur-
posewehavedevisedasimpleandrobusttestthatcanbe
used to estimate the degree of contaminationin an image
Fig.2. Intensitydistributions forthreesimulatedimages and operated in an automatic way on large amounts of
obtainedinjecting 10,200and500stars(fromtopto bot- data.
tom) on a background I =5000 e− with a Poissonian
sky The first step in the application of our method is the
noise distribution. The solid curves are Gaussian profiles
mode estimate, which we describe in the next section.
centred on the distribution mode and having the σ ex-
pected for a Poissonian noise (see text for more details).
distribution. In fact, if n(I) is the noise distribution, the
observed intensity distribution is simply given by: 3. Mode estimate
I
Themode I ofadistributionf(I)isdefinedasthemost
f(I)= F(I z)n(I,z)dz (1)
h i
− probable value of I (see for example Lupton 1993) or, in
Z−∞
otherwords,the valueofI wheref(I)takesitsmaximum
i.e. the convolution F n of the signal F(I) with a
⊗ value. For moderately skewed distributions the mode can
variableresponsefunctionn(I,z),whichcanbeexpressed
as n(I,z) = 1/√2π z exp[−21(z2/(I −z)], provided that bmeeaanpp(rKoxeinmdaatllely&cSomtupaurtte1d9a9s7)h.IUin≃fo3rt×unmateedlyi,anth−e 2ob×-
the read-out noise is negligible. From these simple con-
servedintensitydistributionsinFORS1imagesoftenshow
siderations it is clear that what really matters for the
very extended tails, due to several effects like saturated
degradation of the resulting peak sharpness and the con-
stars,cosmic-rayevents and so on.While this tailusually
sequentuncertaintyintheestimateofitsoriginalposition
does not affect the mode, it does perturb the mean and
istheshapeoftheinputdistributionfunctionF(I)within
theaboveformulawouldleadtowrongresults.Apossible
5 I fromits peak, while the behaviour at higher in-
sky
∼ solution is given by the application of this approximation
tensities is unrelevant (see Fig. 1, right plots). With re-
p after an iterative clipping of the distribution around its
spect to these effects, it is quite instructive to look at
median, as it is done for example in Sextractor (Bertin &
the results of some numerical simulations. In Fig. 2 we
Arnouts1996).Inpracticeoneisforcedtousethismethod
have plotted the relative intensity distributions of three
whenonehastocomputethebackgroundmapofanimage
300 300px artificial frames with a fixed value of the sky
× usingsmallsub-windows.Infactinthatcasethestatistics
background (I =5000 electrons) on which we have ran-
sky
would be too poor to compute the mode in a reliable way
domlyinjected10,200and500starsrespectively.Wehave
directly using the distribution shape.
usedaMoffatprofileforthestarswithβ=4andaFWHM
of 5.3 px. In the case of FORS1 and SR collimator this This is not the case here, since we are more interested
would correspond to a field of 1 arcmin2 and a seeing of inanaveragevalueratherthaninamapofthebackground
1′′. The peak intensity of the stars was randomly gener- withinthesameimage.Forthisreasonwecanperformour
ated in the range 0-104 electrons. analysis in muchlargersub-windows so that we canbuild
Basicallythreeeffectsarevisibleasthenumberofstars up a very good signal-to-noise distribution function. This
grows: A) the mode I of the distribution steadily in- allows us to compute the mode just using its definition,
h i
creases; B) the distribution core width increases; C) the without any loss of generality and in a very robust way.
4 FerdinandoPatat: Skybackground computation
3.1. Introducing the Optimal Binning Technique
(OBT)
Since we are dealing with discrete distribution functions,
finding the mode implies that the data have to be binned
and the maximum of the distribution has to be found. To
reduce the noise due to the finite (and possibly large)bin
size∆I,wehavechosentorefinethedirectmodeestimate
using a quadratic interpolation on the modal bin and the
twoadjacentbins(seeforinstanceSpiegel1988).IfI and
m
f = f(I ) are the values of the intensity and intensity
m m
distribution in the modal bin and (I ,f ),(I ,f ) are the
r r l l
correspondingvaluesinthetwoadjacentbinsrespectively,
then the mode I is estimated as the value of I where
h i
the interpolatingparabolareachesits maximum.One can
show that this is given by the following expression:
f f ∆I
l r
I =I + − (2)
m
h i f 2f +f 2
l m r
−
This correctionhas the effect of moving the estimated
mode I within the whole bin ∆I according to the local
h i
shape ofthe distribution function. The modalbin bound-
Fig.3. Upperpanel:Signal-to-noiseratioSNR reached
m
aries are reached when f = f or f = f (provided
m l m r in the modal bin for the optimal value of ∆I as a func-
that f =f ). Numerical simulations haveshownthat the
l 6 r tion of npix. The solid line represents a linear fit to the
parabolicinterpolationisveryeffectiveinreducingtheer-
datawithn 90.Middlepanel:normalisedoptimalbin
pix
roronthemodeestimate(seebelow),atleastforthekind ≥
width. The solid line traces Eq. 8. Lower panel: optimal
of distributions considered here. For more general cases,
values of the RMS error on the mode. The dotted line is
the interpolation option will improve the mode accuracy a best fit to the data (ǫ n−0.75) while the dashed
onlyiftheasymmetryofthedistributionaroundthemode RMS ∝ pix
line traces the law ǫ n−1. In all panels each point
perturbs (fl fr)significantlyless thanthe averagevalue RMS ∝ pix
− is the result of 5000 simulations.
of the correction corresponding to an offset of the central
bin.
The choice of bin size is critical. If the bin ∆I is too This implies that if f(I) is the number of pixels falling in
small the resulting histogram will be too noisy to give a a given bin, the RMS uncertainty on f(I) is simply given
robust estimate of I . On the other hand, if ∆I is too by f(I) and hence we can define a signal-to-noise ratio
h i
large, then the mode estimate will be affected by a large SNR = f(I). In particular, we can introduce SNR ,
p m
uncertainty due to the small signal-to-noise in the two i.e.thesignal-to-noiseratioreachedinthe modalbin.Iff
p
outer bins. For this reason we need to find an optimal is the fraction of pixels falling in the modal bin, we have
value for ∆I which minimises the error on the mode. obviously
Forthispurposewehaverunaseriesofnumericalsim-
f =(SNR /n )2 (3)
ulations of uncontaminated windows with a known sky m pix
levelontopofwhichwehaveaddedaPoissoniannoiseand
Having this in mind, we can now have a look at the
atypicalread-outnoise(6electrons).Foreachofthe sim-
behaviourofSNR asafunctionofn whentheoptimal
ulatedframesthemodewasestimatedusingtheparabolic m pix
bin size is used to estimate the mode in our simulations.
interpolationdescribedabove.This hasbeen donefordif-
This is portrayed in the upper panel of Fig. 3, where one
ferent values of the bin size ∆I and the number of pixels
caneasilyseethatSNR dependsalmostlinearlyonn .
N = n2 included in each testing window. For each pair m pix
pix Intherange90 n 300,thebestfittothesimulated
(∆I,npix)wehaveperformed5000simulations.Whatone data (see Fig. 3≤, upppiexr≤panel) gives
sees is that for a given value of n , the RMS deviation
pix
of the estimated mode from the known input value Isky SNRm 6.2+0.48npix (4)
decreases as ∆I increases, reaches a minimum and then ≃
grows again, as expected from the above considerations For large values of n this relation can be approx-
pix
(see also Fig. 4). For this reason it makes sense to adopt imated as SNR n /2, which means that optimal
m pix
≈
the value of ∆I where the RMS error ǫ reaches its results are obtained when 25% of the pixels fall in the
RMS
∼
minimum as the optimum bin size for the mode estimate. modal bin.
Beforeweproceedwiththepresentationoftheresults, Having this result, it is easy to compute the optimal
we wantto discuss a few morepoints. The pixelcountsin fraction f of pixels by means of Eq. 3 and Eq. 4 (or its
the intensity distribution bins obey Poissonian statistics. approximate expression). Then, assuming the underlying
FerdinandoPatat: Skybackground computation 5
distribution f(I) to be a Gaussian with σ = σ , one can In the previous section we have seen that, for a suf-
p
derive the corresponding bin size in terms of σ as ficiently large value of n , the optimal bin size can be
p pix
expressed as:
∆I =2k(f)σ (5)
p
SNR2
∆I =2.6 m I (9)
wherekisimplicitlydefinedbythefollowingequation: sky
N
p
1 hIi+kσp (I I )2 Of course this requires to know the value of Isky, or
f = √2πσ exp − −2σh2i dI (6) at least to have a rough estimate of it. For this purpose,
p ZhIi−kσp (cid:20) p (cid:21) one can approximate it with the median of the distribu-
This equation can be solved numerically and the so- tion Imed, and SNRm can be computed using Eq. 4. For
lution fitted by a low order polynomial. For f 0.33 npix=300 we have SNRm 150.
≤ ≃
(n >60), the solution can be approximated rather ac- Once this guess value for the bin is computed, the
pix
curately by the following expression: histogram of the intensity distribution between I and
min
I is built and the mode is found. Since in the case of
max
k(f) 1.28f (7)
skewed distributions the median is only a rough approxi-
≃
mation of the mode and the distribution is definitely not
Finally,usingEqs.5and7onecaneasilycomputethe
Gaussian, the actual signal-to-noise ratio SNR′ in the
optimalbin width. Forn 90 we canapproximatethe m
pix
≥ modalbinismeasuredandfromthisvalueanewbinisre-
expression for the optimal bin as follows:
computedas∆I′ =∆I (SNR /SNR′ )2.Thehistogram
m m
∆I 6.2 2 is recalculated and a new mode value is estimated. This
2.6 0.48+ (8) procedure has a clear effect: if the distribution is skewed
σ ≃ n
p (cid:18) pix(cid:19) (and hence less peaked), then a larger bin size is required
while, for very large values of npix, the ratio between to achieve the same SNRm.
theoptimalbinsize∆I andσ approachesasymptotically The method has been tested using numerical simu-
p
the value 0.6. The result is shown in the central panel of lations of uncontaminated 300 300 px windows with a
×
Fig.3,wherewehavenormalisedthe optimalbintoσp to known sky level Isky on top of which we have added a
remove the dependency on the sky level. The comparison Poissoniannoiseandatypicalread-outnoise(6electrons).
betweenthepredictedvalues(solidline)andtheresultsof Foreachsimulatedfield,wehavethencomputedthemode
the simulations are in good agreement across all the npix I andthedeviationfromtheinputskybackgroundvalue
h i
explored range. defined as
Finally,theRMSdeviationfromtheinputvalueinour
I I
sky
simulations clearlydecreasesfor increasingvalues ofn , ǫ= h i− (10)
pix I
as shown in the lower panel of Fig. 3. For comparison we sky
have overplotted the n−1 law (dashed line) expected if InFig.4wepresentthe resultswe haveobtainedfora
pix
the error would scale proportionally to the overall signal- very low sky background (100 electrons). The RMS devi-
to-noise ratio. It is interesting to note that the RMS er- ation ǫ (solid line) reaches its minimum value (0.2%)
RMS
ror in the simulations decreases at a slower rate, approxi- for SNR =150, as expected, and so does the maximum
m
mately as n−0.75 (dotted line). We will come back to this error ǫ (dotted line). In the same figure we have also
pix max
pointinthenextsection,whendiscussingtheefficiencyof plotted the RMS error one obtains without using the re-
the method in reducing the error; for the time being, we finement given by Eq. 2 (ǫ′ , dashed line). This shows
RMS
onlynoticethatwithnpix=20theexpectedRMSerrorfor thatthe twomethodsgiveoptimalresultsattwodifferent
Isky=100electronsisabout1%,whichreducesto0.2%for bin sizes. In such conditions the interpolationreduces the
npix=300. RMSerrorbyaboutafactorof6.Forcomparisonwehave
also plotted the bin half width (long-dashed line), which
canbeassumedasanestimatorofthemaximumdeviation
3.2. Application of OBT to the real case
when the parabolic interpolation is not used.
As we have discussed in Sec. 2, the real data show a va- On the basis of these results one would tend to adopt
riety of contaminating effects, which tend to skew the in- SNR =150,butthere isacaveat thatwehavetokeepin
m
tensity distribution and affect in different ways the mode mind.Theseresultsarevalidforuncontaminateddistribu-
estimate. In that sectionwe have also mentionedthat the tions,wherethe signal-to-noiseratioSNR inthe modal
m
optimal window size for the mode estimate is n =300. bin is reached using the bin size given by Eq. 9. Since
pix
Both simulations and tests with real data show that with realdistributionsarecontaminated,alargerbinsizeisre-
such a number of pixels we can achieve RMS errors on quiredinordertoachievethesameSNR andthisinturn
m
the mode smaller than 1% for I 100 electrons, which generates larger errors. The simulations and extensive
sky
≥
we believe is a sufficient accuracy for our purposes. For tests with real data have shown that a good compromise
this reason, from now on we will concentrate on the spe- is reached using SNR =120 which, for mildly contam-
m
cific case of n =300 and discuss several aspects of the inated distributions gives maximum errors smaller than
pix
method application. 1% for a backgroundlevel of 100 electrons (see Tab. 1).
6 FerdinandoPatat: Skybackground computation
fitting Eq. 11 to the simulated data, we get N 1970.
eff
∼
Combining Eq. 9 and Eq. 11 one gets the following rela-
tionbetweenthebinsize∆I andthe expectedRMSerror
σ on the mode:
hIi
N ∆I
σ = (12)
hIi N 2.6SNR2
eff m
p
Substituting the proper values in this expression, we
obtain simply σ 0.05 ∆I, which means that for
hIi
≃
SNR =120 and moderately contaminated distributions,
m
theexpectedRMSerroronthemodeisoftheorderof5%
of the bin size. This is clearly shown in the lower panel
of Fig. 5, where we have plotted the measured error as a
function of the relativebin size derivedfromthe same set
of simulations previously described. The match between
the expected RMS (dotted line) and the observed RMS
(dashed line) is fairly good.
In both panels of Fig. 5 we have plotted the maxi-
mum errors encountered during the simulations. As one
cansee,theyarewellconfinedwithinthe5σ level(dashed
line), which was computed using the RMS error given
Fig.4. EstimatedRMSmodeerrors(solidline)fromsim-
by the simulations. More precisely, maximum errors lie
ulated 300 300 px un-contaminated windows as a func-
× with a good approximation on the 4σ level. From Fig. 5
tion of SNR . The maximum error is traced with a dot-
m we can conclude that the expected RMS errors on the
ted line, while the dashed line indicates the RMS error
mode estimate are below 0.3% at all sky backgroundlev-
when the parabolic interpolation is not used. Finally, the
els larger than 100 electrons, while maximum errors are
long-dashed line indicates the percentage bin half width.
always smaller than 1%.
The values for each SNR level are the result of 5000
m Theintroductionofartificialstarshastheeffectofsys-
simulations.
tematically increasing the mode of the distribution with
respecttotherealvalueoftheskybackground(seeSec.2).
To see how the error behaves as a function of the sky
For this reason, in the general case of contaminated dis-
background,wehaveperformedanothersetofsimilarsim- tributions, we can talk about two different errors. While
ulations,wherewehaveadoptedSNR =120andwehave
m theformerisarandommeasurementerrorintrinsictothe
variedthe input skylevelfrom102 to 104 electrons.A to-
adopted method and to the quality of the data,the latter
talof5000simulationsperlevelwereexecuted.Theresults is a systematic error which depends on the intensity dis-
are presented in Fig. 5 and summarized in Tab. 1. tributionofthecontaminatingobjects.Aswehavesaidin
In the upper panel of Fig. 5 we have plotted the per- Sec.2,wewanttouseonlythecaseswherethesystematic
centage errors ǫ as a function of the backgroundlevel. As error is smaller than some threshold value. The approach
expected, the percentage RMS error, traced with a solid to this problem is describedin Sec. 5, while here we focus
line, decreases for increasing values of I , following very only on the randomerrorsrelatedto the way the mode is
h i
well the usual law for the standard error of the average: estimated.
Toevaluatetheeffectofcontaminationonthemethod
I error, we have performed several sets of simulations in-
σ = h i (11)
hIi sNeff jecting a fixed number of stars N∗ with random positions
on a sky background of intensity I . The stars intensi-
sky
with the difference that Neff is smaller than the total ties I were uniformly generated in the range 0 < I
∗ ∗
number N of pixels. For this reason, Neff can be consid- (Isat Isky), where Isat is the detector’s numerical satu≤-
eredastheeffectivenumberofdatapointsdrawnfromthe ration−level (I 106,000 electrons for the FORS1 de-
sat
originalsetonewouldhavetousetogetthesameRMSer- tector, high gain)∼. Finally, a Poissonian noise was added
ror when adopting the averageas backgroundestimator1. to the artificial 300 300 px frames to simulate the pho-
As we have seen at the end of Sec. 3.1, the RMS error tonshotnoiseandt×he modewasmeasuredwith the OBT
scales as n−pi0x.75, and hence we have Neff = N0.75 when using SNRm=120.Numerical tests using a more realistic
theoptimalSNRmisadopted.Sincewehavechosentouse intensity distribution, drawn from observed star counts,
SNRm=120, we expect Neff to be even smaller. In fact, showthatverysimilarcontaminationeffectsareachieved.
1 Of course this is true only if the distribution is not con- The only difference is that one needs to generate much
taminated, because otherwise the average gives much larger more artificial stars in the non-uniform case, and this
systematic errors makes it numerically less efficient.
FerdinandoPatat: Skybackground computation 7
Table 1. Estimated errors from numerical simulations of 300 300 px uncontaminated windows for three different
×
values of the sky background. Maximum errors are indicated within parenthesis. A total of 5000 simulations per sky
level and SNR value were performed.
m
SNRm f (%) nb(±σ) ∆I/hIi(%) RMS error (%)
102 103 104 102 103 104
120 16.0 5 4.7 1.3 0.4 0.25 (1.00) 0.07 (0.27) 0.02 (0.11)
150 25.0 3 7.5 2.1 0.7 0.17 (0.60) 0.04 (0.16) 0.02 (0.06)
4. The ∆-test
Oncethemode I oftheskybackgrounddistributionand
h i
the corresponding error σ have been computed, one is
hIi
left with the task of recognizing and rejecting contami-
nated sub-windows, which is the topic of this section.
As we have already mentioned, in an uncontaminated
image the intensity distribution f(I) obeys the photon
shotnoisedistributionandhencetheexpectedRMSnoise
is given by σ = I . Since the left wing of f(I) is
p
h i
much less disturbed than the right wing (cfr. Fig. 2), we
p
can estimate the underlying background noise using only
the N pixels for which I I . As a robust estimator
l
≤ h i
we have chosen to use the Median Absolute Deviation
(MAD) from the mode I applied to the left wing of
h i
the distribution:
MAD =median I I I I (13)
l
{| −h i| ≤h i}
Wenotethattheuseofthestandarddeviationinstead
of MAD has the effect of over-estimating the noise in the
real cases. In fact, the lower tail of the intensity distri-
Fig.5. Relative errors on the mode estimate as a func- bution is always affected by defects like bad pixels and
tion of the sky background (upper panel) and the rel- possible vignetting. While this fact has a rather strong
ative bin size (lower panel). The circles represent the influence on the standard deviation, it leaves the MAD
measured RMS deviations, while the dotted line indi- unperturbed, which is much less sensitive to the presence
cates the expected RMS (see Eq. 12). In both panels, the of outliers.
crossesindicate the maximumdeviationencountereddur- One can easily show that for a Gaussian distribution
ing the simulations, while the dashed line traces the 5σ theratiobetweenthestandarddeviationandtheMADis
level. Simulations were performed using SNR =120 and 1.483 (Huber 1981). Hence, to have an estimator which
m
∼
RON=6 electrons. hasthesamemeaningofthestandarddeviationinthenor-
malcase2,wedefineσ =1.483MAD .Thisallowsoneto
l l
comparedirectly σ andσ which, foranuncontaminated
p l
distribution,shouldbe approximatelythe same.Toquan-
tifypossibledeviationsfromthePoissonianbehaviour,we
As expected, the simulations show that the system- introduce the parameter ∆:
atic error grows with N , while the random error keeps
∗
obeyingto Eq.12witha goodapproximation,atleastfor ∆= σl2−RON2−σp (14)
N∗ <200 and I 500 electrons. More precisely, in the p σp
h i ≥
case of N = 200, the simulations give σ 0.08 ∆I
∗ hIi
in the range 102 I 104 electrons. We≃can safely which equals 0 in the Poissonian case, whilst it gets
concludethattheR≤MhSier≤rorsintroducedbythemodede- largerandlargerasthedeviationfromthePoissoniandis-
termination method we have described in this section are tribution grows. The idea is to try to estimate the un-
smaller than 0.3% for I 5 102-104 electrons. Finally, knownsystematicerrorfrom∆,whichisameasurablepa-
the RMSerrorcanbechonis≥erva×tivelyassumedtobe 8%of rameter. To illustrate this concept we use a realexample,
the bin size ∆I in the same intensity range, which corre- 2 Infact,sincethecentralvaluehIioftheintensitydistribu-
spondstoNeff 1000.Thesevalueshavebeenusedinthe tionswearedealingwithisalwayslargerthanseveralhundred
∼
remaining sections of this work and for all sky brightness counts, the noise Poissonian distribution can be very well ap-
measurements discussed in Patat (2003). proximated with a Gaussian with σ= hIi.
p
8 FerdinandoPatat: Skybackground computation
Fig.6. An image of NGC 3521 obtained by FORS1 on
Fig.7. Intensity distribution for the sub-window placed
04-04-2000 (R, 30 seconds, standard resolution collima-
on a outer region of NGC 3521 (see Fig. 6). The vertical
tor). The box indicates the test sub-window used for the
linesareplacedatthedistributionmode(dashed)andthe
exampledescribedinthetext.Itincludes300 300pixels,
which correspond to 1′ 1′ on the sky. × realskybackground(dashed-dotted).Thelatterwasman-
× uallymeasuredintheupperleftcorneroftheinputimage,
in an object-free area. The solid and dotted lines plotted
onthe leftside ofthe distributionaretwoGaussianswith
where we have performed the analysis on a 300 300 px
sub-window centred on the outer parts of the lar×ge spiral σ = σp and σ = σl respectively. The figure insert shows
an intensity contour plot of the region.
galaxy NGC 3521,as shown in Fig. 6. The corresponding
intensity distribution is shown in Fig. 7, where the mode
(indicatedbytheverticaldashedline)wascomputedwith i.e. 102, 103 and 104 electrons. The circles represent ǫ ,
theOptimalBinningTechnique(OBT)wehavedescribed whilethedashedanddottedlinesindicatetwodifferenthesi-
in Sec. 3. timatesoftherandomerrorσ ;oneistheRMSdeviation
hIi
A superficial inspection of the input image shows al- directly computed from the simulated data and the other
ready that this sub-window is strongly contaminated by is estimated using Eq. 11. While the systematic error is
diffuse emission from the galaxy’s spiral arms. The mode clearly due to the presence of contaminating objects, the
of the intensity distribution is in fact 34% higher than randomerrorσ isrelatedtotheaccuracyofthemethod
∼ hIi
the background level manually measured on an object- one adopts to estimate the mode I (see Sec. 3).
h i
free areaof the same image close to the upper left corner. A clear result emerges from Fig. 8: the ∆ parameter
However, the presence of strong fluctuations within the is effective in estimating the systematic error and hence
sub-windows produces a significant increase in the distri- gives a statistically significantcriterionto decide whether
bution core width, which becomes 80% larger than the a sub-window can be considered as unperturbed or not.
∼
expected photon shot noise (∆=0.78). This suggests that In order to have a simple description of the behaviour of
∆ may provide a tool to estimate the systematic error ǫ as a function of ∆, we have performed a linear fitting
induced by contaminating sources and this, coupled with hofithe simulated results in the range 0 ∆ 10 (see
some threshold criterion, would allow us to recognise and Fig. 8, continuous lines). Actually ǫ ten≤ds to≤bend for
disregard the critical cases. To study this possibility, we larger values of ∆ and this has theheffiect that our linear
haveexecuteda seriesofsimulationsofcontaminateddis- fittingslightlyoverestimatesthevalueof ǫ .Thisisnota
tributions of the same type of those described in Sec. 3.2 problem, since in this sense the criterionhwie are going to
for several values of the input sky background intensity establishisjustmorerestrictive.Asexpected,theslopeof
Isky. For each simulation one can then compute the er- the ∆, ǫ relation is inversely proportional to the signal-
ror ǫ of the mode estimate (see Eq. 10) and measure ∆. to-noisehriatiooftheskybackground.Thebestfitgivesthe
The resulting range in the ∆ parameter can be binned to following result:
some suitable value and the corresponding average error
ǫ computed. In Fig. 8 we have plotted the results one 3.3
h i ǫ = (∆+0.96) (15)
obtains for three different values of the sky background, h i I
sky
p
FerdinandoPatat: Skybackground computation 9
intensity. For a typical value of ǫ =1%, ∆ is 1.1,
max max
7.7 and 28.4% for sky levels of 102, 103 and 104 electrons
respectively.
Once all sub-windows in an input image have been
analysed with the criterion we have just discussed, a first
sky brightness guess can be obtained using the median of
all selected values. This has the effect of excluding cases
like those producedbyocculting masks insertedin the fo-
calplane.Infact,toavoidstrongsaturationeffects,FORS
instruments allow the user to place a number of movable
blades on specific positions of the focal plane. The height
of these blades is about 20′′, which correspond to 100 px
whentheSRcollimatorisused.Whentheoccultedregions
are larger than the adopted sub-window size, there is a
chance that such areas pass the ∆-test. Since the counts
in those regions are very low, the median always removes
them successfully, provided that the occulted fraction of
the fieldofviewisnotlargerthan50%,aconditionwhich
is always fulfilled.
After this selection is done, the only background fluc-
tuations which are expected to be left in the remaining
sub-windowsareduetoanonperfectflat-fielding.Infact,
Fig.8. Errors on the estimated sky background I as
sky
as we have mentioned in the introduction, the use of twi-
a function of ∆ (see text) from numerical simulations
light flats introduces large scale gradients,which produce
for three different sky background levels (102, 103 and
maximumpeak-to-peakdeviationsof6%fromperfectflat-
104 electrons). The continuous lines represent linear least
ness (see also next section). For this reason, we have de-
squares fits in the range 0 ∆ 10, while the circles are
≤ ≤ cidedtooperateafurtherrefinement,choosingonlythose
the average (systematic) errors. Each point is the result
n sub-windows which deviate less than 3% from the me-
of5000simulations.The dashedand dottedlines indicate g
dianvalue.Thefinalestimateofthe backgroundintensity
the RMS random errors computed from the simulations
I is eventually obtained computing the weighted mean
and Eq. 11 (N = 1000) respectively. The horizontal sky
eff I = ng I w / ng w , (w = 1/σ2 ). To allow
dotted line is placed at ǫmax=1%. sky j=1h ij j j=1 j j hIi
for a statistically significant result, we have imposed that
P P
n 5inourautomaticprocedures.Ifthisconditionisnot
g
≥
where I is given in electrons and ∆,ǫ are expressed met, then the input frame is considered as unsuitable for
sky
as percentage values. The fact that for ∆=0 ǫ =0 is skybackgroundmeasuresandrejected.Thishastheeffect
h i 6
due to the following effect. When the sky background is ofoperatingafirstroughfilteringontheinputdata.Tobe
contaminated by faint stars, whose maximum intensity is conservative,n isloggedtogetherwiththeotherrelevant
g
comparable to I , the mode and the width of the in- parameters,andafurtherandmorerestrictiveselectionis
sky
tensity distribution increase in such a way that ∆ tends always possible.
p
to be small, while the effective error is not. Intensive tests with realdata,as discussed in the next
We can now use the total maximum error ǫ to section, have shown that the set of selection criteria we
max
fix a limit on ∆ below which we can consider a given havediscussedmakethemethodreliable,robustandsuit-
intensity distributions as practically unperturbed. Since able to be implemented in a fully automatic procedure.
the distribution of random errors around the system-
atic deviation is Gaussian, we can conservatively assume
5. Testing the method on real data
ǫ ǫ +3σ . Using Eq. 15 for ∆ one gets
max hIi
≃|h i|
The method we have just outlined has been tested on
ǫ
∆ =0.9 I max σ 0.96 (16) a sample of 4678 FORS1 reduced frames, obtained with
max h i 3 − hIi − different filters between April 1, 2000 and September 30,
p h i
2001. For each input frame the results relative to all 36
whereσ isgivenbytheadoptedmodeestimator.As
hIi
sub-windows have been logged and this has allowed us to
it is shown in Sec. 3, this can be approximated by Eq. 11
build-up a sample with more than 168,000 entries, each
with N =1000, which gives the following expression:
eff
of which has been flagged according to the results of the
∆-test. The basic results are shown in Fig. 9, where we
I
∆max = h i ǫmax 1.9 (17) have plotted the measured noise (corrected for the read-
3.3 −
p out noise) as a function of the expected Poissoniannoise.
Thecriticalvaluefor∆dependsofcourseontheaccu- The solid line traces the locus where the two noises have
racyonewantstoreachintheestimateofthebackground thesamevalue,whilethedashedoneindicatesthelimiton
10 FerdinandoPatat: Skybackground computation
Fig.9. Comparisonbetween the expected and measured Fig.10. Maximum deviations from the sky background
noises on a sample of more than 4600 FORS1 images. weightedmeaninallframeswhichpassedthe∆-test.Solid
Each point refers to a single sub-window. The solid line line indicates the data from all sub-windows, the dotted
indicates the locus where the expected noise equals the line refers to the selected windows and the dashed line
measured one, while the dashed line traces the limit set corresponds to the cases where n =36. The upper right
g
by Eq. 17. The dashed-dotted line indicates the expected insert shows the cumulative functions with the same line
global noise generated by the photon statistics and the codings.
flat-fielding correction for a typical FORS1 configuration
(seetextformoredetails).Theupperpanelshowsthefrac-
ally negligible with respect to the noise present in the
tionofsub-windowswhichpassedthe∆-testasafunction
science frame to be corrected. When the background ex-
of the expected noise.
posure level reaches high values, this is no longer true
and the noise added by the flat-fielding process becomes
the measured noise imposed by Eq. 17. All sub-windows significant. The expected global noise is in fact given by
lying below the latter line wouldbe selectedfor sky back- σ2 = σ2 +I2 /(N I ), where I is the background
p sky F F sky
ground estimates. level in the input science frame. This law gives a fair re-
The fraction of selected windows is plotted in the up- production of the observed behaviour, as it is shown in
per panel of the same figure, again as a function of the Fig.9(dashed-dottedline),wherewehaveusedthetypical
expected noise. As one can see, for noise values of less values of I and N quoted above. The few data points
F F
than 20 electrons (or background intensities smaller than lying below this line are generated by observations ob-
400 electrons), this fraction is very low. For larger in- tained with the low gain mode ( 0.3 ADU electron−1),
∼ ∼
tensities, itreachesa roughlyconstantvalue of80%.This for which I is abouta factor of2 largerthan in the high
F
means that basically all frames with I 400 electrons gainmode( 0.6ADUelectron−1),whichisalsothestan-
h i ≤ ∼
will be rejected. In the adopted test sample, this however dardforFORS1imaging.Inprinciple,whencomputing∆,
accounts for 11% of total number of frames only. From one could correct the measured noise for the flat-fielding
the test sample we can conclude that, on average, 86.5% effect. Inpractice,atanexposurelevelof5 104 electrons
×
of FORS1 frames are suitable for sky-background mea- thecorrectiononthemeasurednoiseisoftheorderof20%
surements, 95% of which have n 12. and hence a very small number of sub-windows which are
g
≥
An interesting feature visible in Fig. 9 is the system- rejected by the ∆ criterion would move to the safe re-
atic trend shown by the minimum measured noise. In gion (see Fig. 9). Furthermore, 90% of all sub-windows
fact, this tends to deviate more and more from the ex- included in the test sample have a background intensity
pectedPoissoniannoiseforbackgroundvalueslargerthan smallerthan2 104 electrons.Atthis level,the correction
104 electrons. This effect is explained by the following amounts to ab×out 10% only. For this reason we have de-
∼
considerations. FORS1 master sky flat-fields are usually cidedto ignorethe flat-fieldingeffect whenevaluating the
obtained by the combination of N =3 twilight sky flats, deviation from the pure Poissoniannoise.
F
which have typical exposure levels of I 3.5 104 elec- Aswehavementioned(seealsoSec.3formoredetails),
F
∼ ×
trons. The resulting noise on the combined frame is usu- the maximum formal errors on the mode estimate within