ConvPhot: A Profile-Matching Algorithm for Precision Photometry C. De Santis, A. Grazian, A. Fontana, P. Santini INAF - Osservatorio Astronomico Roma Via Frascati, 33 – 00040 Monte Porzio 7 Catone (RM), Italy 0 0 2 n Abstract a J 9 We describe in this paper a new, public software for accurate “PSF-matched” 1 multiband photometry for images of different resolution and depth, that we have v namedConvPhot,ofwhichweanalyseperformancesandlimitations.Itisdesignedto 2 work when a high resolution image is available to identify and extract the objects, 3 2 and colours or variations in luminosity are to be measured in another image of 1 lower resolution but comparable depth. To maximise the usability of this software, 0 we explicitly usetheoutputs of thepopularSExtractor code, that is used to extract 7 0 all objects from the high resolution “detection” image. The technique adopted by h/ the code is essentially to convolve each object to the PSF of the lower resolution p “measure” image, and to obtain the flux of each object by a global χ2 minimisation - o on such measure image. We remark that no a priori assumption is done on the r shape of the objects. In this paper we provide a full description of the algorithm, t s a discussion of the possible systematic effects involved and the results of a set of a : simulations andvalidation tests that wehave performedon real as well as simulated v i images. The source code of ConvPhot,written in C language underthe GNU Public X License, is released worldwide. r a Key words: methods: data analysis – techniques: image processing – techniques: photometric 1 Introduction The availability of efficient imagers, operating with good imaging quality over a large range of wavelengths, has opened a new era in astronomy. Multi- Email address: desantis,grazian,fontana,[email protected](C. De Santis, A. Grazian, A. Fontana, P. Santini). Preprint submitted to Elsevier Science 4 February 2008 wavelength imaging surveys have been executed or planned, and are provid- ing major breakthroughs in many fields of modern astronomy. These surveys often collect images of different quality and depth, typically resulting from the combined effort of ground and space–based facilities. In this context, the difficulties originating in the analysis of these often inhomogeneous imaging databases have hampered the proper exploitation of these data sets, especially in the field of faint, high redshift galaxies. On the one side, images at varying wavelengths may provide a surprisingly different glimpse of the Universe, with objects fading or emerging from the background. On the other side, the image resolution is usually not constant overthewavelengths,duetodifferentinstrumentcharacteristics.Atypicalcase is the combination of high resolution HST images with lower resolution images obtained by ground based telescopes or by the Spitzer Space telescope. In the latter case, the blending among the objects in the lower resolution images often prevents a full exploitation of the multicolour informations contained in the data. The difficulties involved in the analysis of this kind of data has led to the development of several techniques. The first emphasis was based on refine- ments of the usual detection algorithms (Szalay, Connolly & Szokoly 1999). The algorithm here discussed, that is designed to work especially for faint galaxies, allows instead to accurately measure colours in relatively crowded fields, making full use of the spatial and morphological information contained in the highest quality images. This approach has already been proposed and used in previous works (Fernandez–Soto et al., 1999; Papovich et al., 2001; Grazian et al., 2006): here, we discuss a specific implementation of the soft- ware code, named ConvPhot, that we developed and make publicly available to analyse data with inhomogeneous image quality. Although we focus in the following on multi-wavelength observations, i.e. the case where different im- ages of the same portion of sky are available in different bands, the approach followed here can be adopted also for variability studies, where images in the same band but taken at different epochs are used. In this case, the highest qualityimagecanbeusedtoidentifytheobjects,andthemagnitudevariations of all the objects can be measured with no systematic effects. The plan of the paper is the following: in Sect. 2, we outline the technique adopted; in Sect. 3, we discuss more exhaustively the algorithm, of which we provide more details. We also describe the problematics involving the deter- mination of the isophotal area and the optimisations we adopted. In Sect. 4 we comment on possible issues of ConvPhot, taking into account the cases of blended sources in the detection image. In Sect. 5 we describe the systematics that may be cause of inappropriate results from the algorithm. In Sect. 6, we discuss the validation tests that we performed on simulated as well as on real images, discussing in particular the usage of this software and the future 2 prospects for improving ConvPhot. A brief summary and conclusions of the paper are given in Sect. 7. 2 The basic technique The technique that we discuss here has been described and adopted for the first time by the Stony-Brook group to optimise the analysis of the J, H and K images of the Hubble Deep Field North (HDF-N). The method is described in Fernandez–Soto et al. (1999, hereafter FSLY99) and the catalog obtained has been used in several scientific analysis of the HDF-N, both by the Stony- Brook group (Lanzetta et al., 1999; Phillipps et al., 2000), as well as by other groups, including our own (Fontana et al., 2000; Poli et al., 2001). The same method has been adopted by (Papovich et al., 2001; Dickinson et al., 2003), to deal with the similar problems existing in the HDFs data set and, recently, to derive a photometric catalog of the GOODS fields in the 24 micron band of MIPS (PSF is 5 arcseconds) using the 3.6 micron band of IRAC (1.6 arcsec of resolution) to deblend (Chary et al., in prep.) the MIPS sources1. Although in both these cases the method has been applied to the combination between HST and ground–based, near–IR data, the procedure is much more general and can be applied to any combination of data sets, provided that the following assumptions are satisfied: 1. A high resolution image (hereafter named “detection image”) is available, that is used to detect objects and isolate their area; this image should be deep enough to allow all the sources in the “measure image” to be detected. Ideally, such image should be well sampled and allow a proper resolution for most sources. 2. Colours (or temporal variations) are to be measured in a lower resolution image (hereafter named “measure image”); 3. The PSF is accurately estimated in both images, and a convolution kernel hasbeenobtainedtosmooththe“detection image”tothePSFofthe“measure image”. Conceptually, the method is quite straightforward, and can be better under- stood by looking at Fig. 1. In the upper panel, we plot the case of two objects that are clearly detected in the “detection image”, but severely blended in the “measure” one. The technique described here can remind that adopted by DOPHOT (Schechter et al., 1993), a software for PSF fitting in crowded stellar fields. There are otherwise some differences, since ConvPhot is thought especially for galaxies or extended objects in general, and the minimisation of the model image to the measure image is done simultaneously, reducing the 1 http://data.spitzer.caltech.edu/popular/goods/Documents/goods dr3.html 3 possibility of wrong fits in severely blended objects, an usual fact when the fit is carried out into different steps. Indeed, the ConvPhot software needs a model of the convolution kernel in input, while in DOPHOT the PSF is fitted analytically. The procedure followed by ConvPhot consists of the following steps: a) Object identification is done relying the parameters and area obtained by a SExtractor (Bertin & Arnouts, 1996) run on the “detection image”: in prac- tice, we use the “Segmentation” image produced by SExtractor, that is ob- tained from the isophotal area. We extract small images centered around each object, on which all the subsequent operations will be performed. Such small images (dynamically stored in the computer’s memory) are named “thumb- nails” in the following. Details are described in Section 3.1, including the treatment of border effects. b) Since the isophotal area is typically an underestimate of the actual object size, such that a fraction of the object flux is lost in the tails outside the contour of the last isophote, we have developed a specific software (named dilate) to expand the SExtractor area of an amount f proportional to DIL the object size. This procedure is described in Section 3.2 c) Another required input is a convolution kernel that converts the detection image to the same resolution of the “measure” image. Although this is not providedwithinConvPhot,ashortdescriptionoftherequirements andpossible methods is given in Section 3.3. d) The background of each object (both in the “detection” and in the “mea- sure” images) is computed, as described in Section 3.4; e) Each object, along with its segmentation image, is (individually) smoothed, and then normalised to unit total flux: we refer to it as the “model profiles” of the objects, as it represents the expected shape of the object in the “measure” image (Section 3.5). f) Finally, the intensity of all “model” objects are simultaneously adjusted in order tomatchtheintensity oftheobjects inthe“measure”image.This global minimisation provides the final normalisation (i.e. intensity) of each object in the “model” image, and hence the fundamental output of ConvPhot. This is described in Section 3.6 In this scaling procedure there are as many free parameters (that we named 4 Fig. 1. A schematic representation of the ConvPhot algorithm. a). Two objects are clearly detected and separated in the high resolution detection image (blue, solid-thin line). The same two objects are blended in the low resolution measure image (red, solid-thick line) and have quite different colours. b). The two objects are isolated in the high resolution detection image and are individually smoothed to the PSF of the measure image, to obtain the “model” images. c). The intensity of each object is scaled to match the global profile of the measure image. The scaling factors are found with a global χ2 minimisation over the object areas. F , the fitted flux of the i-th galaxy) as the number of objects in the detec- i tion image, that are typically thousands. The free parameters F are com- i puted with a χ2 minimisation over all the pixels of the images, and all ob- jects are fitted simultaneously to take into account the effects of blending between nearby objects. Conceptually, this approach is somewhat similar to the line fitting of complex absorption metal systems in high redshift quasars (Fontana & Ballester, 1995). Although in practical cases the number of free parameters can be quite large, the resulting linear system is very sparse (FSLY99), since non null terms represent only the rare overlapping/blended sources, and the minimisation can be performed in a quite efficient way by using standard numerical techniques. 5 As can be appreciated from the example plotted in Fig. 1, the main advan- tage of the method is that it relies on the accurate spatial and morphological information contained in the “detection” image to measure colours in rela- tively crowded fields, even in the case that the colours of blended objects are markedly different. Still, the method relies on a few assumptions that must be well understood and taken into account. First of all, it is assumed that the objects have no measurable proper motion. This is not a particular concern in small, deep extragalactic fields like HDF or GOODS, but might be important in other applications. A more general concern is that morphology and positions of the objects should not change significantly between the two bandwidths. Also, the depth and central bandpass of the “detection” image must ensure that most of the objects detected in the “measure” image are contained in the input catalog. Finally, the objects should be well separated on the detection image, although they may be blended in the measure image. In practice, it is unlikely that all these conditions are satisfied in real cases. In the case ofthe match between ACS andgroundbased Ks images, forinstance, very red objects may be detected in the Ks band with no counterpart in the optical images, and some morphological change is expected due to the increasing contribution of the bulge in the near–IR bands. This mis-match may be particularly significant for high redshift sources, where the knots of star formation in the UV rest frame are redshifted to the optical wavelengths used as model images. Also, in the case that the pixel-size of the “measure” image (i.e. ISAAC or VIMOS) is much larger than that of the “input” one (ACS for example), the actual limitations due to intrinsic inaccuracy in image aligning may lead to systematic errors. We will show below that these effects may be minimised or corrected with reasonable accuracy: at this purpose, we have included in the code several options and fine–tuning parameters to minimise the systematics involved. 3 The ConvPhot algorithm Inthisparagraph,weillustrateindetailthealgorithmadoptedfortheConvPhot software. Particular emphasis is given to resolve critical issues in this profile matching algorithm, such as a correct estimation of the background, an un- biased reproduction of the object profile and a fast and robust minimisation technique. We analyse in detail these critical points. 6 Fig. 2. This figure shows how ConvPhot defines the overlap region between the detection and the measure image in case of different sizes and offsets, provided by the input parameters CPHXOFF and CPHYOFF. The objects whosesegmentation is partly inside the overlap region are flagged as bound. 3.1 Input images and catalog As mentioned above, ConvPhot relies on the output of the popular SExtrac- tor code to detect and deblend the objects and to define the object position and main morphological parameters in the detection image. In the following, we shall also use the SExtractor naming convention to specify several image characteristics. We remind that a “segmentation” image is an image where background is set to 0 and all the pixels assigned to a given objects are filled with the catalog number of the object. The user must first obtain a reliable SExtractor catalog of the detection image prior to running ConvPhot. In practice, the required output from SExtractor consists of the SEGMENTATION image and of a catalog with the center and boundary coordinates of the objects in pixel units and the value of the background, as shown in Tab. 1. 7 Table 1 Example of input catalog for ConvPhot. NUMBER X Y XMIN XMAX YMIN YMAX BACK D BACK M 1 2319.71 216.86 2134 2507 34 384 1.350E-05 0.0 2 8849.97 79.61 8820 8881 32 124 -1.300E-06 0.0 3 5269.25 16.99 5245 5293 1 39 -6.372E-06 0.0 4 3596.99 21.14 3581 3613 5 37 -2.710E-05 0.0 5 4878.07 11.14 4852 4905 1 37 -6.634E-06 0.0 NUMBER is the identification number of the source in the segmentation image created by SExtractor, X, Y, XMIN, XMAX, YMIN, YMAX are the center and boundarycoordinatesoftheobjects(inSExtractortheyareindicatedasX IMAGE, Y IMAGE, XMIN IMAGE, XMAX IMAGE, YMIN IMAGE, YMAX IMAGE, respectively). BACK D is the value of the background for the detection image, while BACK M is the background for the measure image. It is possible to provide every object with its local background estimate, for example that produced by SExtractor. This example is taken from the GOODS-MUSIC catalog produced using the z band provided by ACS in the GOODS South field (Grazian et al., 2006). In addition to this, the user must provide the detection image and the measure withitsMAP RMS,i.e.animagecontaining theabsoluter.m.s.ofthemeasure image, in the same units. It is important to remind that the measure image may be shifted and/or have a different size with respect to the detection one: it is only required that the images are aligned (i.e. have been rebinned to the same pixel size and be astrometrically registered; for problems related to the rebinning of the measure image see Sect. 6). This offsetting option may be useful when dealing with the follow–up of large mosaics (i.e. GOODS, COSMOS). Fig. 2 shows the case of a measure image having a different size and shifted respect to the detection one. In this case, the relative offset (in pixels) between the two images must be provided by the input parameters CPHXOFF and CPHYOFF. The measure image and its RMS must have the same size. As shown in Fig. 2, ConvPhot works on objects whose segmentation lies in the overlapping region between the detection and the measure image. In the detection image, objects which partly lie inside the overlap region are flagged as bound and will be processed only if the fraction between the flux inside the overlap region and the total flux in the detection image is equal or larger than a user-defined parameter f . flux 3.2 Small segmentation for objects The “object segmentation” defines the area (i.e. the set of pixels) over which each object is distributed. This area, which is different in the measure and 8 in the detection image because of the typical different image quality (PSF), is used to extract the object profile in the detection image and to define the fitting region in the measure one. To establish the area of each object in the detection image, ConvPhot relies on the “segmentation” image produced by SExtractor. In this image, each pixel contained within the isophotal area of an object is filled with the relevant object ID, and can be used to reconstruct the shape of the object itself. This estimate, however, is not very robust since the size of the isophotal area de- pends on the isophotal threshold adopted in the detection image, and in any case misses a significant fraction of the flux of faint objects. This effect is particularly dramatic as the galaxy becomes either very faint or very large, as in space based images like those provided by HST (WFPC2, ACS),withresulting lowsurfacebrightness oftheexternalpixels. To overcome this problem, we have developed and distributed a dedicated software (named dilate), that expand the segmentation produced by SExtractor, though pre- venting merging between neighbouring objects and preserving the original shape of the segmentation. The segmentation of each object is individually smoothed with an adaptive kernel which preserves its shape and is then renor- malised to the original number value. This software allows to fix the dilate factor f for the magnification of the segmentation, that is defined as the DIL ratio between the output and input isophotal area. While this procedure is useful for relatively bright objects, of which the increase of size is doubled for an f = 4,it provides toosmall an enlargement for the faintest objects, that DIL are often detected over very few pixels (this is particularly true in the case of faintsmallgalaxies observed withHST,asinthecaseoftheGOODSdata).At this purpose, we allow to define a minimum area m for the segmentation, AREA such that objects smaller thanthis area (even after dilation) are forcedto have this area. A minimum value for the area of extended segmentation is useful, for instance, in the ACS images, where extended but low surface brightness objects are detected by SExtractor only due to the bright nucleus and their isophotal area is limited only to the brighter knots. These values should be tuned accurately, because too small values for f or m may cause a DIL AREA distortion of the profile for the detection image, which alters the photometry in the measure image, while too large values for f or m are not useful DIL AREA since enhance the noise in the model profile. Typical values for f are 3-4, DIL in order to double the area of the original segmentation image. The m AREA parameter is set in order to match 2-3 times the typical size of galaxies in the field: for deep imaging surveys with HST, faint galaxies have an half light radius of Rh = 0.2−0.3 arcsec, which translates to a minimum area of around 800 pixels, if the pixel scale of ACS (0.03 arcsec) is used. Moreover, ConvPhot also compute the so-called “detection magnitude”, that is the magnitude of the object in the detection image for the dilated segmentation. This quantity should be compared with the total detection magnitude of the input catalog 9 in order to correct the model magnitude for the fraction of flux missed by the limited segmentation. In addition, the possibility of restricting the fit to the central region of the objects (as described in Sect 3.5) further suppresses the effects of tails in the fit itself. 3.3 PSF matching A key process within ConvPhot is the smoothing of the “detection” high res- olution image to the PSF of the lower resolution “measure” one. The original feature of the method is that this step is not performed on the global im- age, but rather on each individual object, after it has been extracted from the “detection” image, in order to prevent the confusion to to blending in low resolution images. Such a smoothing is performed by applying a convolution kernel, that is the filter required to transform the PSF of the detection image into the PSF of the “measure” image. Filtering is performed on the thumbnails of the extracted objects, as well as on the corresponding segmentation and auxiliary images, whenever required. The derivation of the convolution kernel is not done by ConvPhot and must be provided by the users. Several techniques have been developed to derive an ac- curate convolution kernel, in particular those by Alard & Lupton (1998) and Alard (2000), and we advice to use such sophisticated techniques whenever possible. In the case of the GOODS–MUSIC sample, where we had to obtain more than 70 independent kernels (one for each original image of the J, H and K mosaics), we have adopted a simpler method, based on analysis in the the Fourier space. Such method can be implemented and automatized with standard astronomical tools for image analysis. We summarise here the basic recipe to compute it. First of all, we derive the PSF of the detection and measure image, P and P respectively, e.g. by summing up several stars in 1 2 the field, in order to gain in Signal to Noise ratio. These two PSFs must be normalised, in order to have total flux equal to unity. The convolution kernel K is, by definition, P (x,y) = K(x,y)⊗P (x,y) . (1) 2 1 The derivation of the exact shape of the convolution kernel is done in the Fourier space: taking ℵ = FT(K), ℘ = TF(P ) and ℘ = TF(P ), the 1 1 2 2 10