Package astLib :: Module astImages
[hide private]
[frames] | no frames]

Source Code for Module astLib.astImages

   1  """module for simple .fits image tasks (rotation, clipping out sections, making .pngs etc.) 
   2   
   3  (c) 2007-2013 Matt Hilton  
   4   
   5  U{http://astlib.sourceforge.net} 
   6   
   7  Some routines in this module will fail if, e.g., asked to clip a section from a .fits image at a 
   8  position not found within the image (as determined using the WCS). Where this occurs, the function 
   9  will return None. An error message will be printed to the console when this happens if 
  10  astImages.REPORT_ERRORS=True (the default). Testing if an astImages function returns None can be 
  11  used to handle errors in scripts.  
  12   
  13  """ 
  14   
  15  REPORT_ERRORS=True 
  16   
  17  import os 
  18  import sys 
  19  import math 
  20  from astLib import astWCS 
  21  import pyfits 
  22  try: 
  23      from scipy import ndimage 
  24      from scipy import interpolate 
  25  except: 
  26      print("WARNING: astImages: failed to import scipy.ndimage - some functions will not work.") 
  27  import numpy 
  28  try: 
  29      import matplotlib 
  30      from matplotlib import pylab 
  31      matplotlib.interactive(False) 
  32  except: 
  33      print("WARNING: astImages: failed to import matplotlib - some functions will not work.") 
  34   
  35  #--------------------------------------------------------------------------------------------------- 
36 -def clipImageSectionWCS(imageData, imageWCS, RADeg, decDeg, clipSizeDeg, returnWCS = True):
37 """Clips a square or rectangular section from an image array at the given celestial coordinates. 38 An updated WCS for the clipped section is optionally returned, as well as the x, y pixel 39 coordinates in the original image corresponding to the clipped section. 40 41 Note that the clip size is specified in degrees on the sky. For projections that have varying 42 real pixel scale across the map (e.g. CEA), use L{clipUsingRADecCoords} instead. 43 44 @type imageData: numpy array 45 @param imageData: image data array 46 @type imageWCS: astWCS.WCS 47 @param imageWCS: astWCS.WCS object 48 @type RADeg: float 49 @param RADeg: coordinate in decimal degrees 50 @type decDeg: float 51 @param decDeg: coordinate in decimal degrees 52 @type clipSizeDeg: float or list in format [widthDeg, heightDeg] 53 @param clipSizeDeg: if float, size of square clipped section in decimal degrees; if list, 54 size of clipped section in degrees in x, y axes of image respectively 55 @type returnWCS: bool 56 @param returnWCS: if True, return an updated WCS for the clipped section 57 @rtype: dictionary 58 @return: clipped image section (numpy array), updated astWCS WCS object for 59 clipped image section, and coordinates of clipped section in imageData in format 60 {'data', 'wcs', 'clippedSection'}. 61 62 """ 63 64 imHeight=imageData.shape[0] 65 imWidth=imageData.shape[1] 66 xImScale=imageWCS.getXPixelSizeDeg() 67 yImScale=imageWCS.getYPixelSizeDeg() 68 69 if type(clipSizeDeg) == float: 70 xHalfClipSizeDeg=clipSizeDeg/2.0 71 yHalfClipSizeDeg=xHalfClipSizeDeg 72 elif type(clipSizeDeg) == list or type(clipSizeDeg) == tuple: 73 xHalfClipSizeDeg=clipSizeDeg[0]/2.0 74 yHalfClipSizeDeg=clipSizeDeg[1]/2.0 75 else: 76 raise Exception("did not understand clipSizeDeg: should be float, or [widthDeg, heightDeg]") 77 78 xHalfSizePix=xHalfClipSizeDeg/xImScale 79 yHalfSizePix=yHalfClipSizeDeg/yImScale 80 81 cPixCoords=imageWCS.wcs2pix(RADeg, decDeg) 82 83 cTopLeft=[cPixCoords[0]+xHalfSizePix, cPixCoords[1]+yHalfSizePix] 84 cBottomRight=[cPixCoords[0]-xHalfSizePix, cPixCoords[1]-yHalfSizePix] 85 86 X=[int(round(cTopLeft[0])),int(round(cBottomRight[0]))] 87 Y=[int(round(cTopLeft[1])),int(round(cBottomRight[1]))] 88 89 X.sort() 90 Y.sort() 91 92 if X[0] < 0: 93 X[0]=0 94 if X[1] > imWidth: 95 X[1]=imWidth 96 if Y[0] < 0: 97 Y[0]=0 98 if Y[1] > imHeight: 99 Y[1]=imHeight 100 101 clippedData=imageData[Y[0]:Y[1],X[0]:X[1]] 102 103 # Update WCS 104 if returnWCS == True: 105 try: 106 oldCRPIX1=imageWCS.header['CRPIX1'] 107 oldCRPIX2=imageWCS.header['CRPIX2'] 108 clippedWCS=imageWCS.copy() 109 clippedWCS.header.update('NAXIS1', clippedData.shape[1]) 110 clippedWCS.header.update('NAXIS2', clippedData.shape[0]) 111 clippedWCS.header.update('CRPIX1', oldCRPIX1-X[0]) 112 clippedWCS.header.update('CRPIX2', oldCRPIX2-Y[0]) 113 clippedWCS.updateFromHeader() 114 115 except KeyError: 116 117 if REPORT_ERRORS == True: 118 119 print("WARNING: astImages.clipImageSectionWCS() : no CRPIX1, CRPIX2 keywords found - not updating clipped image WCS.") 120 121 clippedData=imageData[Y[0]:Y[1],X[0]:X[1]] 122 clippedWCS=imageWCS.copy() 123 else: 124 clippedWCS=None 125 126 return {'data': clippedData, 'wcs': clippedWCS, 'clippedSection': [X[0], X[1], Y[0], Y[1]]}
127 128 #---------------------------------------------------------------------------------------------------
129 -def clipImageSectionPix(imageData, XCoord, YCoord, clipSizePix):
130 """Clips a square or rectangular section from an image array at the given pixel coordinates. 131 132 @type imageData: numpy array 133 @param imageData: image data array 134 @type XCoord: float 135 @param XCoord: coordinate in pixels 136 @type YCoord: float 137 @param YCoord: coordinate in pixels 138 @type clipSizePix: float or list in format [widthPix, heightPix] 139 @param clipSizePix: if float, size of square clipped section in pixels; if list, 140 size of clipped section in pixels in x, y axes of output image respectively 141 @rtype: numpy array 142 @return: clipped image section 143 144 """ 145 146 imHeight=imageData.shape[0] 147 imWidth=imageData.shape[1] 148 149 if type(clipSizePix) == float or type(clipSizePix) == int: 150 xHalfClipSizePix=int(round(clipSizePix/2.0)) 151 yHalfClipSizePix=xHalfClipSizePix 152 elif type(clipSizePix) == list or type(clipSizePix) == tuple: 153 xHalfClipSizePix=int(round(clipSizePix[0]/2.0)) 154 yHalfClipSizePix=int(round(clipSizePix[1]/2.0)) 155 else: 156 raise Exception("did not understand clipSizePix: should be float, or [widthPix, heightPix]") 157 158 cTopLeft=[XCoord+xHalfClipSizePix, YCoord+yHalfClipSizePix] 159 cBottomRight=[XCoord-xHalfClipSizePix, YCoord-yHalfClipSizePix] 160 161 X=[int(round(cTopLeft[0])),int(round(cBottomRight[0]))] 162 Y=[int(round(cTopLeft[1])),int(round(cBottomRight[1]))] 163 164 X.sort() 165 Y.sort() 166 167 if X[0] < 0: 168 X[0]=0 169 if X[1] > imWidth: 170 X[1]=imWidth 171 if Y[0] < 0: 172 Y[0]=0 173 if Y[1] > imHeight: 174 Y[1]=imHeight 175 176 return imageData[Y[0]:Y[1],X[0]:X[1]]
177 178 #---------------------------------------------------------------------------------------------------
179 -def clipRotatedImageSectionWCS(imageData, imageWCS, RADeg, decDeg, clipSizeDeg, returnWCS = True):
180 """Clips a square or rectangular section from an image array at the given celestial coordinates. 181 The resulting clip is rotated and/or flipped such that North is at the top, and East appears at 182 the left. An updated WCS for the clipped section is also returned. Note that the alignment 183 of the rotated WCS is currently not perfect - however, it is probably good enough in most 184 cases for use with L{ImagePlot} for plotting purposes. 185 186 Note that the clip size is specified in degrees on the sky. For projections that have varying 187 real pixel scale across the map (e.g. CEA), use L{clipUsingRADecCoords} instead. 188 189 @type imageData: numpy array 190 @param imageData: image data array 191 @type imageWCS: astWCS.WCS 192 @param imageWCS: astWCS.WCS object 193 @type RADeg: float 194 @param RADeg: coordinate in decimal degrees 195 @type decDeg: float 196 @param decDeg: coordinate in decimal degrees 197 @type clipSizeDeg: float 198 @param clipSizeDeg: if float, size of square clipped section in decimal degrees; if list, 199 size of clipped section in degrees in RA, dec. axes of output rotated image respectively 200 @type returnWCS: bool 201 @param returnWCS: if True, return an updated WCS for the clipped section 202 @rtype: dictionary 203 @return: clipped image section (numpy array), updated astWCS WCS object for 204 clipped image section, in format {'data', 'wcs'}. 205 206 @note: Returns 'None' if the requested position is not found within the image. If the image 207 WCS does not have keywords of the form CD1_1 etc., the output WCS will not be rotated. 208 209 """ 210 211 halfImageSize=imageWCS.getHalfSizeDeg() 212 imageCentre=imageWCS.getCentreWCSCoords() 213 imScale=imageWCS.getPixelSizeDeg() 214 215 if type(clipSizeDeg) == float: 216 xHalfClipSizeDeg=clipSizeDeg/2.0 217 yHalfClipSizeDeg=xHalfClipSizeDeg 218 elif type(clipSizeDeg) == list or type(clipSizeDeg) == tuple: 219 xHalfClipSizeDeg=clipSizeDeg[0]/2.0 220 yHalfClipSizeDeg=clipSizeDeg[1]/2.0 221 else: 222 raise Exception("did not understand clipSizeDeg: should be float, or [widthDeg, heightDeg]") 223 224 diagonalHalfSizeDeg=math.sqrt((xHalfClipSizeDeg*xHalfClipSizeDeg) \ 225 +(yHalfClipSizeDeg*yHalfClipSizeDeg)) 226 227 diagonalHalfSizePix=diagonalHalfSizeDeg/imScale 228 229 if RADeg>imageCentre[0]-halfImageSize[0] and RADeg<imageCentre[0]+halfImageSize[0] \ 230 and decDeg>imageCentre[1]-halfImageSize[1] and decDeg<imageCentre[1]+halfImageSize[1]: 231 232 imageDiagonalClip=clipImageSectionWCS(imageData, imageWCS, RADeg, 233 decDeg, diagonalHalfSizeDeg*2.0) 234 diagonalClip=imageDiagonalClip['data'] 235 diagonalWCS=imageDiagonalClip['wcs'] 236 237 rotDeg=diagonalWCS.getRotationDeg() 238 imageRotated=ndimage.rotate(diagonalClip, rotDeg) 239 if diagonalWCS.isFlipped() == 1: 240 imageRotated=pylab.fliplr(imageRotated) 241 242 # Handle WCS rotation 243 rotatedWCS=diagonalWCS.copy() 244 rotRadians=math.radians(rotDeg) 245 246 if returnWCS == True: 247 try: 248 249 CD11=rotatedWCS.header['CD1_1'] 250 CD21=rotatedWCS.header['CD2_1'] 251 CD12=rotatedWCS.header['CD1_2'] 252 CD22=rotatedWCS.header['CD2_2'] 253 if rotatedWCS.isFlipped() == 1: 254 CD11=CD11*-1 255 CD12=CD12*-1 256 CDMatrix=numpy.array([[CD11, CD12], [CD21, CD22]], dtype=numpy.float64) 257 258 rotRadians=rotRadians 259 rot11=math.cos(rotRadians) 260 rot12=math.sin(rotRadians) 261 rot21=-math.sin(rotRadians) 262 rot22=math.cos(rotRadians) 263 rotMatrix=numpy.array([[rot11, rot12], [rot21, rot22]], dtype=numpy.float64) 264 newCDMatrix=numpy.dot(rotMatrix, CDMatrix) 265 266 P1=diagonalWCS.header['CRPIX1'] 267 P2=diagonalWCS.header['CRPIX2'] 268 V1=diagonalWCS.header['CRVAL1'] 269 V2=diagonalWCS.header['CRVAL2'] 270 271 PMatrix=numpy.zeros((2,), dtype = numpy.float64) 272 PMatrix[0]=P1 273 PMatrix[1]=P2 274 275 # BELOW IS HOW TO WORK OUT THE NEW REF PIXEL 276 CMatrix=numpy.array([imageRotated.shape[1]/2.0, imageRotated.shape[0]/2.0]) 277 centreCoords=diagonalWCS.getCentreWCSCoords() 278 alphaRad=math.radians(centreCoords[0]) 279 deltaRad=math.radians(centreCoords[1]) 280 thetaRad=math.asin(math.sin(deltaRad)*math.sin(math.radians(V2)) + \ 281 math.cos(deltaRad)*math.cos(math.radians(V2))*math.cos(alphaRad-math.radians(V1))) 282 phiRad=math.atan2(-math.cos(deltaRad)*math.sin(alphaRad-math.radians(V1)), \ 283 math.sin(deltaRad)*math.cos(math.radians(V2)) - \ 284 math.cos(deltaRad)*math.sin(math.radians(V2))*math.cos(alphaRad-math.radians(V1))) + \ 285 math.pi 286 RTheta=(180.0/math.pi)*(1.0/math.tan(thetaRad)) 287 288 xy=numpy.zeros((2,), dtype=numpy.float64) 289 xy[0]=RTheta*math.sin(phiRad) 290 xy[1]=-RTheta*math.cos(phiRad) 291 newPMatrix=CMatrix - numpy.dot(numpy.linalg.inv(newCDMatrix), xy) 292 293 # But there's a small offset to CRPIX due to the rotatedImage being rounded to an integer 294 # number of pixels (not sure this helps much) 295 #d=numpy.dot(rotMatrix, [diagonalClip.shape[1], diagonalClip.shape[0]]) 296 #offset=abs(d)-numpy.array(imageRotated.shape) 297 298 rotatedWCS.header.update('NAXIS1', imageRotated.shape[1]) 299 rotatedWCS.header.update('NAXIS2', imageRotated.shape[0]) 300 rotatedWCS.header.update('CRPIX1', newPMatrix[0]) 301 rotatedWCS.header.update('CRPIX2', newPMatrix[1]) 302 rotatedWCS.header.update('CRVAL1', V1) 303 rotatedWCS.header.update('CRVAL2', V2) 304 rotatedWCS.header.update('CD1_1', newCDMatrix[0][0]) 305 rotatedWCS.header.update('CD2_1', newCDMatrix[1][0]) 306 rotatedWCS.header.update('CD1_2', newCDMatrix[0][1]) 307 rotatedWCS.header.update('CD2_2', newCDMatrix[1][1]) 308 rotatedWCS.updateFromHeader() 309 310 except KeyError: 311 312 if REPORT_ERRORS == True: 313 print("WARNING: astImages.clipRotatedImageSectionWCS() : no CDi_j keywords found - not rotating WCS.") 314 315 imageRotated=diagonalClip 316 rotatedWCS=diagonalWCS 317 318 imageRotatedClip=clipImageSectionWCS(imageRotated, rotatedWCS, RADeg, decDeg, clipSizeDeg) 319 320 if returnWCS == True: 321 return {'data': imageRotatedClip['data'], 'wcs': imageRotatedClip['wcs']} 322 else: 323 return {'data': imageRotatedClip['data'], 'wcs': None} 324 325 else: 326 327 if REPORT_ERRORS==True: 328 print("""ERROR: astImages.clipRotatedImageSectionWCS() : 329 RADeg, decDeg are not within imageData.""") 330 331 return None
332 333 #---------------------------------------------------------------------------------------------------
334 -def clipUsingRADecCoords(imageData, imageWCS, RAMin, RAMax, decMin, decMax, returnWCS = True):
335 """Clips a section from an image array at the pixel coordinates corresponding to the given 336 celestial coordinates. 337 338 @type imageData: numpy array 339 @param imageData: image data array 340 @type imageWCS: astWCS.WCS 341 @param imageWCS: astWCS.WCS object 342 @type RAMin: float 343 @param RAMin: minimum RA coordinate in decimal degrees 344 @type RAMax: float 345 @param RAMax: maximum RA coordinate in decimal degrees 346 @type decMin: float 347 @param decMin: minimum dec coordinate in decimal degrees 348 @type decMax: float 349 @param decMax: maximum dec coordinate in decimal degrees 350 @type returnWCS: bool 351 @param returnWCS: if True, return an updated WCS for the clipped section 352 @rtype: dictionary 353 @return: clipped image section (numpy array), updated astWCS WCS object for 354 clipped image section, and corresponding pixel coordinates in imageData in format 355 {'data', 'wcs', 'clippedSection'}. 356 357 @note: Returns 'None' if the requested position is not found within the image. 358 359 """ 360 361 imHeight=imageData.shape[0] 362 imWidth=imageData.shape[1] 363 364 xMin, yMin=imageWCS.wcs2pix(RAMin, decMin) 365 xMax, yMax=imageWCS.wcs2pix(RAMax, decMax) 366 xMin=int(round(xMin)) 367 xMax=int(round(xMax)) 368 yMin=int(round(yMin)) 369 yMax=int(round(yMax)) 370 X=[xMin, xMax] 371 X.sort() 372 Y=[yMin, yMax] 373 Y.sort() 374 375 if X[0] < 0: 376 X[0]=0 377 if X[1] > imWidth: 378 X[1]=imWidth 379 if Y[0] < 0: 380 Y[0]=0 381 if Y[1] > imHeight: 382 Y[1]=imHeight 383 384 clippedData=imageData[Y[0]:Y[1],X[0]:X[1]] 385 386 # Update WCS 387 if returnWCS == True: 388 try: 389 oldCRPIX1=imageWCS.header['CRPIX1'] 390 oldCRPIX2=imageWCS.header['CRPIX2'] 391 clippedWCS=imageWCS.copy() 392 clippedWCS.header.update('NAXIS1', clippedData.shape[1]) 393 clippedWCS.header.update('NAXIS2', clippedData.shape[0]) 394 clippedWCS.header.update('CRPIX1', oldCRPIX1-X[0]) 395 clippedWCS.header.update('CRPIX2', oldCRPIX2-Y[0]) 396 clippedWCS.updateFromHeader() 397 398 except KeyError: 399 400 if REPORT_ERRORS == True: 401 402 print("WARNING: astImages.clipUsingRADecCoords() : no CRPIX1, CRPIX2 keywords found - not updating clipped image WCS.") 403 404 clippedData=imageData[Y[0]:Y[1],X[0]:X[1]] 405 clippedWCS=imageWCS.copy() 406 else: 407 clippedWCS=None 408 409 return {'data': clippedData, 'wcs': clippedWCS, 'clippedSection': [X[0], X[1], Y[0], Y[1]]}
410 411 #---------------------------------------------------------------------------------------------------
412 -def scaleImage(imageData, imageWCS, scaleFactor):
413 """Scales image array and WCS by the given scale factor. 414 415 @type imageData: numpy array 416 @param imageData: image data array 417 @type imageWCS: astWCS.WCS 418 @param imageWCS: astWCS.WCS object 419 @type scaleFactor: float or list or tuple 420 @param scaleFactor: factor to resize image by - if tuple or list, in format 421 [x scale factor, y scale factor] 422 @rtype: dictionary 423 @return: image data (numpy array), updated astWCS WCS object for image, in format {'data', 'wcs'}. 424 425 """ 426 427 if type(scaleFactor) == int or type(scaleFactor) == float: 428 scaleFactor=[float(scaleFactor), float(scaleFactor)] 429 scaledData=ndimage.zoom(imageData, scaleFactor) 430 431 # Take care of offset due to rounding in scaling image to integer pixel dimensions 432 properDimensions=numpy.array(imageData.shape)*scaleFactor 433 offset=properDimensions-numpy.array(scaledData.shape) 434 435 # Rescale WCS 436 try: 437 oldCRPIX1=imageWCS.header['CRPIX1'] 438 oldCRPIX2=imageWCS.header['CRPIX2'] 439 CD11=imageWCS.header['CD1_1'] 440 CD21=imageWCS.header['CD2_1'] 441 CD12=imageWCS.header['CD1_2'] 442 CD22=imageWCS.header['CD2_2'] 443 except KeyError: 444 # Try the older FITS header format 445 try: 446 oldCRPIX1=imageWCS.header['CRPIX1'] 447 oldCRPIX2=imageWCS.header['CRPIX2'] 448 CD11=imageWCS.header['CDELT1'] 449 CD21=0 450 CD12=0 451 CD22=imageWCS.header['CDELT2'] 452 except KeyError: 453 if REPORT_ERRORS == True: 454 print("WARNING: astImages.rescaleImage() : no CDij or CDELT keywords found - not updating WCS.") 455 scaledWCS=imageWCS.copy() 456 return {'data': scaledData, 'wcs': scaledWCS} 457 458 CDMatrix=numpy.array([[CD11, CD12], [CD21, CD22]], dtype=numpy.float64) 459 scaleFactorMatrix=numpy.array([[1.0/scaleFactor[0], 0], [0, 1.0/scaleFactor[1]]]) 460 scaledCDMatrix=numpy.dot(scaleFactorMatrix, CDMatrix) 461 462 scaledWCS=imageWCS.copy() 463 scaledWCS.header.update('NAXIS1', scaledData.shape[1]) 464 scaledWCS.header.update('NAXIS2', scaledData.shape[0]) 465 scaledWCS.header.update('CRPIX1', oldCRPIX1*scaleFactor[0]+offset[1]) 466 scaledWCS.header.update('CRPIX2', oldCRPIX2*scaleFactor[1]+offset[0]) 467 scaledWCS.header.update('CD1_1', scaledCDMatrix[0][0]) 468 scaledWCS.header.update('CD2_1', scaledCDMatrix[1][0]) 469 scaledWCS.header.update('CD1_2', scaledCDMatrix[0][1]) 470 scaledWCS.header.update('CD2_2', scaledCDMatrix[1][1]) 471 scaledWCS.updateFromHeader() 472 473 return {'data': scaledData, 'wcs': scaledWCS}
474 475 #---------------------------------------------------------------------------------------------------
476 -def intensityCutImage(imageData, cutLevels):
477 """Creates a matplotlib.pylab plot of an image array with the specified cuts in intensity 478 applied. This routine is used by L{saveBitmap} and L{saveContourOverlayBitmap}, which both 479 produce output as .png, .jpg, etc. images. 480 481 @type imageData: numpy array 482 @param imageData: image data array 483 @type cutLevels: list 484 @param cutLevels: sets the image scaling - available options: 485 - pixel values: cutLevels=[low value, high value]. 486 - histogram equalisation: cutLevels=["histEq", number of bins ( e.g. 1024)] 487 - relative: cutLevels=["relative", cut per cent level (e.g. 99.5)] 488 - smart: cutLevels=["smart", cut per cent level (e.g. 99.5)] 489 ["smart", 99.5] seems to provide good scaling over a range of different images. 490 @rtype: dictionary 491 @return: image section (numpy.array), matplotlib image normalisation (matplotlib.colors.Normalize), in the format {'image', 'norm'}. 492 493 @note: If cutLevels[0] == "histEq", then only {'image'} is returned. 494 495 """ 496 497 oImWidth=imageData.shape[1] 498 oImHeight=imageData.shape[0] 499 500 # Optional histogram equalisation 501 if cutLevels[0]=="histEq": 502 503 imageData=histEq(imageData, cutLevels[1]) 504 anorm=pylab.normalize(imageData.min(), imageData.max()) 505 506 elif cutLevels[0]=="relative": 507 508 # this turns image data into 1D array then sorts 509 sorted=numpy.sort(numpy.ravel(imageData)) 510 maxValue=sorted.max() 511 minValue=sorted.min() 512 513 # want to discard the top and bottom specified 514 topCutIndex=len(sorted-1) \ 515 -int(math.floor(float((100.0-cutLevels[1])/100.0)*len(sorted-1))) 516 bottomCutIndex=int(math.ceil(float((100.0-cutLevels[1])/100.0)*len(sorted-1))) 517 topCut=sorted[topCutIndex] 518 bottomCut=sorted[bottomCutIndex] 519 anorm=pylab.normalize(bottomCut, topCut) 520 521 elif cutLevels[0]=="smart": 522 523 # this turns image data into 1Darray then sorts 524 sorted=numpy.sort(numpy.ravel(imageData)) 525 maxValue=sorted.max() 526 minValue=sorted.min() 527 numBins=10000 # 0.01 per cent accuracy 528 binWidth=(maxValue-minValue)/float(numBins) 529 histogram=ndimage.histogram(sorted, minValue, maxValue, numBins) 530 531 # Find the bin with the most pixels in it, set that as our minimum 532 # Then search through the bins until we get to a bin with more/or the same number of 533 # pixels in it than the previous one. 534 # We take that to be the maximum. 535 # This means that we avoid the traps of big, bright, saturated stars that cause 536 # problems for relative scaling 537 backgroundValue=histogram.max() 538 foundBackgroundBin=False 539 foundTopBin=False 540 lastBin=-10000 541 for i in range(len(histogram)): 542 543 if histogram[i]>=lastBin and foundBackgroundBin==True: 544 545 # Added a fudge here to stop us picking for top bin a bin within 546 # 10 percent of the background pixel value 547 if (minValue+(binWidth*i))>bottomBinValue*1.1: 548 topBinValue=minValue+(binWidth*i) 549 foundTopBin=True 550 break 551 552 if histogram[i]==backgroundValue and foundBackgroundBin==False: 553 bottomBinValue=minValue+(binWidth*i) 554 foundBackgroundBin=True 555 556 lastBin=histogram[i] 557 558 if foundTopBin==False: 559 topBinValue=maxValue 560 561 #Now we apply relative scaling to this 562 smartClipped=numpy.clip(sorted, bottomBinValue, topBinValue) 563 topCutIndex=len(smartClipped-1) \ 564 -int(math.floor(float((100.0-cutLevels[1])/100.0)*len(smartClipped-1))) 565 bottomCutIndex=int(math.ceil(float((100.0-cutLevels[1])/100.0)*len(smartClipped-1))) 566 topCut=smartClipped[topCutIndex] 567 bottomCut=smartClipped[bottomCutIndex] 568 anorm=pylab.normalize(bottomCut, topCut) 569 else: 570 571 # Normalise using given cut levels 572 anorm=pylab.normalize(cutLevels[0], cutLevels[1]) 573 574 if cutLevels[0]=="histEq": 575 return {'image': imageData.copy()} 576 else: 577 return {'image': imageData.copy(), 'norm': anorm}
578 579 #---------------------------------------------------------------------------------------------------
580 -def resampleToTanProjection(imageData, imageWCS, outputPixDimensions=[600, 600]):
581 """Resamples an image and WCS to a tangent plane projection. Purely for plotting purposes 582 (e.g., ensuring RA, dec. coordinate axes perpendicular). 583 584 @type imageData: numpy array 585 @param imageData: image data array 586 @type imageWCS: astWCS.WCS 587 @param imageWCS: astWCS.WCS object 588 @type outputPixDimensions: list 589 @param outputPixDimensions: [width, height] of output image in pixels 590 @rtype: dictionary 591 @return: image data (numpy array), updated astWCS WCS object for image, in format {'data', 'wcs'}. 592 593 """ 594 595 RADeg, decDeg=imageWCS.getCentreWCSCoords() 596 xPixelScale=imageWCS.getXPixelSizeDeg() 597 yPixelScale=imageWCS.getYPixelSizeDeg() 598 xSizeDeg, ySizeDeg=imageWCS.getFullSizeSkyDeg() 599 xSizePix=int(round(outputPixDimensions[0])) 600 ySizePix=int(round(outputPixDimensions[1])) 601 xRefPix=xSizePix/2.0 602 yRefPix=ySizePix/2.0 603 xOutPixScale=xSizeDeg/xSizePix 604 yOutPixScale=ySizeDeg/ySizePix 605 cardList=pyfits.CardList() 606 cardList.append(pyfits.Card('NAXIS', 2)) 607 cardList.append(pyfits.Card('NAXIS1', xSizePix)) 608 cardList.append(pyfits.Card('NAXIS2', ySizePix)) 609 cardList.append(pyfits.Card('CTYPE1', 'RA---TAN')) 610 cardList.append(pyfits.Card('CTYPE2', 'DEC--TAN')) 611 cardList.append(pyfits.Card('CRVAL1', RADeg)) 612 cardList.append(pyfits.Card('CRVAL2', decDeg)) 613 cardList.append(pyfits.Card('CRPIX1', xRefPix+1)) 614 cardList.append(pyfits.Card('CRPIX2', yRefPix+1)) 615 cardList.append(pyfits.Card('CDELT1', -xOutPixScale)) 616 cardList.append(pyfits.Card('CDELT2', xOutPixScale)) # Makes more sense to use same pix scale 617 cardList.append(pyfits.Card('CUNIT1', 'DEG')) 618 cardList.append(pyfits.Card('CUNIT2', 'DEG')) 619 newHead=pyfits.Header(cards=cardList) 620 newWCS=astWCS.WCS(newHead, mode='pyfits') 621 newImage=numpy.zeros([ySizePix, xSizePix]) 622 623 tanImage=resampleToWCS(newImage, newWCS, imageData, imageWCS, highAccuracy=True, 624 onlyOverlapping=False) 625 626 return tanImage
627 628 #---------------------------------------------------------------------------------------------------
629 -def resampleToWCS(im1Data, im1WCS, im2Data, im2WCS, highAccuracy = False, onlyOverlapping = True):
630 """Resamples data corresponding to second image (with data im2Data, WCS im2WCS) onto the WCS 631 of the first image (im1Data, im1WCS). The output, resampled image is of the pixel same 632 dimensions of the first image. This routine is for assisting in plotting - performing 633 photometry on the output is not recommended. 634 635 Set highAccuracy == True to sample every corresponding pixel in each image; otherwise only 636 every nth pixel (where n is the ratio of the image scales) will be sampled, with values 637 in between being set using a linear interpolation (much faster). 638 639 Set onlyOverlapping == True to speed up resampling by only resampling the overlapping 640 area defined by both image WCSs. 641 642 @type im1Data: numpy array 643 @param im1Data: image data array for first image 644 @type im1WCS: astWCS.WCS 645 @param im1WCS: astWCS.WCS object corresponding to im1Data 646 @type im2Data: numpy array 647 @param im2Data: image data array for second image (to be resampled to match first image) 648 @type im2WCS: astWCS.WCS 649 @param im2WCS: astWCS.WCS object corresponding to im2Data 650 @type highAccuracy: bool 651 @param highAccuracy: if True, sample every corresponding pixel in each image; otherwise, sample 652 every nth pixel, where n = the ratio of the image scales. 653 @type onlyOverlapping: bool 654 @param onlyOverlapping: if True, only consider the overlapping area defined by both image WCSs 655 (speeds things up) 656 @rtype: dictionary 657 @return: numpy image data array and associated WCS in format {'data', 'wcs'} 658 659 """ 660 661 resampledData=numpy.zeros(im1Data.shape) 662 663 # Find overlap - speed things up 664 # But have a border so as not to require the overlap to be perfect 665 # There's also no point in oversampling image 1 if it's much higher res than image 2 666 xPixRatio=(im2WCS.getXPixelSizeDeg()/im1WCS.getXPixelSizeDeg())/2.0 667 yPixRatio=(im2WCS.getYPixelSizeDeg()/im1WCS.getYPixelSizeDeg())/2.0 668 xBorder=xPixRatio*10.0 669 yBorder=yPixRatio*10.0 670 if highAccuracy == False: 671 if xPixRatio > 1: 672 xPixStep=int(math.ceil(xPixRatio)) 673 else: 674 xPixStep=1 675 if yPixRatio > 1: 676 yPixStep=int(math.ceil(yPixRatio)) 677 else: 678 yPixStep=1 679 else: 680 xPixStep=1 681 yPixStep=1 682 683 if onlyOverlapping == True: 684 overlap=astWCS.findWCSOverlap(im1WCS, im2WCS) 685 xOverlap=[overlap['wcs1Pix'][0], overlap['wcs1Pix'][1]] 686 yOverlap=[overlap['wcs1Pix'][2], overlap['wcs1Pix'][3]] 687 xOverlap.sort() 688 yOverlap.sort() 689 xMin=int(math.floor(xOverlap[0]-xBorder)) 690 xMax=int(math.ceil(xOverlap[1]+xBorder)) 691 yMin=int(math.floor(yOverlap[0]-yBorder)) 692 yMax=int(math.ceil(yOverlap[1]+yBorder)) 693 xRemainder=(xMax-xMin) % xPixStep 694 yRemainder=(yMax-yMin) % yPixStep 695 if xRemainder != 0: 696 xMax=xMax+xRemainder 697 if yRemainder != 0: 698 yMax=yMax+yRemainder 699 # Check that we're still within the image boundaries, to be on the safe side 700 if xMin < 0: 701 xMin=0 702 if xMax > im1Data.shape[1]: 703 xMax=im1Data.shape[1] 704 if yMin < 0: 705 yMin=0 706 if yMax > im1Data.shape[0]: 707 yMax=im1Data.shape[0] 708 else: 709 xMin=0 710 xMax=im1Data.shape[1] 711 yMin=0 712 yMax=im1Data.shape[0] 713 714 for x in range(xMin, xMax, xPixStep): 715 for y in range(yMin, yMax, yPixStep): 716 RA, dec=im1WCS.pix2wcs(x, y) 717 x2, y2=im2WCS.wcs2pix(RA, dec) 718 x2=int(round(x2)) 719 y2=int(round(y2)) 720 if x2 >= 0 and x2 < im2Data.shape[1] and y2 >= 0 and y2 < im2Data.shape[0]: 721 resampledData[y][x]=im2Data[y2][x2] 722 723 # linear interpolation 724 if highAccuracy == False: 725 for row in range(resampledData.shape[0]): 726 vals=resampledData[row, numpy.arange(xMin, xMax, xPixStep)] 727 index2data=interpolate.interp1d(numpy.arange(0, vals.shape[0], 1), vals) 728 interpedVals=index2data(numpy.arange(0, vals.shape[0]-1, 1.0/xPixStep)) 729 resampledData[row, xMin:xMin+interpedVals.shape[0]]=interpedVals 730 for col in range(resampledData.shape[1]): 731 vals=resampledData[numpy.arange(yMin, yMax, yPixStep), col] 732 index2data=interpolate.interp1d(numpy.arange(0, vals.shape[0], 1), vals) 733 interpedVals=index2data(numpy.arange(0, vals.shape[0]-1, 1.0/yPixStep)) 734 resampledData[yMin:yMin+interpedVals.shape[0], col]=interpedVals 735 736 # Note: should really just copy im1WCS keywords into im2WCS and return that 737 # Only a problem if we're using this for anything other than plotting 738 return {'data': resampledData, 'wcs': im1WCS.copy()}
739 740 #---------------------------------------------------------------------------------------------------
741 -def generateContourOverlay(backgroundImageData, backgroundImageWCS, contourImageData, contourImageWCS, \ 742 contourLevels, contourSmoothFactor = 0, highAccuracy = False):
743 """Rescales an image array to be used as a contour overlay to have the same dimensions as the 744 background image, and generates a set of contour levels. The image array from which the contours 745 are to be generated will be resampled to the same dimensions as the background image data, and 746 can be optionally smoothed using a Gaussian filter. The sigma of the Gaussian filter 747 (contourSmoothFactor) is specified in arcsec. 748 749 @type backgroundImageData: numpy array 750 @param backgroundImageData: background image data array 751 @type backgroundImageWCS: astWCS.WCS 752 @param backgroundImageWCS: astWCS.WCS object of the background image data array 753 @type contourImageData: numpy array 754 @param contourImageData: image data array from which contours are to be generated 755 @type contourImageWCS: astWCS.WCS 756 @param contourImageWCS: astWCS.WCS object corresponding to contourImageData 757 @type contourLevels: list 758 @param contourLevels: sets the contour levels - available options: 759 - values: contourLevels=[list of values specifying each level] 760 - linear spacing: contourLevels=['linear', min level value, max level value, number 761 of levels] - can use "min", "max" to automatically set min, max levels from image data 762 - log spacing: contourLevels=['log', min level value, max level value, number of 763 levels] - can use "min", "max" to automatically set min, max levels from image data 764 @type contourSmoothFactor: float 765 @param contourSmoothFactor: standard deviation (in arcsec) of Gaussian filter for 766 pre-smoothing of contour image data (set to 0 for no smoothing) 767 @type highAccuracy: bool 768 @param highAccuracy: if True, sample every corresponding pixel in each image; otherwise, sample 769 every nth pixel, where n = the ratio of the image scales. 770 771 """ 772 773 # For compromise between speed and accuracy, scale a copy of the background 774 # image down to a scale that is one pixel = 1/5 of a pixel in the contour image 775 # But only do this if it has CDij keywords as we know how to scale those 776 if ("CD1_1" in backgroundImageWCS.header) == True: 777 xScaleFactor=backgroundImageWCS.getXPixelSizeDeg()/(contourImageWCS.getXPixelSizeDeg()/5.0) 778 yScaleFactor=backgroundImageWCS.getYPixelSizeDeg()/(contourImageWCS.getYPixelSizeDeg()/5.0) 779 scaledBackground=scaleImage(backgroundImageData, backgroundImageWCS, (xScaleFactor, yScaleFactor)) 780 scaled=resampleToWCS(scaledBackground['data'], scaledBackground['wcs'], 781 contourImageData, contourImageWCS, highAccuracy = highAccuracy) 782 scaledContourData=scaled['data'] 783 scaledContourWCS=scaled['wcs'] 784 scaledBackground=True 785 else: 786 scaled=resampleToWCS(backgroundImageData, backgroundImageWCS, 787 contourImageData, contourImageWCS, highAccuracy = highAccuracy) 788 scaledContourData=scaled['data'] 789 scaledContourWCS=scaled['wcs'] 790 scaledBackground=False 791 792 if contourSmoothFactor > 0: 793 sigmaPix=(contourSmoothFactor/3600.0)/scaledContourWCS.getPixelSizeDeg() 794 scaledContourData=ndimage.gaussian_filter(scaledContourData, sigmaPix) 795 796 # Various ways of setting the contour levels 797 # If just a list is passed in, use those instead 798 if contourLevels[0] == "linear": 799 if contourLevels[1] == "min": 800 xMin=contourImageData.flatten().min() 801 else: 802 xMin=float(contourLevels[1]) 803 if contourLevels[2] == "max": 804 xMax=contourImageData.flatten().max() 805 else: 806 xMax=float(contourLevels[2]) 807 nLevels=contourLevels[3] 808 xStep=(xMax-xMin)/(nLevels-1) 809 cLevels=[] 810 for j in range(nLevels+1): 811 level=xMin+j*xStep 812 cLevels.append(level) 813 814 elif contourLevels[0] == "log": 815 if contourLevels[1] == "min": 816 xMin=contourImageData.flatten().min() 817 else: 818 xMin=float(contourLevels[1]) 819 if contourLevels[2] == "max": 820 xMax=contourImageData.flatten().max() 821 else: 822 xMax=float(contourLevels[2]) 823 if xMin <= 0.0: 824 raise Exception("minimum contour level set to <= 0 and log scaling chosen.") 825 xLogMin=math.log10(xMin) 826 xLogMax=math.log10(xMax) 827 nLevels=contourLevels[3] 828 xLogStep=(xLogMax-xLogMin)/(nLevels-1) 829 cLevels=[] 830 prevLevel=0 831 for j in range(nLevels+1): 832 level=math.pow(10, xLogMin+j*xLogStep) 833 cLevels.append(level) 834 835 else: 836 cLevels=contourLevels 837 838 # Now blow the contour image data back up to the size of the original image 839 if scaledBackground == True: 840 scaledBack=scaleImage(scaledContourData, scaledContourWCS, (1.0/xScaleFactor, 1.0/yScaleFactor))['data'] 841 else: 842 scaledBack=scaledContourData 843 844 return {'scaledImage': scaledBack, 'contourLevels': cLevels}
845 846 #---------------------------------------------------------------------------------------------------
847 -def saveBitmap(outputFileName, imageData, cutLevels, size, colorMapName):
848 """Makes a bitmap image from an image array; the image format is specified by the 849 filename extension. (e.g. ".jpg" =JPEG, ".png"=PNG). 850 851 @type outputFileName: string 852 @param outputFileName: filename of output bitmap image 853 @type imageData: numpy array 854 @param imageData: image data array 855 @type cutLevels: list 856 @param cutLevels: sets the image scaling - available options: 857 - pixel values: cutLevels=[low value, high value]. 858 - histogram equalisation: cutLevels=["histEq", number of bins ( e.g. 1024)] 859 - relative: cutLevels=["relative", cut per cent level (e.g. 99.5)] 860 - smart: cutLevels=["smart", cut per cent level (e.g. 99.5)] 861 ["smart", 99.5] seems to provide good scaling over a range of different images. 862 @type size: int 863 @param size: size of output image in pixels 864 @type colorMapName: string 865 @param colorMapName: name of a standard matplotlib colormap, e.g. "hot", "cool", "gray" 866 etc. (do "help(pylab.colormaps)" in the Python interpreter to see available options) 867 868 """ 869 870 cut=intensityCutImage(imageData, cutLevels) 871 872 # Make plot 873 aspectR=float(cut['image'].shape[0])/float(cut['image'].shape[1]) 874 pylab.figure(figsize=(10,10*aspectR)) 875 pylab.axes([0,0,1,1]) 876 877 try: 878 colorMap=pylab.cm.get_cmap(colorMapName) 879 except AssertionError: 880 raise Exception(colorMapName+" is not a defined matplotlib colormap.") 881 882 if cutLevels[0]=="histEq": 883 pylab.imshow(cut['image'], interpolation="bilinear", origin='lower', cmap=colorMap) 884 885 else: 886 pylab.imshow(cut['image'], interpolation="bilinear", norm=cut['norm'], origin='lower', 887 cmap=colorMap) 888 889 pylab.axis("off") 890 891 pylab.savefig("out_astImages.png") 892 pylab.close("all") 893 894 try: 895 from PIL import Image 896 except: 897 raise Exception("astImages.saveBitmap requires the Python Imaging Library to be installed.") 898 im=Image.open("out_astImages.png") 899 im.thumbnail((int(size),int(size))) 900 im.save(outputFileName) 901 902 os.remove("out_astImages.png")
903 904 #---------------------------------------------------------------------------------------------------
905 -def saveContourOverlayBitmap(outputFileName, backgroundImageData, backgroundImageWCS, cutLevels, \ 906 size, colorMapName, contourImageData, contourImageWCS, \ 907 contourSmoothFactor, contourLevels, contourColor, contourWidth):
908 """Makes a bitmap image from an image array, with a set of contours generated from a 909 second image array overlaid. The image format is specified by the file extension 910 (e.g. ".jpg"=JPEG, ".png"=PNG). The image array from which the contours are to be generated 911 can optionally be pre-smoothed using a Gaussian filter. 912 913 @type outputFileName: string 914 @param outputFileName: filename of output bitmap image 915 @type backgroundImageData: numpy array 916 @param backgroundImageData: background image data array 917 @type backgroundImageWCS: astWCS.WCS 918 @param backgroundImageWCS: astWCS.WCS object of the background image data array 919 @type cutLevels: list 920 @param cutLevels: sets the image scaling - available options: 921 - pixel values: cutLevels=[low value, high value]. 922 - histogram equalisation: cutLevels=["histEq", number of bins ( e.g. 1024)] 923 - relative: cutLevels=["relative", cut per cent level (e.g. 99.5)] 924 - smart: cutLevels=["smart", cut per cent level (e.g. 99.5)] 925 ["smart", 99.5] seems to provide good scaling over a range of different images. 926 @type size: int 927 @param size: size of output image in pixels 928 @type colorMapName: string 929 @param colorMapName: name of a standard matplotlib colormap, e.g. "hot", "cool", "gray" 930 etc. (do "help(pylab.colormaps)" in the Python interpreter to see available options) 931 @type contourImageData: numpy array 932 @param contourImageData: image data array from which contours are to be generated 933 @type contourImageWCS: astWCS.WCS 934 @param contourImageWCS: astWCS.WCS object corresponding to contourImageData 935 @type contourSmoothFactor: float 936 @param contourSmoothFactor: standard deviation (in pixels) of Gaussian filter for 937 pre-smoothing of contour image data (set to 0 for no smoothing) 938 @type contourLevels: list 939 @param contourLevels: sets the contour levels - available options: 940 - values: contourLevels=[list of values specifying each level] 941 - linear spacing: contourLevels=['linear', min level value, max level value, number 942 of levels] - can use "min", "max" to automatically set min, max levels from image data 943 - log spacing: contourLevels=['log', min level value, max level value, number of 944 levels] - can use "min", "max" to automatically set min, max levels from image data 945 @type contourColor: string 946 @param contourColor: color of the overlaid contours, specified by the name of a standard 947 matplotlib color, e.g., "black", "white", "cyan" 948 etc. (do "help(pylab.colors)" in the Python interpreter to see available options) 949 @type contourWidth: int 950 @param contourWidth: width of the overlaid contours 951 952 """ 953 954 cut=intensityCutImage(backgroundImageData, cutLevels) 955 956 # Make plot of just the background image 957 aspectR=float(cut['image'].shape[0])/float(cut['image'].shape[1]) 958 pylab.figure(figsize=(10,10*aspectR)) 959 pylab.axes([0,0,1,1]) 960 961 try: 962 colorMap=pylab.cm.get_cmap(colorMapName) 963 except AssertionError: 964 raise Exception(colorMapName+" is not a defined matplotlib colormap.") 965 966 if cutLevels[0]=="histEq": 967 pylab.imshow(cut['image'], interpolation="bilinear", origin='lower', cmap=colorMap) 968 969 else: 970 pylab.imshow(cut['image'], interpolation="bilinear", norm=cut['norm'], origin='lower', 971 cmap=colorMap) 972 973 pylab.axis("off") 974 975 # Add the contours 976 contourData=generateContourOverlay(backgroundImageData, backgroundImageWCS, contourImageData, \ 977 contourImageWCS, contourLevels, contourSmoothFactor) 978 979 pylab.contour(contourData['scaledImage'], contourData['contourLevels'], colors=contourColor, 980 linewidths=contourWidth) 981 982 pylab.savefig("out_astImages.png") 983 pylab.close("all") 984 985 try: 986 from PIL import Image 987 except: 988 raise Exception("astImages.saveContourOverlayBitmap requires the Python Imaging Library to be installed") 989 990 im=Image.open("out_astImages.png") 991 im.thumbnail((int(size),int(size))) 992 im.save(outputFileName) 993 994 os.remove("out_astImages.png")
995 996 #---------------------------------------------------------------------------------------------------
997 -def saveFITS(outputFileName, imageData, imageWCS = None):
998 """Writes an image array to a new .fits file. 999 1000 @type outputFileName: string 1001 @param outputFileName: filename of output FITS image 1002 @type imageData: numpy array 1003 @param imageData: image data array 1004 @type imageWCS: astWCS.WCS object 1005 @param imageWCS: image WCS object 1006 1007 @note: If imageWCS=None, the FITS image will be written with a rudimentary header containing 1008 no meta data. 1009 1010 """ 1011 1012 if os.path.exists(outputFileName): 1013 os.remove(outputFileName) 1014 1015 newImg=pyfits.HDUList() 1016 1017 if imageWCS!=None: 1018 hdu=pyfits.PrimaryHDU(None, imageWCS.header) 1019 else: 1020 hdu=pyfits.PrimaryHDU(None, None) 1021 1022 hdu.data=imageData 1023 newImg.append(hdu) 1024 newImg.writeto(outputFileName) 1025 newImg.close()
1026 1027 #---------------------------------------------------------------------------------------------------
1028 -def histEq(inputArray, numBins):
1029 """Performs histogram equalisation of the input numpy array. 1030 1031 @type inputArray: numpy array 1032 @param inputArray: image data array 1033 @type numBins: int 1034 @param numBins: number of bins in which to perform the operation (e.g. 1024) 1035 @rtype: numpy array 1036 @return: image data array 1037 1038 """ 1039 1040 imageData=inputArray 1041 1042 # histogram equalisation: we want an equal number of pixels in each intensity range 1043 sortedDataIntensities=numpy.sort(numpy.ravel(imageData)) 1044 median=numpy.median(sortedDataIntensities) 1045 1046 # Make cumulative histogram of data values, simple min-max used to set bin sizes and range 1047 dataCumHist=numpy.zeros(numBins) 1048 minIntensity=sortedDataIntensities.min() 1049 maxIntensity=sortedDataIntensities.max() 1050 histRange=maxIntensity-minIntensity 1051 binWidth=histRange/float(numBins-1) 1052 for i in range(len(sortedDataIntensities)): 1053 binNumber=int(math.ceil((sortedDataIntensities[i]-minIntensity)/binWidth)) 1054 addArray=numpy.zeros(numBins) 1055 onesArray=numpy.ones(numBins-binNumber) 1056 onesRange=list(range(binNumber, numBins)) 1057 numpy.put(addArray, onesRange, onesArray) 1058 dataCumHist=dataCumHist+addArray 1059 1060 # Make ideal cumulative histogram 1061 idealValue=dataCumHist.max()/float(numBins) 1062 idealCumHist=numpy.arange(idealValue, dataCumHist.max()+idealValue, idealValue) 1063 1064 # Map the data to the ideal 1065 for y in range(imageData.shape[0]): 1066 for x in range(imageData.shape[1]): 1067 # Get index corresponding to dataIntensity 1068 intensityBin=int(math.ceil((imageData[y][x]-minIntensity)/binWidth)) 1069 1070 # Guard against rounding errors (happens rarely I think) 1071 if intensityBin<0: 1072 intensityBin=0 1073 if intensityBin>len(dataCumHist)-1: 1074 intensityBin=len(dataCumHist)-1 1075 1076 # Get the cumulative frequency corresponding intensity level in the data 1077 dataCumFreq=dataCumHist[intensityBin] 1078 1079 # Get the index of the corresponding ideal cumulative frequency 1080 idealBin=numpy.searchsorted(idealCumHist, dataCumFreq) 1081 idealIntensity=(idealBin*binWidth)+minIntensity 1082 imageData[y][x]=idealIntensity 1083 1084 return imageData
1085 1086 #---------------------------------------------------------------------------------------------------
1087 -def normalise(inputArray, clipMinMax):
1088 """Clips the inputArray in intensity and normalises the array such that minimum and maximum 1089 values are 0, 1. Clip in intensity is specified by clipMinMax, a list in the format 1090 [clipMin, clipMax] 1091 1092 Used for normalising image arrays so that they can be turned into RGB arrays that matplotlib 1093 can plot (see L{astPlots.ImagePlot}). 1094 1095 @type inputArray: numpy array 1096 @param inputArray: image data array 1097 @type clipMinMax: list 1098 @param clipMinMax: [minimum value of clipped array, maximum value of clipped array] 1099 @rtype: numpy array 1100 @return: normalised array with minimum value 0, maximum value 1 1101 1102 """ 1103 clipped=inputArray.clip(clipMinMax[0], clipMinMax[1]) 1104 slope=1.0/(clipMinMax[1]-clipMinMax[0]) 1105 intercept=-clipMinMax[0]*slope 1106 clipped=clipped*slope+intercept 1107 1108 return clipped
1109