diff --git a/cscripts/cssrcdetect.par b/cscripts/cssrcdetect.par index d3caa9f..d3c9d96 100755 --- a/cscripts/cssrcdetect.par +++ b/cscripts/cssrcdetect.par @@ -23,13 +23,16 @@ outds9file, f, a, ds9.reg,,, "Output DS9 region file" # # Script parameters #================== -srcmodel, s, a, POINT,POINT,, "Source model type" -bkgmodel, s, a, NONE,NONE|IRF|AEFF|CUBE,, "Background model type" -threshold, r, a, 5.0,,, "Detection threshold (Gaussian sigma)" -maxsrcs, i, h, 20,1,1000, "Maximum number of sources" -exclrad, r, h, 0.2,,, "Exclusion radius (deg)" -fit_pos, b, h, yes,,, "Fit source positions?" -fit_shape, b, h, yes,,, "Fit source shapes?" +srcmodel, s, a, POINT,POINT,, "Source model type" +bkgmodel, s, a, NONE,NONE|IRF|AEFF|CUBE,, "Background model type" +threshold, r, a, 5.0,,, "Detection threshold (Gaussian sigma)" +maxsrcs, i, h, 20,1,1000, "Maximum number of sources" +exclrad, r, h, 0.2,,, "Exclusion radius (deg)" +fit_pos, b, h, yes,,, "Fit source positions?" +fit_shape, b, h, yes,,, "Fit source shapes?" +momentradius, r, a, 200,,, "Radius (deg) for computing map mean and std-deviation values" +smoothkernel, s, a, NONE,NONE|DISK|GAUSSIAN,,"Smoothing kernel to use for smoothing the input sky map" +smoothparam, r, a, 0.05,,,"Smoothing kernel parameters (sigma for GAUS, radius for DISK)" # # Standard parameters diff --git a/cscripts/cssrcdetect.py b/cscripts/cssrcdetect.py index a5dbc9d..f1ee6ae 100755 --- a/cscripts/cssrcdetect.py +++ b/cscripts/cssrcdetect.py @@ -51,6 +51,9 @@ class cssrcdetect(ctools.cscript): else: self._map = gammalib.GSkyMap() + # Initialise other members + self._map_dirs = [] + # Return return @@ -73,10 +76,16 @@ class cssrcdetect(ctools.cscript): self['fit_pos'].boolean() self['fit_shape'].boolean() + # Query the smoothing parameters + self['smoothkernel'].string() + if self['smoothkernel'].string().upper() != 'NONE': + self['smoothparam'].real() + # Query ahead output model filename if self._read_ahead(): self['outmodel'].filename() self['outds9file'].filename() + self['momentradius'].real() # Write input parameters into logger self._log_parameters(gammalib.TERSE) @@ -100,20 +109,13 @@ class cssrcdetect(ctools.cscript): self._log_header3(gammalib.NORMAL, 'Iteration '+str(i+1)) # Get map moments - mean, std = self._map_moments(counts) - - # Log map moments - self._log_value(gammalib.NORMAL, 'Map mean', mean) - self._log_value(gammalib.NORMAL, 'Map standard deviation', std) + mean, std = self._map_moments(counts, self['momentradius'].real()) # Compute threshold - threshold = mean + self['threshold'].real() * std - - # Log map threshold - self._log_value(gammalib.NORMAL, 'Map threshold', threshold) - + sigmap = (counts - mean)/std + # Get maximum value and corresponding sky direction - value, pos = self._find_maximum(counts, threshold=threshold) + value, pos = self._find_maximum(sigmap) # If maximum found then log maximum and add model if pos is not None: @@ -129,7 +131,7 @@ class cssrcdetect(ctools.cscript): self._add_model(pos, name) # Remove maximum from map - counts = self._remove_maximum(counts, pos, mean, + counts = self._remove_maximum(counts, pos, mean(pos), radius=self['exclrad'].real()) # ... otherwise log that no maximum was found and break iterations @@ -145,7 +147,7 @@ class cssrcdetect(ctools.cscript): # Return return - def _find_maximum(self, skymap, threshold=0.0): + def _find_maximum(self, skymap): """ Find maximum in a sky map @@ -153,8 +155,6 @@ class cssrcdetect(ctools.cscript): ---------- skymap : `~gammalib.GSkyMap()` Sky map - threshold : float, optional - Threshold for maximum value Returns ------- @@ -162,14 +162,14 @@ class cssrcdetect(ctools.cscript): Maximum sky map value and corresponding sky direction """ # Initialise maximum pixel value and sky direction - value = threshold + value = self['threshold'].real() pos = None # Loop over all pixels and find maximum for i in range(skymap.npix()): if skymap[i] > value: value = skymap[i] - pos = skymap.inx2dir(i) + pos = self._map_dirs[i] # Return sky direction of maximum return value, pos @@ -197,16 +197,19 @@ class cssrcdetect(ctools.cscript): # Copy skymap skymap_copy = skymap.copy() + # Cache the cosine of the radius + cos_radius = math.cos(math.radians(radius)) + # Loop over all pixels and find maximum for i in range(skymap_copy.npix()): - skymap_dir = skymap_copy.inx2dir(i) - if skymap_dir.dist_deg(pos) < radius: + skymap_dir = self._map_dirs[i] + if skymap_dir.cos_dist(pos) > cos_radius: skymap_copy[i] = value # Return copy of map return skymap_copy - def _map_moments(self, skymap): + def _map_moments(self, skymap, radius): """ Determine moments of sky map pixels @@ -214,26 +217,26 @@ class cssrcdetect(ctools.cscript): ---------- skymap : `~gammalib.GSkyMap()` Sky map + radius : float + radius (deg) for pixel consideration Returns ------- - mean, std : tuple of float - Mean and standard deviation of pixel values + mean, std : tuple of GSkyMap + Mean and standard deviation of pixel values within a given radius """ - # Initialise mean and standard deviation - mean = 0.0 - std = 0.0 + # Copy the input skymap + mean = gammalib.GSkyMap(skymap) + std = gammalib.GSkyMap(skymap) + std *= std - # Loop over all pixels - for i in range(skymap.npix()): - mean += skymap[i] - std += skymap[i]*skymap[i] + # Convolve by disk to get bin-by-bin mean + mean.smooth('DISK', radius) + std.smooth('DISK', radius) - # Compute mean and standard deviation - mean /= float(skymap.npix()) - std /= float(skymap.npix()) - std -= mean * mean - std = math.sqrt(std) + # Compute the standard deviation for each pixel + std = std - (mean*mean) + std = std.sqrt() # Return mean and standard deviation return mean, std @@ -338,6 +341,33 @@ class cssrcdetect(ctools.cscript): return model + def _smooth_skymap(self, skymap): + """ + Smooth the input sky map if a valid kernel was supplied + + Parameters + ---------- + skymap : `~gammalib.GSkyMap()` + Sky map + """ + # Log header + self._log_header3(gammalib.NORMAL, 'Smoothing Skymap') + + # Make sure the smoothing kernel is not 'NONE' + if self['smoothkernel'].string().upper() != 'NONE': + self._log_value(gammalib.NORMAL, 'Kernel', + self['smoothkernel'].string()) + self._log_value(gammalib.NORMAL, 'Parameter', + self['smoothparam'].real()) + skymap.smooth(self['smoothkernel'].string(), + self['smoothparam'].real()) + else: + self._log_string(gammalib.NORMAL, + 'Smoothing kernel is NONE, smoothing will be ignored.') + + return + + def _load_skymap(self): """ Load sky map @@ -360,7 +390,18 @@ class cssrcdetect(ctools.cscript): fits.close() # Return - return skymap + return skymap.extract(0) + + + def _cache_map_dirs(self, skymap): + """ + Cache map pixel positions to save some computation time + """ + # Setup the skymap directions + if not skymap.is_empty(): + self._map_dirs = [skymap.inx2dir(i) for i in range(skymap.npix())] + else: + self._map_dirs = [] # Public methods @@ -378,13 +419,18 @@ class cssrcdetect(ctools.cscript): # Get parameters self._get_parameters() - # Write header - self._log_header1(gammalib.NORMAL, 'Source detection') - # If sky map is empty then load it from input parameter if self._map.is_empty(): self._map = self._load_skymap() + self._cache_map_dirs(self._map) + + # Smooth the map + self._smooth_skymap(self._map) + + # Write header + self._log_header1(gammalib.NORMAL, 'Source detection') + # Detect sources self._detect_sources(self._map) diff --git a/doc/source/users/reference_manual/cssrcdetect.rst b/doc/source/users/reference_manual/cssrcdetect.rst index 749a922..be43076 100755 --- a/doc/source/users/reference_manual/cssrcdetect.rst +++ b/doc/source/users/reference_manual/cssrcdetect.rst @@ -48,6 +48,16 @@ General parameters Detection threshold (Gaussian sigma). Only sources above this detection threshold will be added to the model definiton XML file. +``momentradius [real]`` + Radius (deg) for computing map mean and std-deviation values. + +``smoothkernel [string]`` + Smoothing kernel to use for smoothing the input sky map. + +``smoothparam [real]`` + Smoothing kernel parameters (sigma for GAUS, radius for DISK, both in + degrees). + ``(maxsrcs = 20) [integer]`` Maximum number of sources that should be detected. diff --git a/test/test_cssrcdetect.py b/test/test_cssrcdetect.py index 96ee1e9..81db56a 100755 --- a/test/test_cssrcdetect.py +++ b/test/test_cssrcdetect.py @@ -75,6 +75,8 @@ class Test(test): cmd = cssrcdetect+' inmap="'+self._map+'"'+ \ ' outmodel="cssrcdetect_cmd1.xml"'+ \ ' outds9file="cssrcdetect_cmd1.reg"'+ \ + ' momentradius=200.0'+ \ + ' smoothkernel=NONE smoothparam=0.1'+ \ ' srcmodel="POINT" bkgmodel="NONE" threshold=10.0'+ \ ' logfile="cssrcdetect_cmd1.log" chatter=1' @@ -89,6 +91,8 @@ class Test(test): cmd = cssrcdetect+' inmap="skymap_that_does_not_exist"'+ \ ' outmodel="cssrcdetect_cmd1.xml"'+ \ ' outds9file="cssrcdetect_cmd1.reg"'+ \ + ' momentradius=200.0'+ \ + ' smoothkernel=NONE smoothparam=0.1'+ \ ' srcmodel="POINT" bkgmodel="NONE" threshold=10.0'+ \ ' logfile="cssrcdetect_cmd2.log" debug=yes chatter=2' @@ -130,6 +134,9 @@ class Test(test): srcdetect['srcmodel'] = 'POINT' srcdetect['bkgmodel'] = 'NONE' srcdetect['threshold'] = 10.0 + srcdetect['momentradius'] = 200.0 + srcdetect['smoothkernel'] = "NONE" + srcdetect['smoothparam'] = 0.1 srcdetect['logfile'] = 'cssrcdetect_py1.log' srcdetect['chatter'] = 2