From 31d36e03b23f0fa0affc88cd7654be8bc0219c59 Mon Sep 17 00:00:00 2001
From: Andreas Specovius <andreas.specovius@fau.de>
Date: Thu, 8 Nov 2018 14:54:15 +0100
Subject: [PATCH] Piping for ctskymap implemented.

---
 pyext/ctskymap.i          |  3 +++
 src/ctskymap/ctskymap.cpp | 83 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------
 src/ctskymap/ctskymap.hpp | 16 ++++++++++++++++
 3 files changed, 84 insertions(+), 18 deletions(-)

diff --git a/pyext/ctskymap.i b/pyext/ctskymap.i
index 15cc9f5..b62bf11 100755
--- a/pyext/ctskymap.i
+++ b/pyext/ctskymap.i
@@ -48,6 +48,9 @@ public:
     void           run(void);
     void           save(void);
     void           publish(const std::string& name = "");
+
+    void           exclmap(const GSkyMap& exclusionmap);
+    const GSkyMap& exclmap(void) const;
     const GSkyMap& skymap(void) const;
 };
 
diff --git a/src/ctskymap/ctskymap.cpp b/src/ctskymap/ctskymap.cpp
index b3b4f79..1ba2ed6 100755
--- a/src/ctskymap/ctskymap.cpp
+++ b/src/ctskymap/ctskymap.cpp
@@ -337,6 +337,21 @@ void ctskymap::publish(const std::string& name)
 }
 
 
+/***********************************************************************//**
+ * @brief Set the exclusion map
+ *
+ * @param[in] exclmap Exclusion map object.
+ ***************************************************************************/
+void ctskymap::exclmap(const GSkyMap& exclusionmap)
+{
+    // Assign exclusion map
+    m_exclmap = exclusionmap;
+
+    // Return
+    return;
+}
+
+
 /*==========================================================================
  =                                                                         =
  =                             Private methods                             =
@@ -522,6 +537,7 @@ void ctskymap::get_parameters(void)
         if (m_iterations > 0) {
             m_threshold  = (*this)["threshold"].real();
         }
+
         (*this)["inexclusion"].query();
         m_usefft      = (*this)["usefft"].boolean();
 
@@ -576,7 +592,12 @@ void ctskymap::setup_maps(void)
     if (m_bkgsubtract != "NONE") {
 
         // Setup the exclusions map
-        setup_exclusion_map();
+        if (m_exclmap.is_empty()) {
+            setup_exclusion_map();
+        }
+        else {
+            adjust_exclusion_map();
+        }
 
         // Cache the pixel solid angles and sky directions
         m_solidangle.reserve(m_counts.npix());
@@ -648,24 +669,11 @@ void ctskymap::setup_exclusion_map(void)
  ***************************************************************************/
 void ctskymap::setup_exclusion_map_fits(const GFilename& filename)
 {
-    // Load the fits image
-    GSkyMap inmap(filename);
-
-    // Loop through the individual pixels in the exclusion map
-    for (int i = 0; i < m_exclmap.npix(); ++i) {
-
-        // Get the pixel direction
-        GSkyDir dir = m_exclmap.pix2dir(i);
+    // Load the fits image to the exclusion map
+    m_exclmap = GSkyMap(filename);
 
-        // Check this sky position in the fits map
-        if (inmap.contains(dir) && (inmap(inmap.dir2inx(dir)) != 0.0)) {
-
-            // Set the pixel to 1
-            m_exclmap(i) = 1.0;
-
-        } // endif: pixel,region overlap check
-
-    } // endfor: looped over exclusion map pixels
+    // Verify exclusion map
+    adjust_exclusion_map();
 
     // Return
     return;
@@ -704,6 +712,45 @@ void ctskymap::setup_exclusion_map_region(const GFilename& filename)
 
 
 /***********************************************************************//**
+ * @brief Verifys that exclusion map points in fov of counts
+ *
+ * Adjust exclusion map. Reinitialise exclusion map via the counts object 
+ * and transfer the masked areas provided by the current exclusion map to the
+ * updated fov.
+ ***************************************************************************/
+void ctskymap::adjust_exclusion_map(void)
+{
+    // Initialise working map
+    GSkyMap dest_map = m_counts;
+
+    // Set all pixels to 0 (no pixel excluded)
+    dest_map = 0.0;
+
+    // Loop through the individual pixels in the exclusion map
+    for (int i = 0; i < dest_map.npix(); ++i) {
+
+        // Get the pixel direction
+        GSkyDir dir = dest_map.pix2dir(i);
+
+        // Check this sky position in the exclusion map
+        if (m_exclmap.contains(dir) && (m_exclmap(m_exclmap.dir2inx(dir)) != 0.0)) {
+
+            // Set the pixel to 1
+            dest_map(i) = 1.0;
+
+        } // endif: pixel,region overlap check
+
+    } // endfor: looped over exclusion map pixels
+
+    // Set exclusion map
+    m_exclmap = dest_map;
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
  * @brief Fill maps from observation container
  *
  * Fill counts and optionally acceptance maps from data in the observation
diff --git a/src/ctskymap/ctskymap.hpp b/src/ctskymap/ctskymap.hpp
index c30b265..d5bdd48 100755
--- a/src/ctskymap/ctskymap.hpp
+++ b/src/ctskymap/ctskymap.hpp
@@ -60,6 +60,9 @@ public:
     void           run(void);
     void           save(void);
     void           publish(const std::string& name = "");
+
+    void           exclmap(const GSkyMap& exclusionmap);
+    const GSkyMap& exclmap(void) const;
     const GSkyMap& skymap(void) const;
 
 protected:
@@ -72,6 +75,7 @@ protected:
     void     setup_exclusion_map(void);
     void     setup_exclusion_map_fits(const GFilename& filename);
     void     setup_exclusion_map_region(const GFilename& filename);
+    void     adjust_exclusion_map(void);
     void     fill_maps(void);
     void     fill_maps_counts(GCTAObservation* obs);
     void     fill_maps_acceptance(GCTAObservation* obs);
@@ -131,6 +135,18 @@ protected:
 
 
 /***********************************************************************//**
+ * @brief Return exclusion map
+ *
+ * @return Reference to exclusion map
+ ***************************************************************************/
+inline
+const GSkyMap& ctskymap::exclmap(void) const
+{
+    return m_exclmap;
+}
+
+
+/***********************************************************************//**
  * @brief Return observation container
  *
  * @return Reference to observation container
--
libgit2 0.21.4