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