1#!/usr/local/bin/python3.8 2# This file is part of the Astrometry.net suite. 3# Licensed under a 3-clause BSD style license - see LICENSE 4from __future__ import print_function 5from __future__ import absolute_import 6import matplotlib 7matplotlib.use('Agg') 8import os 9import sys 10import logging 11 12from pylab import * 13from optparse import OptionParser 14 15from astrometry.util.removelines import hist_remove_lines 16from astrometry.util.fits import pyfits_writeto 17 18if __name__ == '__main__': 19 try: 20 try: 21 import pyfits 22 except ImportError: 23 try: 24 from astropy.io import fits as pyfits 25 except ImportError: 26 raise ImportError( 27 "Cannot import either pyfits or astropy.io.fits" 28 ) 29 import astrometry 30 from astrometry.util.shell import shell_escape 31 from astrometry.util.filetype import filetype_short 32 except ImportError: 33 me = sys.argv[0] 34 #print 'i am', me 35 path = os.path.realpath(me) 36 #print 'my real path is', path 37 utildir = os.path.dirname(path) 38 assert(os.path.basename(utildir) == 'util') 39 andir = os.path.dirname(utildir) 40 #assert(os.path.basename(andir) == 'astrometry') 41 rootdir = os.path.dirname(andir) 42 #print 'adding path', rootdir 43 #sys.path += [rootdir] 44 sys.path.insert(1, andir) 45 sys.path.insert(2, rootdir) 46 47import numpy 48try: 49 import pyfits 50except ImportError: 51 try: 52 from astropy.io import fits as pyfits 53 except ImportError: 54 raise ImportError("Cannot import either pyfits or astropy.io.fits") 55from numpy import * 56from numpy.random import rand 57 58''' 59# Returns a numpy array of booleans: True for points that should be kept (are not part of lines) 60if False: 61 # def hist_remove_lines(x, binwidth, binoffset, nsig, plots=False): 62 bins = -binoffset + arange(0, max(x)+binwidth, binwidth) 63 (counts, thebins) = histogram(x, bins) 64 65 # We're ignoring empty bins. 66 occupied = nonzero(counts)[0] 67 noccupied = len(occupied) 68 #k = (counts[occupied] - 1) 69 #mean = sum(k) / float(noccupied) 70 k = counts[occupied] 71 mean = sum(k) / ((max(x) - min(x)) / binwidth) * sqrt(2.) 72 thresh = mean + nsig * sqrt(mean) 73 74 if plots: 75 hist(x, bins) 76 axhline(mean) 77 axhline(thresh, color='r') 78 79 print 'mean', mean, 'thresh:', thresh, 'max:', max(k) 80 81 #logpoisson = k*log(mean) - mean - array([sum(log(1 + arange(kk))) for kk in k]) 82 #uk = unique(k) 83 #ulogpoisson = uk*log(mean) - mean - array([sum(1+arange(kk)) for kk in uk]) 84 #print 85 #for (uuk,ull) in zip(uk,ulogpoisson): 86 # print uuk,ull 87 88 #badbins = occupied[logpoisson < logcut] 89 badbins = occupied[k > thresh] 90 if len(badbins) == 0: 91 return array([True] * len(x)) 92 93 badleft = bins[badbins] 94 badright = badleft + binwidth 95 96 badpoints = sum(array([(x >= L)*(x < R) for (L,R) in zip(badleft, badright)]), 0) 97 return (badpoints == 0) 98''' 99#''' 100 101def removelines(infile, outfile, xcol='X', ycol='Y', plots=False, cut=None, **kwargs): 102 p = pyfits.open(infile) 103 xy = p[1].data 104 hdr = p[1].header 105 x = xy.field(xcol) 106 y = xy.field(ycol) 107 108 NX = max(x) - min(x) 109 NY = max(y) - min(y) 110 nangle = int(ceil(sqrt(NX*NY)/4.)) 111 112 if cut is None: 113 cut = 20 114 115 if plots: 116 clf() 117 plot(x, y, 'r.') 118 119 I = array([True]*len(x)) 120 for i,angle in enumerate(0.75 + linspace(0, pi/2., nangle, endpoint=False)): 121 cost = cos(angle) 122 sint = sin(angle) 123 xx = x*cost + y*sint 124 yy = x*-sint + y*cost 125 xx -= min(xx) 126 yy -= min(yy) 127 128 if plots: 129 print() 130 clf() 131 subplot(2,2,1) 132 plot(xx, yy, 'r.') 133 subplot(2,2,3) 134 135 136 #ix = hist_remove_lines(xx, 0.5, 0.5, 5, plots=plots) 137 ix = hist_remove_lines(xx, 1, 0.5, logcut=-cut) 138 139 if plots: 140 subplot(2,2,4) 141 142 #iy = hist_remove_lines(yy, 0.5, 0.5, 5, plots=plots) 143 iy = hist_remove_lines(yy, 1, 0.5, logcut=-cut) 144 145 I *= ix * iy 146 147 removed = (ix * iy == False) 148 if plots: 149 if sum(removed): 150 plot([min(x[removed]), max(x[removed])], 151 [min(y[removed]), max(y[removed])], 'k-', alpha=0.5) 152 subplot(2,2,1) 153 plot(xx[removed], yy[removed], 'b-', alpha=0.5) 154 plot(xx[removed], yy[removed], 'b.') 155 savefig('rot-%04i.png' % i) 156 157 print('angle', angle, 'removed', (len(x) - sum(ix*iy))) 158 159 xc = x[I] 160 yc = y[I] 161 print('removelines.py: Removed %i sources' % (len(x) - len(xc))) 162 163 if plots: 164 plot(xc, yc, 'o', mec='r', mfc='none') 165 # axes('equal') 166 savefig('after.png') 167 168 p[1].header.add_history('This xylist was filtered by the "removelines_rotate.py" program') 169 p[1].header.add_history('to remove lines of sources') 170 p[1].header.update('REMLINEN', len(x) - len(xc), 'Number of sources removed by "removelines"') 171 172 p[1].data = p[1].data[I] 173 pyfits_writeto(p, outfile) 174 175 return 0 176 177if __name__ == '__main__': 178 parser = OptionParser(usage='%prog <input-file> <output-file>') 179 parser.add_option('-X', '--x-column', dest='xcol', help='X column name') 180 parser.add_option('-Y', '--y-column', dest='ycol', help='Y column name') 181 parser.set_defaults(xcol='X', ycol='Y') 182 opt,args = parser.parse_args() 183 if len(args) != 2: 184 parser.print_help() 185 sys.exit(0) 186 187 kwargs = {} 188 for k in 'xcol','ycol': 189 kwargs[k] = getattr(opt, k) 190 191 infile = args[0] 192 outfile = args[1] 193 rtncode = removelines(infile, outfile, **kwargs) 194 sys.exit(rtncode) 195 196