From a2619925b4cce4f113ab3a1f89dbfb2e82f2701d Mon Sep 17 00:00:00 2001
From: Josh Cardenzana <jcardenzana@irap.omp.eu>
Date: Mon, 20 Aug 2018 17:11:25 +0200
Subject: [PATCH] Add computation of pixel by pixel mean and stddev. Add sources by pixel with maximum stddevs above its respective mean Maximum is removed by setting pixel values to mean associated with peak pixel.

---
 cscripts/cssrcdetect.py | 63 +++++++++++++++++++++++++++++----------------------------------
 1 file changed, 29 insertions(+), 34 deletions(-)

diff --git a/cscripts/cssrcdetect.py b/cscripts/cssrcdetect.py
index e70c1e0..f804c1f 100755
--- a/cscripts/cssrcdetect.py
+++ b/cscripts/cssrcdetect.py
@@ -72,13 +72,17 @@ class cssrcdetect(ctools.cscript):
         self['exclrad'].real()
         self['fit_pos'].boolean()
         self['fit_shape'].boolean()
+
+        # Query the smoothing parameters
         self['smoothkernel'].string()
-        self['smoothparam'].real()
+        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)
@@ -102,20 +106,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:
@@ -131,7 +128,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
@@ -147,7 +144,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
 
@@ -155,8 +152,6 @@ class cssrcdetect(ctools.cscript):
         ----------
         skymap : `~gammalib.GSkyMap()`
             Sky map
-        threshold : float, optional
-            Threshold for maximum value
 
         Returns
         -------
@@ -164,7 +159,7 @@ 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
@@ -208,7 +203,7 @@ class cssrcdetect(ctools.cscript):
         # Return copy of map
         return skymap_copy
 
-    def _map_moments(self, skymap):
+    def _map_moments(self, skymap, radius):
         """
         Determine moments of sky map pixels
 
@@ -216,26 +211,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
-
-        # Loop over all pixels
-        for i in range(skymap.npix()):
-            mean += skymap[i]
-            std  += skymap[i]*skymap[i]
-
-        # Compute mean and standard deviation
-        mean /= float(skymap.npix())
-        std  /= float(skymap.npix())
-        std  -= mean * mean
-        std   = math.sqrt(std)
+        # Copy the input skymap
+        mean = gammalib.GSkyMap(skymap)
+        std  = gammalib.GSkyMap(skymap)
+        std *= std
+
+        # Convolve by disk to get bin-by-bin mean
+        mean.smooth("DISK", radius)
+        std.smooth("DISK", radius)
+
+        # Compute the standard deviation for each pixel
+        std = std - (mean*mean)
+        std = std.sqrt()
 
         # Return mean and standard deviation
         return mean, std
--
libgit2 0.21.4