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