Astronomical Image Subtraction by Cross-Convolution

0.27 MB·English

by  Fang Yuan

Draftversion February2,2008
PreprinttypesetusingLATEXstyleemulateapjv.10/09/06
ASTRONOMICAL IMAGE SUBTRACTION BY CROSS-CONVOLUTION
Fang Yuan1 and Carl W. Akerlof1
Draft versionFebruary 2, 2008
ABSTRACT
In recent years, there has been a proliferation of wide-field sky surveys to search for a variety of
transientobjects. Using relatively shortfocal lengths, the optics of these systems produce undersam-
pled stellar images often marred by a variety of aberrations. As participants in such activities, we
have developed a new algorithm for image subtraction that no longer requires high quality reference
8 imagesfor comparison. The computationalefficiency is comparablewith similar procedurescurrently
0 in use. The general technique is cross-convolution: two convolution kernels are generated to make a
0 test image and a reference image separately transform to match as closely as possible. In analogy to
2 the optimization technique for generating smoothing splines, the inclusion of an RMS width penalty

n term constrains the diffusion of stellar images. In addition, by evaluating the convolution kernels
a on uniformly spaced subimages across the total area, these routines can accomodate point spread
J functions that vary considerably across the focal plane.

2 Subject headings: techniques: photometric, methods: statistical

]
h
1. INTRODUCTION ThischallengewasaddressedbyRobertQuimbywhen
p
hebecameagraduatestudentattheUniversityofTexas.
- The advent of low-noise megapixel electronic image
o AsanundergradattheUniversityofCaliforniaatBerke-
sensors,cheapfastcomputers,andterabytedatastorage
r ley, he had worked extensively with the Supernova Cos-
t systems has enabled searches for rare astrophysicalphe-
s nomenathatwouldotherwisebeeffectivelyundetectable. mologyProject(SCP)andwasquitefamiliarwiththeSN
a discoveryprocess. Quimby adapted the SCP image sub-
ExamplesincludethediscoveriesofMACHOs,transitsof
[
extra-solarplanets,andvastlygreaternumbersofsuper- tractioncode(Perlmutter et al.1999)foruseasthebasic
1 novae. Commontoalloftheseeffortsisaneedtorepeat- tool for finding 30 SNe over a period of two years from
v edly image large portions of the sky to uncover rare and observations with the University of Texas 30% alloca-
6 tion of ROTSE-IIIb time at the McDonald Observatory.
subtle changes of brightness. This forces the observerto
3 Among those discoveries are SN 2005ap (Quimby et al.
contend with images taken under less than ideal condi-
3 2007) and SN 2006gy (Smith et al. 2007) which appear
tions such as poor weather and crowded star fields. To
0 to be the intrinsically brightest SNe ever identified.
solvetheproblemofcomparingimagestakenwithdiffer-
1. ent seeing conditions, a number of research groups have In view of the evident success of the Texas Super-
0 developed image subtraction algorithms which compen- novaSearch(TSS)(Quimby2006),ourgroupattheUni-
8 sateforimageblurringeffectspriortoimagedifferencing. versity of Michigan probed the image subtraction prob-
0 lem with the goal of applying this to the considerably
In this paper, we describe a variation of this technique
: more extensive image data available to the entire suite
v which treats pairs of images in a symmetric fashion, re-
of ROTSE-III telescopes. The originalhope of using the
i ducing the requirements of first obtaining an ideal refer-
X SCP code was abandoned following the realization that
ence image. The software code is relatively simple and
the programwould not be made freely available. We at-
r has been made freely available.
a tempted to adopt the ISIS image subtraction package2
Followingourinitialdiscoveryofpromptopticalradia-
but were discouraged by the initial results. The signifi-
tionfromagamma-rayburstin1999,ourROTSEcollab-
cantundersamplingofROTSE-IIIstellarimagescoupled
oration set out to construct a set of four identical wide-
withasymmetricpointspreadfunctionsacrosstheimage
fieldtelescopes(Akerlof et al.2003)toexplorethesephe-
plane created a severe challenge for making clean sub-
nomena more deeply. The resulting instruments, called
tractions. These issues are not satisfactorily addressed
ROTSE-III, were installed in Australia, Texas, Namibia
by the algorithms described by Alard & Lupton (1998)
and Turkey in 2003 and 2004. By explicitly choosing a
and Alard (2000) for two reasons: (1) we do not always
fast (f/1.9) optical system and short focal length (850
have the luxury of a substantially higher quality refer-
mm), we provided the option to search for possible or-
ence image and (2), the point spread functions (PSFs)
phan GRB afterglowswithout unduly compromisingthe
are often approximately elliptical with the axes oriented
potential for mapping GRB optical light curves at early
at any angle in the image plane. For a variety of rea-
times. Although we alwaysintended to use these instru-
sons, the ROTSE-III PSFs can vary with temperature
mentsforgenericastrophysicalopticaltransientsearches,
and telescope orientation. Thus, the possibility that one
itwasrecognizedthattheplatescale(3.3”perpixel)was
cansimplyconvolveanewimagetoanidealreferenceim-
too coarse for easy identification of a supernova embed-
age is not always viable. With this in mind, we sought
ded in a normal galaxy.
to develop a more symmetric algorithm that would be
1University of Michigan, Randall Laboratory of Physics, 450 robust enough to handle less pristine observations. We
Church St., Ann Arbor, MI, 48109-1040, [email protected],
[email protected] 2 http://www2.iap.fr/users/alard/package.html
2 Yuan & Akerlof
should emphasize that the aim of this project is primar- poorly represent the trend of the data. The solution to
ily for the reliable identification of transients in a very this problem is to add a curvature penalty term to the
large database, not precision photometry. least squares residuals so that a trade-off is reached be-
tweenadequately fitting the data andinserting unneces-
2. MATHEMATICALMETHOD sarilycomplexbehaviorintothesmoothedinterpolation.
The basic technique for image subtraction presented The coefficient that scales the curvature term is a mea-
by Alardand Lupton depends onfinding a suitable PSF sureofthestiffnessofthespline. Thereexistsanelegant
smearing kernel, K(u,v) that when convoluted with the method called cross-validationthat determines the stiff-
referenceimage,R(x,y), generatesatransformedimage, ness parameter from the standard deviation errors for
R∗(x,y), that can be compared on a pixel-by-pixelbasis each data point.
with a new test image of lesser quality, T(x,y): For image subtraction, the degradation of the signal-
to-noise ratio is proportional to the effective number of
R(x,y)⊗K(u,v)=R∗(x,y) (1) pixels that are summed by the convolutionkernel. Since
each pixel has an associated variance, σ2 , the variance
pix
The kernel, K(u,v), is constructed by a linear superpo- for a signal diffused over N pixels will be N σ2 .
sition of basis functions of the form: pix pix pix
Assuming Gaussian distributions, we can estimate N
pix
fn,p,q(u,v)=upvqe−(u2+v2)/2σn2 (2) ftrioonm, tNhe wi∼dth4πof(σth2e eff+ectσi2ve),stwelhlaerrepoσint spirseatdhefubnac--
pix PSF K PSF
Alard and Lupton recommend a three-fold ensemble of sic stellar PSF width (in pixels) and σK is the diffusive
terms with σ values spanning a 9-fold range. These width of the convolution, K. If σ2 is evaluated by the
n K
functionsarepoorchoicestosynthesizeanellipticalPSF normal formula, σ2 = K r2 where K are the kernel
K P i i i
at an arbitrary angle with respect to the imager sensor element amplitudes, the sum can be deceptively small
axes although they will be satisfactory for PSFs with when the values for K alternate in sign. Noting that
i
pcloasnedtqomazuimstubtehadlestyermmmineetrdy.byTahde hspoceccifiomc vpaalruiseosnfsorwσitnh, hbKyii1 ∝ σK1K22,r4thwehvicahlueequfoarllyσK2pencaanlizbesebboetthteproesisttiivmeaatnedd
the characteristic PSFs associated with the particular 2π P i i
negative contributions to the kernel elements. Although
instrument in use. The set of amplitudes for the basis
the scaling behavior of the penalty coefficient is under-
functions is computed by the least squares technique to stood in terms of the image size, σ2 , and σ2 , we
minimize the pixel-by-pixel differences between R∗ and pix PSF
havenot investigatedwhether there is anelegantwayto
T.
evaluatethisquantityanalogouslytothecross-validation
From this starting point, we decided to symmetrize
technique for splines.
theAlard-Luptonprocedurebycreatingtwoconvolution
operators so that: 3. COMPUTATIONALMETHODS
The image preprocessing that we require is similar to
R(x,y)⊗K (u,v)≈T(x,y)⊗K (u,v) (3)
R T thatdescribedbyAlardandLupton. Flat-fieldedimages
areprocessedby SExtractor(Bertin & Arnouts 1996)to
Inthelimitthatthereferenceimageissubstantiallybet-
ter than the test image, the K operatorsmears R with create object lists with precise stellar coordinates. The
the point spread function chaRracteristic of T while K IDL routine3, POLYWARP, is used to warp the new test
T
image to overlay stellar objects in the reference image
will be essentially equal to the identity operator so that
itseffectonT willbenegligible. Underthese conditions, as closely as possible. A valid pixel map is generated
to avoid pixels close to the image perimeter and screen
the computation becomes functionally equivalent to the
against saturated values, etc. At this point, the fun-
procedure adopted by Alard and Lupton. However, in
damental image subtraction code is invoked as an IDL
general, the addition of a second convolution operator
routine which first performs image flux normalizationto
injects new mathematical degrees of freedom that must
be constrained. The most obvious is that K and K equalize the mean values of the two images under com-
R T
can be multiplied by an arbitrary constant without vi- parison.
olating image convolution equality. This can be conve- The most basic choice that the user must make is the
niently resolved by demanding that one or both kernels representation of the convolution kernels. We have re-
be flux-conserving, ie: strictedthemton×narrayswithnodd. Thispermitsa
simple representationforthe convolutionidentity opera-
tor: K = 1 while all other elements of K are
XK(u,v)=1 (4) [n/2],[n/2] i,j
zero. For the ROTSE-III images, n = 9 appears to pro-
Thesecond,morecomplex,problemarisesfromthedif- videmorethanadequatecoverageofstellarpointspread
fusionofastellarimageifK andK approximatebroad functions under the worst conditions.
R T
Gaussian distributions. The convoluted image equality The values for the convolutionkernel elements are de-
will be maintained but the signal-to-noise ratio of the rived from the difference image:
subtracted image will drastically diminish. The solution
to this was found by analogy to a similar problem in D(x,y)={R(x,y)⊗K (u,v)}−{T(x,y)⊗K (u,v)}
R T
the application of smoothing splines to drawing curves (5)
through data with errors. In the latter case, one could Invokingthe criterionthat D(x,y)2 shouldbe a mini-
P
trivially create a spline curve that ran exactly through mum subject to the requirements for K andK gener-
R T
each and every data point (as long as the abscissas are
distinct). Such a curve would appear very wiggily and 3 ITTVisualInformationSolutions,ITTIndustries,Inc.
IMAGE SUBTRACTION BY CROSS-CONVOLUTION 3
ates2(n2−1)linearequationsviatheusualleastsquares what Robert Quimby had obtained using the SCP code
procedure to solve for the independent coefficients for as adapted for his Texas Supernova Search.
K (u,v) and K (u,v) after imposing the kernel unitar- MostoftheimagesubtractioncodewaswritteninIDL
R T
ityconstraints. Asdescribedearlier,theseequationswill withtheexceptionoftheevaluationoftheregressionma-
notprovideuniquesolutionsforK andK becausethe trix. Since this is the core of the computational burden
R T
effective width of the two convolution kernels can still andIDL is notparticularly efficient in handling the nec-
be radiallyscaledwithoutsubstantiallyaffectingthe dif- essary array indexing, this portion was coded in C and
ference image, D. Thus, the quantity to be minimized linked to the rest of the IDL programs using the IDL
must include a penalty term for radially diffusing the standardexternalcallinginterface. Acrucialdetail,par-
convolutedimagesanyfurtherthannecessary. Following ticularly important for the ROTSE-III telescopes, is the
fromearlierremarks,thisfigure-of-meritfunctioncanbe variation of the stellar PSF across the image plane. To
represented as: accomodate this problem, each 2045×2049 image was
subdividedinto36sub-imagesofroughlyequalsize. The
setof50kernelamplitudecoefficientswascalculated,one
Q=XD(x,y)2+λX(u2+v2)2{KR(u,v)2+KT(u,v)2} by one, for each of these sub-images. Cross-convolved
(6) images were obtained by bilinear interpolation for ev-
where λ is a constant selected to balance the contribu- ery pixel using the four nearestneighbor coefficient sets.
tionsofthetwocompetingerrorterms. Fromthediscus- Theunitarityoftheconvolutionkernelsisguaranteedby
sion given above, the value for λ should scale as: the linearity ofthe interpolationmethod with respectto
the gridded coefficient reference values. Although this
(σ2 +σ2)
λ=2πN R T λ′ (7) sounds somewhat complicated, the calculation was ex-
image σ2 tremely efficient.
PSF
Another significant concern is the estimation of the
where N is the totalnumber of pixels in the image,
image backgroundskyintensity. Initially,wereliedonSExtrac-
σ2 and σ2 are the pixel amplitude variances, σ is
thRe characTteristic stellar PSF width and λ′ is a coPnSstFant tor to remove this background before subtraction. How-
ever,when applying our code to images containing large
of order unity.
galaxies such as M31 and M33, we realized that these
Withtwo9×9convolutionkernels,thenumberoffree
backgrounds are poorly estimated around the cores of
parameters is 160 and the size of the regression matrix
bright galaxies. The solution we adopted only removes
becomes problematic. The main concern is that if the
the background difference between the two images in-
images are essentially featureless (ie. no stars), the ma-
stead of the individual background for each image sepa-
trix elements become indistinguishable and the inverse
rately. After the images are scaled so that stellar fluxes
matrix will be singular. To avoid these effects as well as
match,askydifferencemapisgeneratedbyfirstperform-
variousothercomputationalissues,abinaryvaluedmask
ingpixelbypixelsubtraction. Thelowfrequencyspatial
array is created to eliminate sampling around the image
variationofthisdifferenceimageisobtainedbyaprocess
perimeter, saturated pixels and all featureless areas not
similartotheoneusedbySExtractor. Thedifferenceim-
associatedwithstellarobjectsasdeterminedbySExtrac-
age is divided into 32×32 pixel subimages and median
tor. This approach was quite successful: the degree of
pixelvaluesarerecursivelyevaluated,subjecttothecon-
singularityof the regressionmatrix was determined dur-
straint that pixels with 3σ excursions from the median
ing the inversion process using the IDL singular value
are ignored. Saturated pixels are also excluded from the
decomposition routines SVDC and SVSOL, codes derived
mediancomputation. Theresultantslowlyvaryingback-
from Numerical Recipes in C(Press et al. 1992).
groundissubtractedfromoneoftheinputimagesbefore
For our ROTSE project, computational efficiency is
invokingthecross-convolutionalgorithm. Theremaining
criticalbecausewetypicallyacquire400imagespernight
common non-zero background doesn’t affect the estima-
with each telescope and these must be reduced in situ
tion of the convolution kernels and the final subtracted
comfortably within 24 hours. It was easily verified that
image will be background-free.
most of the image subtractioncalculations were devoted
A comparison of the results of the cross-convolution
to computing the convolution kernels regression matrix
and the single convolution algorithms is shown in Fig-
described above. Examination of the two-dimensional
ure 2. In the limiting case where the PSFs of the refer-
structureofthe kernelsshowedthattheamplitudes near
ence image are azimuthally symmetric, the two meth-
the edges of the 9×9 arrays were always small and sug-
ods should produce rather similar results. However,
gested that the representation could be significantly re-
whenthatconditionisnotsatisfied,thecross-convolution
ducedfrom81values to 25by assuminga mapping from
method is more appropriate.
a reduced number of wavelet functions. Thus, each con-
For anyone wishing to employ the cross-
volutionkernelwasrepresentedbyalinearsuperposition:
convolution technique described in this paper, the
K(u,v)=XAiBi(u,v) source code can be downloaded from the URL,
http://hdl.handle.net/2027.42/57484, which points to
with the basis functions, B , chosen as discrete approx- a set of files stored in Deep Blue, the University of
i
imations to bicubic spline functions with characteristic Michigan institutional repository. Note: This Web
widths of 1, 3 and 2 pixels, centered as shown in Fig- site is currently inaccessible. Please contact the
2
ure 1. Using this technique shrankthe regressionmatrix authors directly.
from 160×160 to 48×48 with a consequent reduction
in processing time of about anorder of magnitude. This
broughtthecomputationthroughputtovaluessimilarto 4. OPERATIONALEXPERIENCE
4 Yuan & Akerlof
Subtraction of a 2045× 2049 ROTSE image from a higher for those corresponding to stellar objects. This
reference frame using the method described above takes later criterion is intended to suppress variable stars.
approximately 4 minutes with a 2.0 GHz personal com- In the 5-month period to date, the pipeline has iden-
puter. Ifthesamereferenceimageisusedmultipletimes, tified all 13 reported supernovae that lie within our
it needs to be convolvedwith the base kernels just once, searched fields. One of these initially escaped but was
savingcomputationaltime. Subtractionsofthreetypical detected following modification of the mask size to pro-
ROTSE images from the same reference frame takes ∼ vide better performance during bad seeing. Also due
10 minutes on the same processor. It should be noted to our early inexperience, two of these SNe in relatively
that the memory allocation for the process, mainly for bright galaxies were initially missed during hand scan-
storingthebasekernelconvolutedimages,scaleswiththe ning. In addition, the pipeline detected 7 novae in the
size of the image and number of kernel basis sets. For fields of M31 and M33. Two novae rather close to the
our choice of 25 kernel base vectors and a 2045×2049 center of M31 were missed before the background eval-
image size, ∼ 1 GB memory is required. This is not a uation problem was addressed as described in section 3.
serious handicap for a modern desktop computer. In terms of transient recoveryefficiency, both real-world
Since August 2007, a supernova search pipeline us- andlimitedMonteCarlocomparisonsshowthatoursub-
ing this subtraction code has been running on images tractioncodeiscomparabletothemodifiedversionofthe
takenby the ROTSE-IIIbtelescope. Selected fields with Supernova Cosmology Project search code employed by
nearby rich clusters and a high density of known galax- the Texas Supernova Search.
ies are monitored on a daily basis, weather permitting,
to a typical limiting magnitude of 18.5. For each field, 5. SUMMARY
two sets of four 60-sec exposures (20-sec for fields with
The algorithm described in this paper can be adapted
bright target galaxies) are taken with a 30-minute ca-
for a wide variety of photometric searches for transient
dence. Following the method developed by the Texas
objects. Its performance appears to be at least as good
Supernova Search,the imagesfor each4-exposureepoch
as other codes presently in use. Since the method is de-
are co-added as are the total 8 images for the night. All
signed to handle images with significantly varying qual-
three co-additions are subtracted by the same reference
ity,itshouldremaineffectivewhen alternativeprograms
image. ThedifferenceimagesareprocessedthroughSEx-
may fail.
tractortofindresidualobjects. Torejectfalsedetections
duetobadpixels,cosmicrays,asteroidsandsubtraction
noise, further filtering is applied. The signal-to-noisera- 6. ACKNOWLEDGEMENTS
tio of a candidate has to be above5 inthe nightly 8-fold The authors especially thank Robert Quimby for his
sumand2.5forthesumofsingleepoch. Thepositionsof advice, suggestions and encouragement for implement-
a candidate in 3 subtractions must match to within one ing this image analysis code and validating the results.
pixelfordetectionsabove15SNRand1.5pixelsforthose We also appreciate the contributions of a number of
with SNR below 15. The FWHM of the candidate has ROTSEcollaborators,particularlyJamesAretakis,Tim-
tolie withinthe rangeofone pixelandtwicethe median othyMcKay,EliRykoffandHeatherSwan. Thisresearch
FWHM for stars in the convolved reference image. Fi- was supported by NASA grant NNG-04WC41G. F. Y.
nally,minimumfluxchangecutsareappliedwithalower wassupportedbyNASASwift GuestInvestigatorgrants
thresholdfordetectionsembeddedinknowngalaxiesand NNG-06GI90G and NNX-07AF02G.
REFERENCES
Akerlof,C.W.,etal. 2003,PASP,115,132 Quimby,R.M. 2006,PhDthesis,UniversityofTexas,USA
Alard,C.&Lupton,R.H. 1998,ApJ,503,325 Quimby,R.M.,etal. 2007,ApJ,668,L99
Alard,C. 2000,A&AS,144,509 Smith,N.,etal. 2007, ApJ,666,1116
Bertin,E.&Arnouts,S. 1996,A&AS,117,393
Perlmutter,S.,Aldering,G.,Goldhaber,G.,etal. 1999,ApJ,517,
565
Press,W.H.,Teukolsky,S.A.,Vetterling,W.T.,&Flannery,B.P.
1992 Numerical recipes in C: The art of scientific computing.
CambridgeUniversityPress,Cambridge,2dedition
IMAGE SUBTRACTION BY CROSS-CONVOLUTION 5
Fig. 1.— Diagram of the location of the 25 bicubic B-splinesused to construct the convolution kernels. The 9 circles, 8squares and 8
hexagons markthecenters oftheB-splinemaximawithwidthsof1, 3 and2pixelsrespectively. Thedashedlinesindicatethe9×9grid
2
oftheunderlyingconvolutionkernels.
6 Yuan & Akerlof
Fig. 2.— Comparison of image subtractions using the cross-convolution method described in this paper and the single convolution
methoddescribedbyAlardandLuptonandimplementedintheISIScode. TheinitialimageswereobtainedbytheROTSE-IIIbtelescope
atMcDonaldObservatory. Shownhereare260×260pixelsubframescentered onα=16:50:02.21,δ=+23:46:32.88,coveringafield
of0.235◦×0.235◦. Todemonstratetheresults,threeartificial“variable”starswereaddedtothetestimage(a)andthereferenceimage(b)
withPSFs appropriatelymatched totheir respective fields. Thelocations areshownbyblack arrows. The subtracted imageobtained by
cross-convolutionisdepictedin(c)andtheAlard-Luptonresultsareshownin(d). Thebrightstarnearthelowerrightcorneroftheimages
hasbeenreplacedwithauniformgraylevelsinceneither subtractiontechniquecanextractusefulinformationfromsaturatedpixels.