multidrizzle.minmed

25 """ Create a median array, rejecting the highest pixel and computing the lowest valid pixel after mask application""" 26 27 """ 28 # In this case we want to calculate two things: 29 # 1) the median array, rejecting the highest pixel (thus running 30 # imcombine with nlow=0, nhigh=1, nkeep=1, using the masks) 31 # 2) the lowest valid pixel after applying the masks (thus running 32 # imcombine with nlow=0, nhigh=3, nkeep=1, using the masks) 33 # 34 # We also calculate the sum of the weight files (to produce the total 35 # effective exposure time for each pixel). 36 # 37 # The total effective background in the final image is calculated as follows: 38 # - convert background for each input image to counts/s (divide by exptime) 39 # - multiply this value by the weight image, to obtain the effective background 40 # counts (in DN) for each pixel, for each image 41 # - Add these images together, to obtain the total effective background 42 # for the combined image. 43 # 44 # Once we've made these two files, then calculate the SNR based on the 45 # median-pixel image, and compare with the minimum. 46 """ 47 48 """ In this version of the mimmed algorithm we assume that the units of all input data is electons."""

49 - def __init__(self, 50 imageList, # list of input data to be combined. 51 weightImageList, # list of input data weight images to be combined. 52 readnoiseList, # list of readnoise values to use for the input images. 53 exposureTimeList, # list of exposure times to use for the input images. 54 backgroundValueList, # list of image background values to use for the input images 55 weightMaskList= None, # list of imput data weight masks to use for pixel rejection. 56 combine_grow = 1, # Radius (pixels) for neighbor rejection 57 combine_nsigma1 = 4, # Significance for accepting minimum instead of median 58 combine_nsigma2 = 3 # Significance for accepting minimum instead of median 59 60 ):

61 62 # Define input variables 63 self.__imageList = imageList 64 self.__weightImageList = weightImageList 65 self.__weightMaskList = weightMaskList 66 self.__exposureTimeList = exposureTimeList 67 self.__readnoiseList = readnoiseList 68 self.__backgroundValueList = backgroundValueList 69 self.__numberOfImages = len( self.__imageList) 70 self.__combine_grow = combine_grow 71 self.__combine_nsigma1 = combine_nsigma1 72 self.__combine_nsigma2 = combine_nsigma2 73 74 75 # Create a different median image based upon the number of images in the input list. 76 __median_file = np.zeros(self.__imageList[0].shape,dtype=self.__imageList[0].dtype) 77 if (self.__numberOfImages == 2): 78 __tmp = numCombine(self.__imageList,numarrayMaskList=self.__weightMaskList, 79 combinationType="mean",nlow=0,nhigh=0, 80 nkeep=1,upper=None,lower=None) 81 __median_file = __tmp.combArrObj 82 else: 83 # The value of NHIGH=1 will cause problems when there is only 1 valid 84 # unmasked input image for that pixel due to a difference in behavior 85 # between 'numcombine' and 'iraf.imcombine'. 86 # This value may need to be adjusted on the fly based on the number of 87 # inputs and the number of masked values/pixel. 88 # 89 __tmp = numCombine(self.__imageList,numarrayMaskList=self.__weightMaskList, 90 combinationType="median",nlow=0,nhigh=1, 91 nkeep=1,upper=None,lower=None) 92 __median_file = __tmp.combArrObj 93 94 # The following section of code will address the problem caused by having 95 # a value of nhigh = 1. This will behave in a way similar to the way the 96 # IRAF task IMCOMBINE behaves. In order to accomplish this, the following 97 # procedure will be followed: 98 # 1) The input masks will be summed. 99 # 2) The science data will be summed. 100 # 3) In the locations of the summed mask where the sum is 1 less than the 101 # the total number of images, the value of that location in the summed 102 # sceince image will be used to replace the existing value in the 103 # existing __median_file. 104 # 105 # This procuedure is being used to prevent too much data from being thrown 106 # out of the image. Take for example the case of 3 input images. In two 107 # of the images the pixel locations have been masked out. Now, if nhigh 108 # is applied there will be no value to use for that position. However, 109 # if this new procedure is used that value in the resulting images will 110 # be the value that was rejected by the nhigh rejection step. 111 # 112 113 # We need to make certain that "bad" pixels in the sci data are set to 0. That way, 114 # when the sci images are summed, the value of the sum will only come from the "good" 115 # pixels. 116 tmpList = [] 117 for image in xrange(len(self.__imageList)): 118 tmp =np.where(self.__weightMaskList[image] == 1, 0, self.__imageList[image]) 119 tmpList.append(tmp) 120 121 # Sum the mask files 122 maskSum = self.__sumImages(self.__weightMaskList) 123 # Sum the science images 124 sciSum = self.__sumImages(tmpList) 125 del(tmpList) 126 # Use the summed sci image values in locations where the maskSum indicates 127 # that there is only 1 good pixel to use. The value will be used in the 128 # __median_file image 129 __median_file = np.where(maskSum == self.__numberOfImages-1,sciSum,__median_file) 130 131 # Sum the weightMaskList elements 132 __maskSum = self.__sumImages(self.__weightMaskList) 133 134 # Create the minimum image from the stack of input images. 135 # Find the maximum pixel value for the image stack. 136 _maxValue = -1e+9 137 for image in self.__imageList: 138 _newMax = image.max() 139 if (_newMax > _maxValue): 140 _maxValue = _newMax 141 142 # For each image, set pixels masked as "bad" to the "super-maximum" value. 143 for image in xrange(len(self.__imageList)): 144 self.__imageList[image] = np.where(self.__weightMaskList[image] == 1,_maxValue+1,self.__imageList[image]) 145 146 # Call numcombine throwing out the highest N - 1 pixels. 147 __tmp = numCombine(self.__imageList,numarrayMaskList=None, 148 combinationType="median",nlow=0,nhigh=self.__numberOfImages-1, 149 nkeep=1,upper=None,lower=None) 150 __minimum_file = __tmp.combArrObj 151 # Reset any pixl at _maxValue + 1 to 0. 152 __minimum_file = np.where(__maskSum == self.__numberOfImages, 0, __minimum_file) 153 154 # Scale the weight images by the background values and add them to the bk 155 __backgroundFileList = [] 156 for image in xrange(len(self.__weightImageList)): 157 __tmp = self.__weightImageList[image] * (self.__backgroundValueList[image]/(self.__exposureTimeList[image])) 158 __backgroundFileList.append(__tmp) 159 160 # Create an image of the total effective background (in DN) per pixel: 161 # (which is the sum of all the background-scaled weight files) 162 # 163 __bkgd_file = self.__sumImages(__backgroundFileList) 164 del(__backgroundFileList) 165 166 # 167 # Scale the weight mask images by the square of the readnoise values 168 # 169 __readnoiseFileList = [] 170 for image in xrange(len(self.__weightMaskList)): 171 __tmp = np.logical_not(self.__weightMaskList[image]) * (self.__readnoiseList[image] * self.__readnoiseList[image]) 172 __readnoiseFileList.append(__tmp) 173 174 # Create an image of the total readnoise**2 per pixel: 175 # (which is the sum of all the input readnoise values) 176 # 177 __readnoise_file = self.__sumImages(__readnoiseFileList) 178 del(__readnoiseFileList) 179 180 # Create an image of the total effective exposure time per pixel: 181 # (which is simply the sum of all the drizzle output weight files) 182 # 183 __weight_file = self.__sumImages(self.__weightImageList) 184 185 186 # Scale up both the median and minimum arrays by the total effective exposure time 187 # per pixel. 188 # 189 __minimum_file_weighted = __minimum_file * __weight_file 190 __median_file_weighted = __median_file * __weight_file 191 del(__weight_file) 192 193 # Calculate the 1-sigma r.m.s.: 194 # variance = median_electrons + bkgd_electrons + readnoise**2 195 # rms = sqrt(variance) 196 # This image has units of electrons. 197 # 198 __rms_file = np.sqrt(__median_file_weighted + __bkgd_file + __readnoise_file) 199 200 del __bkgd_file, __readnoise_file 201 # For the median array, calculate the n-sigma lower threshold to the array 202 # and incorporate that into the pixel values. 203 # 204 __median_rms_file = __median_file_weighted - (__rms_file * self.__combine_nsigma1) 205 206 if self.__combine_grow != 0: 207 # 208 # Do a more sophisticated rejection: For all cases where the minimum pixel will 209 # be accepted instead of the median, set a lower threshold for that pixel and the 210 # ones around it (ie become less conservative in rejecting the median). This is 211 # because in cases of triple-incidence cosmic rays, quite often the low-lying 212 # outliers of the CRs can influence the median for the initial relatively high 213 # value of sigma, so a lower threshold must be used to mnake sure that the minimum 214 # is selected. 215 # 216 # This is done as follows: 217 # 1) make an image which is zero everywhere except where the minimum will be accepted 218 # 2) box-car smooth this image, to make these regions grow. 219 # 3) In the file "median_rms_file_electrons", replace these pixels 220 # by median - combine_nsigma2 * rms 221 # 222 # Then use this image in the final replacement, in the same way as for the 223 # case where this option is not selected. 224 225 __minimum_flag_file = np.where(np.less(__minimum_file_weighted,__median_rms_file),1,0) 226 227 # The box size value must be an integer. This is not a problem since __combine_grow should always 228 # be an integer type. The combine_grow column in the MDRIZTAB should also be an integer type. 229 __boxsize = int(2 * self.__combine_grow + 1) 230 __boxshape = (__boxsize,__boxsize) 231 __minimum_grow_file = np.zeros(self.__imageList[0].shape,dtype=self.__imageList[0].dtype) 232 233 234 # If the boxcar convolution has failed it is potentially for two reasons: 235 # 1) The kernel size for the boxcar is bigger than the actual image. 236 # 2) The grow parameter was specified with a value < 0. This would result 237 # in an illegal boxshape kernel. The dimensions of the kernel box *MUST* 238 # be integer and greater than zero. 239 # 240 # If the boxcar convolution has failed, try to give a meaningfull explanation 241 # as to why based upon the conditionals described above. 242 243 if (__boxsize <= 0): 244 errormsg1 = "############################################################\n" 245 errormsg1 += "# The boxcar convolution in minmed has failed. The 'grow' #\n" 246 errormsg1 += "# parameter must be greater than or equal to zero. You #\n" 247 errormsg1 += "# specified an input value for the 'grow' parameter of: #\n" 248 errormsg1 += " combine_grow: " + str(self.__combine_grow)+'\n' 249 errormsg1 += "############################################################\n" 250 raise ValueError,errormsg1 251 if (__boxsize > self.__imageList[0].shape[0]): 252 errormsg2 = "############################################################\n" 253 errormsg2 += "# The boxcar convolution in minmed has failed. The 'grow' #\n" 254 errormsg2 += "# parameter specified has resulted in a boxcar kernel that #\n" 255 errormsg2 += "# has dimensions larger than the actual image. You #\n" 256 errormsg2 += "# specified an input value for the 'grow' parameter of: #\n" 257 errormsg2 += " combine_grow: " +str(self.__combine_grow)+'\n' 258 errormsg2 += "############################################################\n" 259 print self.__imageList[0].shape 260 raise ValueError,errormsg2 261 262 # Attempt the boxcar convolution using the boxshape based upon the user input value of "grow" 263 NC.boxcar(__minimum_flag_file,__boxshape,output=__minimum_grow_file,mode='constant',cval=0) 264 265 del(__minimum_flag_file) 266 267 __temp1 = (__median_file_weighted - (__rms_file * self.__combine_nsigma1)) 268 __temp2 = (__median_file_weighted - (__rms_file * self.__combine_nsigma2)) 269 __median_rms2_file = np.where(np.equal(__minimum_grow_file,0),__temp1,__temp2) 270 del(__temp1) 271 del(__temp2) 272 del(__rms_file) 273 del(__minimum_grow_file) 274 275 # Finally decide whether to use the minimim or the median (in counts/s), 276 # based on whether the median is more than 3 sigma above the minimum. 277 # 278 self.combArrObj = __tmp 279 self.combArrObj = np.where(np.less(__minimum_file_weighted,__median_rms2_file), 280 __minimum_file, 281 __median_file) 282 283 else: 284 # Finally decide whether to use the minimim or the median (in counts/s), 285 # based on whether the median is more than 3 sigma above the minimum. 286 # 287 self.combArrObj = __tmp 288 self.combArrObj = np.where(np.less(__minimum_file_weighted,__median_rms_file), 289 __minimum_file, 290 __median_file) 291 292 # Set fill regions to a pixel value of 0. 293 self.combArrObj = np.where(__maskSum == self.__numberOfImages, 0, self.combArrObj)

294 295 # self.out_file1 = __median_rms2_file.copy() 296 # self.out_file2 = __minimum_file_weighted.copy() 297 298 299

300 - def __sumImages(self,numarrayObjectList):

301 """ Sum a list of numarray objects. """ 302 __sum = np.zeros(numarrayObjectList[0].shape,dtype=numarrayObjectList[0].dtype) 303 for image in numarrayObjectList: 304 __sum += image 305 return __sum

Source Code for Module multidrizzle.minmed