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