Simon Bonnefoy / ctools

Commit 77884dc896d7a33f755fe4da0f76249dd05d6b2f

Authored by Jürgen Knödlseder
2 parents a5cfc4d8 d6d31814

Merge branch '2498-cssrcdetect_local_selection' into integration

Inline Side-by-side
Showing 4 changed files with 111 additions and 45 deletions
cscripts/cssrcdetect.par
  View file @77884dc
... ... @@ -23,13 +23,16 @@ outds9file, f, a, ds9.reg,,, "Output DS9 region file"
23 23 #
24 24 # Script parameters
25 25 #==================
26   -srcmodel, s, a, POINT,POINT,, "Source model type"
27   -bkgmodel, s, a, NONE,NONE|IRF|AEFF|CUBE,, "Background model type"
28   -threshold, r, a, 5.0,,, "Detection threshold (Gaussian sigma)"
29   -maxsrcs, i, h, 20,1,1000, "Maximum number of sources"
30   -exclrad, r, h, 0.2,,, "Exclusion radius (deg)"
31   -fit_pos, b, h, yes,,, "Fit source positions?"
32   -fit_shape, b, h, yes,,, "Fit source shapes?"
  26 +srcmodel, s, a, POINT,POINT,, "Source model type"
  27 +bkgmodel, s, a, NONE,NONE|IRF|AEFF|CUBE,, "Background model type"
  28 +threshold, r, a, 5.0,,, "Detection threshold (Gaussian sigma)"
  29 +maxsrcs, i, h, 20,1,1000, "Maximum number of sources"
  30 +exclrad, r, h, 0.2,,, "Exclusion radius (deg)"
  31 +fit_pos, b, h, yes,,, "Fit source positions?"
  32 +fit_shape, b, h, yes,,, "Fit source shapes?"
  33 +momentradius, r, a, 200,,, "Radius (deg) for computing map mean and std-deviation values"
  34 +smoothkernel, s, a, NONE,NONE|DISK|GAUSSIAN,,"Smoothing kernel to use for smoothing the input sky map"
  35 +smoothparam, r, a, 0.05,,,"Smoothing kernel parameters (sigma for GAUS, radius for DISK)"
33 36  
34 37 #
35 38 # Standard parameters
... ...
cscripts/cssrcdetect.py
  View file @77884dc
... ... @@ -51,6 +51,9 @@ class cssrcdetect(ctools.cscript):
51 51 else:
52 52 self._map = gammalib.GSkyMap()
53 53  
  54 + # Initialise other members
  55 + self._map_dirs = []
  56 +
54 57 # Return
55 58 return
56 59  
... ... @@ -73,10 +76,16 @@ class cssrcdetect(ctools.cscript):
73 76 self['fit_pos'].boolean()
74 77 self['fit_shape'].boolean()
75 78  
  79 + # Query the smoothing parameters
  80 + self['smoothkernel'].string()
  81 + if self['smoothkernel'].string().upper() != 'NONE':
  82 + self['smoothparam'].real()
  83 +
76 84 # Query ahead output model filename
77 85 if self._read_ahead():
78 86 self['outmodel'].filename()
79 87 self['outds9file'].filename()
  88 + self['momentradius'].real()
80 89  
81 90 # Write input parameters into logger
82 91 self._log_parameters(gammalib.TERSE)
... ... @@ -100,20 +109,13 @@ class cssrcdetect(ctools.cscript):
100 109 self._log_header3(gammalib.NORMAL, 'Iteration '+str(i+1))
101 110  
102 111 # Get map moments
103   - mean, std = self._map_moments(counts)
104   -
105   - # Log map moments
106   - self._log_value(gammalib.NORMAL, 'Map mean', mean)
107   - self._log_value(gammalib.NORMAL, 'Map standard deviation', std)
  112 + mean, std = self._map_moments(counts, self['momentradius'].real())
108 113  
109 114 # Compute threshold
110   - threshold = mean + self['threshold'].real() * std
111   -
112   - # Log map threshold
113   - self._log_value(gammalib.NORMAL, 'Map threshold', threshold)
114   -
  115 + sigmap = (counts - mean)/std
  116 +
115 117 # Get maximum value and corresponding sky direction
116   - value, pos = self._find_maximum(counts, threshold=threshold)
  118 + value, pos = self._find_maximum(sigmap)
117 119  
118 120 # If maximum found then log maximum and add model
119 121 if pos is not None:
... ... @@ -129,7 +131,7 @@ class cssrcdetect(ctools.cscript):
129 131 self._add_model(pos, name)
130 132  
131 133 # Remove maximum from map
132   - counts = self._remove_maximum(counts, pos, mean,
  134 + counts = self._remove_maximum(counts, pos, mean(pos),
133 135 radius=self['exclrad'].real())
134 136  
135 137 # ... otherwise log that no maximum was found and break iterations
... ... @@ -145,7 +147,7 @@ class cssrcdetect(ctools.cscript):
145 147 # Return
146 148 return
147 149  
148   - def _find_maximum(self, skymap, threshold=0.0):
  150 + def _find_maximum(self, skymap):
149 151 """
150 152 Find maximum in a sky map
151 153  
... ... @@ -153,8 +155,6 @@ class cssrcdetect(ctools.cscript):
153 155 ----------
154 156 skymap : `~gammalib.GSkyMap()`
155 157 Sky map
156   - threshold : float, optional
157   - Threshold for maximum value
158 158  
159 159 Returns
160 160 -------
... ... @@ -162,14 +162,14 @@ class cssrcdetect(ctools.cscript):
162 162 Maximum sky map value and corresponding sky direction
163 163 """
164 164 # Initialise maximum pixel value and sky direction
165   - value = threshold
  165 + value = self['threshold'].real()
166 166 pos = None
167 167  
168 168 # Loop over all pixels and find maximum
169 169 for i in range(skymap.npix()):
170 170 if skymap[i] > value:
171 171 value = skymap[i]
172   - pos = skymap.inx2dir(i)
  172 + pos = self._map_dirs[i]
173 173  
174 174 # Return sky direction of maximum
175 175 return value, pos
... ... @@ -197,16 +197,19 @@ class cssrcdetect(ctools.cscript):
197 197 # Copy skymap
198 198 skymap_copy = skymap.copy()
199 199  
  200 + # Cache the cosine of the radius
  201 + cos_radius = math.cos(math.radians(radius))
  202 +
200 203 # Loop over all pixels and find maximum
201 204 for i in range(skymap_copy.npix()):
202   - skymap_dir = skymap_copy.inx2dir(i)
203   - if skymap_dir.dist_deg(pos) < radius:
  205 + skymap_dir = self._map_dirs[i]
  206 + if skymap_dir.cos_dist(pos) > cos_radius:
204 207 skymap_copy[i] = value
205 208  
206 209 # Return copy of map
207 210 return skymap_copy
208 211  
209   - def _map_moments(self, skymap):
  212 + def _map_moments(self, skymap, radius):
210 213 """
211 214 Determine moments of sky map pixels
212 215  
... ... @@ -214,26 +217,26 @@ class cssrcdetect(ctools.cscript):
214 217 ----------
215 218 skymap : `~gammalib.GSkyMap()`
216 219 Sky map
  220 + radius : float
  221 + radius (deg) for pixel consideration
217 222  
218 223 Returns
219 224 -------
220   - mean, std : tuple of float
221   - Mean and standard deviation of pixel values
  225 + mean, std : tuple of GSkyMap
  226 + Mean and standard deviation of pixel values within a given radius
222 227 """
223   - # Initialise mean and standard deviation
224   - mean = 0.0
225   - std = 0.0
  228 + # Copy the input skymap
  229 + mean = gammalib.GSkyMap(skymap)
  230 + std = gammalib.GSkyMap(skymap)
  231 + std *= std
226 232  
227   - # Loop over all pixels
228   - for i in range(skymap.npix()):
229   - mean += skymap[i]
230   - std += skymap[i]*skymap[i]
  233 + # Convolve by disk to get bin-by-bin mean
  234 + mean.smooth('DISK', radius)
  235 + std.smooth('DISK', radius)
231 236  
232   - # Compute mean and standard deviation
233   - mean /= float(skymap.npix())
234   - std /= float(skymap.npix())
235   - std -= mean * mean
236   - std = math.sqrt(std)
  237 + # Compute the standard deviation for each pixel
  238 + std = std - (mean*mean)
  239 + std = std.sqrt()
237 240  
238 241 # Return mean and standard deviation
239 242 return mean, std
... ... @@ -338,6 +341,33 @@ class cssrcdetect(ctools.cscript):
338 341 return model
339 342  
340 343  
  344 + def _smooth_skymap(self, skymap):
  345 + """
  346 + Smooth the input sky map if a valid kernel was supplied
  347 +
  348 + Parameters
  349 + ----------
  350 + skymap : `~gammalib.GSkyMap()`
  351 + Sky map
  352 + """
  353 + # Log header
  354 + self._log_header3(gammalib.NORMAL, 'Smoothing Skymap')
  355 +
  356 + # Make sure the smoothing kernel is not 'NONE'
  357 + if self['smoothkernel'].string().upper() != 'NONE':
  358 + self._log_value(gammalib.NORMAL, 'Kernel',
  359 + self['smoothkernel'].string())
  360 + self._log_value(gammalib.NORMAL, 'Parameter',
  361 + self['smoothparam'].real())
  362 + skymap.smooth(self['smoothkernel'].string(),
  363 + self['smoothparam'].real())
  364 + else:
  365 + self._log_string(gammalib.NORMAL,
  366 + 'Smoothing kernel is NONE, smoothing will be ignored.')
  367 +
  368 + return
  369 +
  370 +
341 371 def _load_skymap(self):
342 372 """
343 373 Load sky map
... ... @@ -360,7 +390,18 @@ class cssrcdetect(ctools.cscript):
360 390 fits.close()
361 391  
362 392 # Return
363   - return skymap
  393 + return skymap.extract(0)
  394 +
  395 +
  396 + def _cache_map_dirs(self, skymap):
  397 + """
  398 + Cache map pixel positions to save some computation time
  399 + """
  400 + # Setup the skymap directions
  401 + if not skymap.is_empty():
  402 + self._map_dirs = [skymap.inx2dir(i) for i in range(skymap.npix())]
  403 + else:
  404 + self._map_dirs = []
364 405  
365 406  
366 407 # Public methods
... ... @@ -378,13 +419,18 @@ class cssrcdetect(ctools.cscript):
378 419 # Get parameters
379 420 self._get_parameters()
380 421  
381   - # Write header
382   - self._log_header1(gammalib.NORMAL, 'Source detection')
383   -
384 422 # If sky map is empty then load it from input parameter
385 423 if self._map.is_empty():
386 424 self._map = self._load_skymap()
387 425  
  426 + self._cache_map_dirs(self._map)
  427 +
  428 + # Smooth the map
  429 + self._smooth_skymap(self._map)
  430 +
  431 + # Write header
  432 + self._log_header1(gammalib.NORMAL, 'Source detection')
  433 +
388 434 # Detect sources
389 435 self._detect_sources(self._map)
390 436  
... ...
doc/source/users/reference_manual/cssrcdetect.rst
  View file @77884dc
... ... @@ -48,6 +48,16 @@ General parameters
48 48 Detection threshold (Gaussian sigma). Only sources above this detection
49 49 threshold will be added to the model definiton XML file.
50 50  
  51 +``momentradius [real]``
  52 + Radius (deg) for computing map mean and std-deviation values.
  53 +
  54 +``smoothkernel <NONE|DISK|GAUSSIAN> [string]``
  55 + Smoothing kernel to use for smoothing the input sky map.
  56 +
  57 +``smoothparam [real]``
  58 + Smoothing kernel parameters (sigma for GAUS, radius for DISK, both in
  59 + degrees).
  60 +
51 61 ``(maxsrcs = 20) [integer]``
52 62 Maximum number of sources that should be detected.
53 63  
... ...
test/test_cssrcdetect.py
  View file @77884dc
... ... @@ -75,6 +75,8 @@ class Test(test):
75 75 cmd = cssrcdetect+' inmap="'+self._map+'"'+ \
76 76 ' outmodel="cssrcdetect_cmd1.xml"'+ \
77 77 ' outds9file="cssrcdetect_cmd1.reg"'+ \
  78 + ' momentradius=200.0'+ \
  79 + ' smoothkernel=NONE smoothparam=0.1'+ \
78 80 ' srcmodel="POINT" bkgmodel="NONE" threshold=10.0'+ \
79 81 ' logfile="cssrcdetect_cmd1.log" chatter=1'
80 82  
... ... @@ -89,6 +91,8 @@ class Test(test):
89 91 cmd = cssrcdetect+' inmap="skymap_that_does_not_exist"'+ \
90 92 ' outmodel="cssrcdetect_cmd1.xml"'+ \
91 93 ' outds9file="cssrcdetect_cmd1.reg"'+ \
  94 + ' momentradius=200.0'+ \
  95 + ' smoothkernel=NONE smoothparam=0.1'+ \
92 96 ' srcmodel="POINT" bkgmodel="NONE" threshold=10.0'+ \
93 97 ' logfile="cssrcdetect_cmd2.log" debug=yes chatter=2'
94 98  
... ... @@ -130,6 +134,9 @@ class Test(test):
130 134 srcdetect['srcmodel'] = 'POINT'
131 135 srcdetect['bkgmodel'] = 'NONE'
132 136 srcdetect['threshold'] = 10.0
  137 + srcdetect['momentradius'] = 200.0
  138 + srcdetect['smoothkernel'] = "NONE"
  139 + srcdetect['smoothparam'] = 0.1
133 140 srcdetect['logfile'] = 'cssrcdetect_py1.log'
134 141 srcdetect['chatter'] = 2
135 142  
... ...