Add new Module: RegionOfInterest

has the same functionality as previously implemented in CLICpix2Analysis

src/modules/CLICpix2Analysis/README.md

@@ -9,7 +9,6 @@ Analysis module for CLICpix2 prototypes. This module is still work in progress,
 * `timeCutFrameEdge`: Parameter to discard telescope tracks at the frame edges (start and end of the current CLICpix2 frame). Defaults to `20ns`.
 * `spatialCut`: Spatial cut for associating a track with a DUT cluster, defaults to `50um`.
 * `chi2ndofCut`: Acceptance criterion for telescope tracks, defaults to a value of `3`.
-* `roi`: A region of interest, given as polygon by specifying a matrix of pixel coordinates. Defaults to no ROI defined.
 
 #### Plots produced
 * 2D Map of associated cluster positions
@@ -34,8 +33,6 @@ Analysis module for CLICpix2 prototypes. This module is still work in progress,
 ```toml
 [CLICpix2Analysis]
 timeCutFrameEdge = 50ns
-# Define a rectangle as region of interest:
-roi = [1, 1], [1, 120], [60, 120], [60, 1]
 ```
 
 Parameters to be used in multiple algorithms can also be defined globally at the top of the configuration file. This is highly encouraged for parameters such as `DUT` and `reference`.

src/modules/RegionOfInterest/CMakeLists.txt

+# Define module and return the generated name as MODULE_NAME
+CORRYVRECKAN_MODULE(MODULE_NAME)
+
+# Add source files to library
+CORRYVRECKAN_MODULE_SOURCES(${MODULE_NAME}
+    RegionOfInterest.cpp
+    # ADD SOURCE FILES HERE...
+)
+
+# Provide standard install target
+CORRYVRECKAN_MODULE_INSTALL(${MODULE_NAME})

src/modules/RegionOfInterest/README.md

+## RegionOfInterest
+**Maintainer**: Simon Spannagel (<simon.spannagel@cern.ch>)  
+**Status**: Functional
+
+#### Description
+This module allows marking tracks to be within a region of interest (ROI) on the detector marked as `DUT`. This information can be used in later analyses to restrict the selection of tracks to certain regions of the DUT.
+
+The ROI is defined as a polynomial in local pixle coordinates of the DUT. E.g. a rectangle could be defined by providing the four corners of the shape via
+
+```toml
+roi = [1, 1], [1, 120], [60, 120], [60, 1]
+```
+
+If this module is used but no `roi` parameter is provided in the configuration, all tracks which have an intercept with the DUT are considered to be within the ROI.
+If this module is not used, all tracks default to be outside the ROI.
+
+#### Parameters
+* `roi`: A region of interest, given as polygon by specifying a matrix of pixel coordinates. Defaults to no ROI defined.
+
+#### Plots produced
+No plots are produced.
+
+#### Usage
+```toml
+[RegionOfInterest]
+# Define a rectangle as region of interest:
+roi = [1, 1], [1, 120], [60, 120], [60, 1]
+```
+Parameters to be used in multiple algorithms can also be defined globally at the top of the configuration file. This is highly encouraged for parameters such as `DUT` and `reference`.

src/modules/RegionOfInterest/RegionOfInterest.cpp

+#include "RegionOfInterest.h"
+
+using namespace corryvreckan;
+using namespace std;
+
+RegionOfInterest::RegionOfInterest(Configuration config, std::vector<Detector*> detectors)
+    : Module(std::move(config), std::move(detectors)) {
+    m_DUT = m_config.get<std::string>("DUT");
+    m_roi = m_config.getMatrix<int>("roi", std::vector<std::vector<int>>());
+}
+
+StatusCode RegionOfInterest::run(Clipboard* clipboard) {
+
+    // Get the telescope tracks from the clipboard
+    Tracks* tracks = (Tracks*)clipboard->get("tracks");
+    if(tracks == NULL) {
+        LOG(DEBUG) << "No tracks on the clipboard";
+        return Success;
+    }
+
+    // Get the DUT detector:
+    auto detector = get_detector(m_DUT);
+
+    // Loop over all tracks
+    for(auto& track : (*tracks)) {
+        // Check if it intercepts the DUT
+        auto globalIntercept = detector->getIntercept(track);
+        if(!detector->hasIntercept(track, 1.)) {
+            LOG(DEBUG) << "Track outside DUT";
+            track->set_within_roi(false);
+            continue;
+        }
+
+        // If it does, we define the track as being inside the ROI
+        track->set_within_roi(false);
+
+        // Check that track is within region of interest using winding number algorithm
+        auto localIntercept = detector->getLocalIntercept(track);
+        auto coordinates = std::make_pair(detector->getColumn(localIntercept), detector->getRow(localIntercept));
+        if(winding_number(coordinates, m_roi) == 0) {
+            LOG(DEBUG) << "Track outside DUT ROI";
+            track->set_within_roi(false);
+        }
+    }
+
+    return Success;
+}
+
+/* isLeft(): tests if a point is Left|On|Right of an infinite line.
+ * via: http://geomalgorithms.com/a03-_inclusion.html
+ *    Input:  three points P0, P1, and P2
+ *    Return: >0 for P2 left of the line through P0 and P1
+ *            =0 for P2  on the line
+ *            <0 for P2  right of the line
+ *    See: Algorithm 1 "Area of Triangles and Polygons"
+ */
+int RegionOfInterest::isLeft(std::pair<int, int> pt0, std::pair<int, int> pt1, std::pair<int, int> pt2) {
+    return ((pt1.first - pt0.first) * (pt2.second - pt0.second) - (pt2.first - pt0.first) * (pt1.second - pt0.second));
+}
+
+/* Winding number test for a point in a polygon
+ * via: http://geomalgorithms.com/a03-_inclusion.html
+ *      Input:   x, y = a point,
+ *               polygon = vector of vertex points of a polygon V[n+1] with V[n]=V[0]
+ *      Return:  wn = the winding number (=0 only when P is outside)
+ */
+int RegionOfInterest::winding_number(std::pair<int, int> probe, std::vector<std::vector<int>> polygon) {
+    // Two points don't make an area
+    if(polygon.size() < 3) {
+        LOG(DEBUG) << "No ROI given.";
+        return 0;
+    }
+
+    int wn = 0; // the  winding number counter
+
+    // loop through all edges of the polygon
+
+    // edge from V[i] to  V[i+1]
+    for(int i = 0; i < polygon.size(); i++) {
+        auto point_this = std::make_pair(polygon.at(i).at(0), polygon.at(i).at(1));
+        auto point_next = (i + 1 < polygon.size() ? std::make_pair(polygon.at(i + 1).at(0), polygon.at(i + 1).at(1))
+                                                  : std::make_pair(polygon.at(0).at(0), polygon.at(0).at(1)));
+
+        // start y <= P.y
+        if(point_this.second <= probe.second) {
+            // an upward crossing
+            if(point_next.second > probe.second) {
+                // P left of  edge
+                if(isLeft(point_this, point_next, probe) > 0) {
+                    // have  a valid up intersect
+                    ++wn;
+                }
+            }
+        } else {
+            // start y > P.y (no test needed)
+
+            // a downward crossing
+            if(point_next.second <= probe.second) {
+                // P right of  edge
+                if(isLeft(point_this, point_next, probe) < 0) {
+                    // have  a valid down intersect
+                    --wn;
+                }
+            }
+        }
+    }
+    return wn;
+}

src/modules/RegionOfInterest/RegionOfInterest.h

+#ifndef RegionOfInterest_H
+#define RegionOfInterest_H 1
+
+#include <iostream>
+#include "TCanvas.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "core/module/Module.hpp"
+#include "objects/Cluster.h"
+#include "objects/Pixel.h"
+#include "objects/Track.h"
+
+namespace corryvreckan {
+    /** @ingroup Modules
+     */
+    class RegionOfInterest : public Module {
+
+    public:
+        // Constructors and destructors
+        RegionOfInterest(Configuration config, std::vector<Detector*> detectors);
+        ~RegionOfInterest() {}
+
+        StatusCode run(Clipboard* clipboard);
+
+    private:
+        static int winding_number(std::pair<int, int> probe, std::vector<std::vector<int>> polygon);
+        inline static int isLeft(std::pair<int, int> pt0, std::pair<int, int> pt1, std::pair<int, int> pt2);
+
+        std::vector<std::vector<int>> m_roi;
+        std::string m_DUT;
+    };
+} // namespace corryvreckan
+#endif // RegionOfInterest_H