Table Of Content

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 operator: 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.