99import utils
1010import numpy
1111import math
12+ from astLib .astWCS import WCS
13+ import pyfits
1214
1315class load (object ):
1416
1517
16- def __init__ (self , imagename , pmodel , nmodel , header , psfname = None , noise = None ,
17- snr_thresh = 40 , local_thresh = 0.4 , high_corr_thresh = 0.5 , negdetec_region = 10 ,
18+ def __init__ (self , imagename , psfname , pmodel , nmodel , local_step = 10 ,
19+ snr_thresh = 40 , high_corr_thresh = 0.5 , negdetec_region = 10 ,
1820 negatives_thresh = 5 , phasecenter_excl_radius = None ,
1921 prefix = None , loglevel = 0 ):
2022
@@ -30,16 +32,9 @@ def __init__(self, imagename, pmodel, nmodel, header, psfname=None, noise=None,
3032
3133 header: The header of the input image
3234
33- noise: float, Default None.
34- The noise of the image.
35-
3635 snr_thresh: float, optional. Default is 40.
3736 Any source with 40 x the minimum SNR.
3837
39- local_thresh: float, optional. Default is 0.4.
40- Sources with local variance > 0.4 * the noise have
41- high local variance.
42-
4338 high_corr_thresh: float, optional. Default is 0.5.
4439 Sources of high PSF correlation have correlation above 0.5.
4540
@@ -70,46 +65,43 @@ def __init__(self, imagename, pmodel, nmodel, header, psfname=None, noise=None,
7065 self .pmodel = pmodel
7166 self .nmodel = nmodel
7267 self .psfname = psfname
73- self .hdr = header
7468
7569 self .log .info (" Loading image data" )
7670
77- self .noise = noise
78- if not self .noise :
79- self .log .info (" No noise value provided."
80- " Setting it to 1e-6. Please provide"
81- " the noise." )
82- self .noise = 1e-6
83-
8471 # tags
8572 self .dd_tag = "dE"
8673
8774 # thresholds
88- self .snr_thresh = snr_thresh
89- self .localthresh = local_thresh
75+ self .snr_factor = snr_thresh
76+ # self.localthresh = local_thresh
9077 self .corrthresh = high_corr_thresh
9178 self .negthresh = negatives_thresh
92-
9379
80+ with pyfits .open (imagename ) as hdu :
81+ self .hdr = hdu [0 ].header
82+ self .data = utils .image_data (hdu [0 ].data )
83+
84+ self .wcs = WCS (self .hdr , mode = "pyfits" )
85+
86+ self .locstep = local_step
87+
9488 #regions
9589 self .phaserad = phasecenter_excl_radius # radius from the phase center
9690 self .negregion = negdetec_region # region to look for negatives
9791
9892 # conversion
99- self .r2d = 180.0 / math .pi
100- self .d2r = math .pi / 180.0
10193 self .bmaj = self .hdr ["BMAJ" ] # in degrees
10294
103- self .ra0 = self .hdr ["CRVAL1" ] * self . d2r
104- self .dec0 = self .hdr ["CRVAL2" ] * self . d2r
105-
95+ self .ra0 = numpy . deg2rad ( self .hdr ["CRVAL1" ])
96+ self .dec0 = numpy . deg2rad ( self .hdr ["CRVAL2" ])
97+
10698
10799 def number_negatives (self , source ):
108100
109101
110102 #sources = filter(lambda src: src.getTag(tag), psources)
111103
112- tolerance = self .negregion * self .bmaj * self . d2r
104+ tolerance = numpy . deg2rad ( self .negregion * self .bmaj )
113105
114106 if self .phaserad :
115107 radius = numpy .deg2rad (self .phaserad * self .bmaj )
@@ -123,26 +115,40 @@ def number_negatives(self, source):
123115 source .setTag (self .dd_tag , True )
124116 else :
125117 source .setTag (self .dd_tag , True )
118+
126119
120+ def local_noise (self , pos ):
121+
122+ # computing the local noise using MAD
123+ x , y = pos
124+
125+ sub_data = self .data [y - self .locstep :y + self .locstep ,
126+ x - self .locstep :x + self .locstep ]
127+ noise = numpy .mean (abs (sub_data - numpy .mean (sub_data )))
128+ return noise
129+
127130
128131 def source_selection (self ):
129132
130133 sources = self .pmodel .sources
131-
132- snr = [ src . flux . I / self . noise for src in sources ]
133-
134- thresh = self .snr_thresh * min ( snr )
135- n = 0
136- for srs , s in zip ( sources , snr ):
137- if s > thresh :
138- local = srs . l
139- if local > self . localthresh * self . noise :
140- if self .psfname :
141- corr = srs . cf
142- if corr > self . corrthresh :
143- self .number_negatives ( srs )
144-
145- if not self .psfname :
134+ noise , mean = utils . negative_noise ( self . data , self . prefix )
135+ for srs in sources :
136+ pos = map ( lambda rad : numpy . rad2deg ( rad ),( srs . pos . ra , srs . pos . dec ))
137+ positions = self .wcs . wcs2pix ( * pos )
138+
139+ # local noise, determining SNR using the local noise
140+ # threshold is determined using the global noise
141+ local_noise = self . local_noise ( positions )
142+ signal_to_noise = srs . flux . I / local_noise
143+ thresh = self .snr_factor * noise
144+
145+ if signal_to_noise > thresh and srs . rel > 0.99 :
146+ if self .psfname :
147+ corr = srs . correlation_factor
148+ if corr > self .corrthresh :
146149 self .number_negatives (srs )
150+
151+ if not self .psfname :
152+ self .number_negatives (srs )
147153
148154 return self .pmodel , self .nmodel
0 commit comments