Commit 6ed2e4de authored by Simon Spannagel's avatar Simon Spannagel
Browse files

Move ROI functionality into Detector class

parent ee0b5a93
Pipeline #434421 passed with stages
in 4 minutes and 10 seconds
......@@ -44,6 +44,7 @@ Detector::Detector(const Configuration& config) {
m_nPixelsX = npixels.x();
m_nPixelsY = npixels.y();
m_timingOffset = config.get<double>("time_offset", 0.0);
m_roi = config.getMatrix<int>("roi", std::vector<std::vector<int>>());
this->initialise();
......@@ -172,6 +173,8 @@ Configuration Detector::getConfiguration() {
config.set("mask_file", m_maskfile_name);
}
config.setMatrix("roi", m_roi);
return config;
}
......@@ -271,3 +274,81 @@ double Detector::inPixelY(const PositionVector3D<Cartesian3D<double>> localPosit
double inPixelY = m_pitch.Y() * (row - floor(row));
return inPixelY;
}
// Check if track position is within ROI:
bool Detector::isWithinROI(const Track* track) {
// Check that track is within region of interest using winding number algorithm
auto localIntercept = this->getLocalIntercept(track);
auto coordinates = std::make_pair(this->getColumn(localIntercept), this->getRow(localIntercept));
if(winding_number(coordinates, m_roi) != 0) {
return true;
}
// Outside ROI:
return false;
}
// Check if cluster is within ROI and/or touches ROI border:
bool Detector::isWithinROI(const Cluster* cluster) {}
/* 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 Detector::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 Detector::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;
}
......@@ -103,6 +103,9 @@ namespace corryvreckan {
ROOT::Math::XYZPoint localToGlobal(ROOT::Math::XYZPoint local) { return m_localToGlobal * local; };
ROOT::Math::XYZPoint globalToLocal(ROOT::Math::XYZPoint global) { return m_globalToLocal * global; };
bool isWithinROI(const Track* track);
bool isWithinROI(const Cluster* cluster);
private:
// Member variables
// Detector information
......@@ -113,6 +116,10 @@ namespace corryvreckan {
int m_nPixelsY;
double m_timingOffset;
std::vector<std::vector<int>> m_roi;
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);
// Displacement and rotation in x,y,z
ROOT::Math::XYZPoint m_displacement;
ROOT::Math::XYZVector m_orientation;
......
# 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})
## 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`.
#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(true);
// 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);
} else {
LOG(DEBUG) << "Track within DUT ROI";
}
}
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;
}
#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
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment