1# This file is part of the Astrometry.net suite. 2# Licensed under a 3-clause BSD style license - see LICENSE 3try: 4 import pyfits 5except ImportError: 6 try: 7 from astropy.io import fits as pyfits 8 except ImportError: 9 raise ImportError("Cannot import either pyfits or astropy.io.fits") 10import math 11from math import exp 12from matplotlib.pylab import imread 13from numpy.oldnumeric.functions import zeros, ravel 14 15I=imread('3.png') 16I=I[:,:,:3] 17(h,w,planes) = I.shape 18XY = pyfits.open('16b.fits')[1].data 19X = XY.field('X') 20Y = XY.field('Y') 21 22psfw = 1.0 23stars = zeros((h,w)).astype(float) 24 25for (x,y) in zip(X,Y): 26 ix = int(round(x)) 27 iy = int(round(y)) 28 for dy in range(-5, 6): 29 yy = iy + dy 30 if yy < 0 or yy >= h: 31 continue 32 for dx in range(-5, 6): 33 xx = ix + dx 34 if xx < 0 or xx >= w: 35 continue 36 dd = (xx - x)**2 + (yy - y)**2 37 stars[yy,xx] += exp(-dd / (2 * psfw**2)) #1./(psfw**2 * 2 * math.pi 38 39#origfrac = 0.5 40#maxorig = I.max() 41#starfrac = (1.0 - origfrac) + (1.0 - maxorig) 42#for p in range(planes): 43# I[:,:,p] = I[:,:,p] * origfrac + stars/stars.max() * starfrac 44 45for p in range(planes): 46 I[:,:,p] = I[:,:,p] * 0.7 + stars/stars.max() * 0.8 47 48f=open('out.ppm', 'wb') 49f.write('P6 %i %i %i\n' % (w, h, 255)) 50#for j in range(h): 51# for i in range(w): 52# for p in range(planes): 53# f.write(chr(int(round(I[j,i,p] * 255.0)))) 54flatI = (I.ravel() * 255.0).round().astype(int) 55f.write("".join([chr(min(i,255)) for i in flatI])) 56f.close() 57