epic/src/GeometryHelpers.cpp

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Chao Peng, Whitney Armstrong
0003
0004 #include "GeometryHelpers.h"
0005
0006 // some utility functions that can be shared
0007 namespace epic::geo {
0008
0009 typedef ROOT::Math::XYPoint Point;
0010
0011 // check if a 2d point is already in the container
0012 bool already_placed(const Point& p, const std::vector<Point>& vec, double xs = 1.0, double ys = 1.0,
0013                     double tol = 1e-6) {
0014   for (auto& pt : vec) {
0015     if ((std::abs(pt.x() - p.x()) / xs < tol) && std::abs(pt.y() - p.y()) / ys < tol) {
0016       return true;
0017     }
0018   }
0019   return false;
0020 }
0021
0022 // check if a square in a ring
0023 inline bool rec_in_ring(const Point& pt, double sx, double sy, double rmin, double rmax,
0024                         double phmin, double phmax) {
0025   if (pt.r() > rmax || pt.r() < rmin) {
0026     return false;
0027   }
0028
0029   // check four corners
0030   std::vector<Point> pts{
0031       Point(pt.x() - sx / 2., pt.y() - sy / 2.),
0032       Point(pt.x() - sx / 2., pt.y() + sy / 2.),
0033       Point(pt.x() + sx / 2., pt.y() - sy / 2.),
0034       Point(pt.x() + sx / 2., pt.y() + sy / 2.),
0035   };
0036   for (auto& p : pts) {
0037     if (p.r() > rmax || p.r() < rmin || p.phi() > phmax || p.phi() < phmin) {
0038       return false;
0039     }
0040   }
0041   return true;
0042 }
0043
0044 // a helper function to recursively fill square in a ring
0045 void add_rectangle(Point p, std::vector<Point>& res, double sx, double sy, double rmin, double rmax,
0046                    double phmin, double phmax, int max_depth = 20, int depth = 0) {
0047   // std::cout << depth << "/" << max_depth << std::endl;
0048   // exceeds the maximum depth in searching or already placed
0049   if ((depth > max_depth) || (already_placed(p, res, sx, sy))) {
0050     return;
0051   }
0052
0053   bool in_ring = rec_in_ring(p, sx, sy, rmin, rmax, phmin, phmax);
0054   if (in_ring) {
0055     res.emplace_back(p);
0056   }
0057
0058   // continue search for a good placement or if no placement found yet
0059   if (in_ring || res.empty()) {
0060     // check adjacent squares
0061     add_rectangle(Point(p.x() + sx, p.y()), res, sx, sy, rmin, rmax, phmin, phmax, max_depth,
0062                   depth + 1);
0063     add_rectangle(Point(p.x() - sx, p.y()), res, sx, sy, rmin, rmax, phmin, phmax, max_depth,
0064                   depth + 1);
0065     add_rectangle(Point(p.x(), p.y() + sy), res, sx, sy, rmin, rmax, phmin, phmax, max_depth,
0066                   depth + 1);
0067     add_rectangle(Point(p.x(), p.y() - sy), res, sx, sy, rmin, rmax, phmin, phmax, max_depth,
0068                   depth + 1);
0069   }
0070 }
0071
0072 // fill squares
0073 std::vector<Point> fillRectangles(Point ref, double sx, double sy, double rmin, double rmax,
0074                                   double phmin, double phmax) {
0075   // convert (0, 2pi) to (-pi, pi)
0076   if (phmax > M_PI) {
0077     phmin -= M_PI;
0078     phmax -= M_PI;
0079   }
0080   // start with a seed square and find one in the ring
0081   // move to center
0082   ref = ref - Point(int(ref.x() / sx) * sx, int(ref.y() / sy) * sy);
0083
0084   std::vector<Point> res;
0085   add_rectangle(ref, res, sx, sy, rmin, rmax, phmin, phmax,
0086                 (int(rmax / sx) + 1) * (int(rmax / sy) + 1) * 2);
0087   return res;
0088 }
0089
0090 // check if a regular polygon is inside a ring
0091 bool poly_in_ring(const Point& p, int nsides, double lside, double rmin, double rmax, double phmin,
0092                   double phmax) {
0093   // outer radius is contained
0094   if ((p.r() + lside <= rmax) && (p.r() - lside >= rmin)) {
0095     return true;
0096   }
0097
0098   // inner radius is not contained
0099   double rin = std::cos(M_PI / nsides) * lside;
0100   if ((p.r() + rin > rmax) || (p.r() - rin < rmin)) {
0101     return false;
0102   }
0103
0104   // in between, check every corner
0105   for (int i = 0; i < nsides; ++i) {
0106     double phi = (i + 0.5) * 2. * M_PI / static_cast<double>(nsides);
0107     Point p2(p.x() + 2. * lside * std::sin(phi), p.y() + 2. * lside * std::cos(phi));
0108     if ((p2.r() > rmax) || (p2.r() < rmin) || p.phi() > phmax || p.phi() < phmin) {
0109       return false;
0110     }
0111   }
0112   return true;
0113 }
0114
0115 // recursively fill square (nside=4) or hexagon (nside=6) in a ring, other polygons won't work
0116 void add_poly(Point p, std::vector<Point>& res, int nsides, double lside, double rmin, double rmax,
0117               double phmin, double phmax, int max_depth = 20, int depth = 0) {
0118   // std::cout << depth << "/" << max_depth << std::endl;
0119   // exceeds the maximum depth in searching or already placed
0120   if ((depth > max_depth) || (already_placed(p, res, lside, lside))) {
0121     return;
0122   }
0123
0124   bool in_ring = poly_in_ring(p, nsides, lside, rmin, rmax, phmin, phmax);
0125   if (in_ring) {
0126     res.emplace_back(p);
0127   }
0128
0129   // recursively add neigbors, continue if it was a good placement or no placement found yet
0130   if (in_ring || res.empty()) {
0131     for (int i = 0; i < nsides; ++i) {
0132       double phi = i * 2. * M_PI / static_cast<double>(nsides);
0133       add_poly(Point(p.x() + 2. * lside * std::sin(phi), p.y() + 2. * lside * std::cos(phi)), res,
0134                nsides, lside, rmin, rmax, phmin, phmax, max_depth, depth + 1);
0135     }
0136   }
0137 }
0138
0139 std::vector<Point> fillHexagons(Point ref, double lside, double rmin, double rmax, double phmin,
0140                                 double phmax) {
0141   // convert (0, 2pi) to (-pi, pi)
0142   if (phmax > M_PI) {
0143     phmin -= M_PI;
0144     phmax -= M_PI;
0145   }
0146   // start with a seed and find one in the ring
0147   // move to center
0148   ref = ref - Point(int(ref.x() / lside) * lside, int(ref.y() / lside) * lside);
0149
0150   std::vector<Point> res;
0151   add_poly(ref, res, 6, lside, rmin, rmax, phmin, phmax, std::pow(int(rmax / lside) + 1, 2) * 2);
0152   return res;
0153 }
0154
0155 bool isPointInsidePolygon(Point p, std::vector<Point> vertices) {
0156   int n      = vertices.size();
0157   bool check = false; // check == false (outside the polygon), check == true (inside the polygon)
0158   const double tolerance = 0.000001;
0159
0160   // When the point overlaps with vertex in the tolerance.
0161   //
0162   for (int i = 0; i < n; i++)
0163     if (std::abs(p.x() - vertices[i].x()) < tolerance &&
0164         std::abs(p.y() - vertices[i].y()) < tolerance)
0165       check = !check;
0166
0167   // When the point is on the line connected two vertices in the tolerance.
0168   //
0169   if (check == false) {
0170     for (int i = 0, j = n - 1; i < n; j = i++)
0171       if (std::abs(p.x() - vertices[i].x()) < tolerance &&
0172           std::abs(p.x() - vertices[j].x()) < tolerance)
0173         if ((vertices[i].y() > p.y()) != (vertices[j].y() > p.y()))
0174           check = !check;
0175   }
0176   if (check == false) {
0177     for (int i = 0, j = n - 1; i < n; j = i++)
0178       if (std::abs(p.y() - vertices[i].y()) < tolerance &&
0179           std::abs(p.y() - vertices[j].y()) < tolerance)
0180         if ((vertices[i].x() > p.x()) != (vertices[j].x() > p.x()))
0181           check = !check;
0182   }
0183
0184   if (check == false) {
0185     for (int i = 0, j = n - 1; i < n; j = i++) {
0186       double ver_i    = vertices[i].y();
0187       double ver_j    = vertices[j].y();
0188       double criteria = (vertices[j].x() - vertices[i].x()) * (p.y() - vertices[i].y()) /
0189                             (vertices[j].y() - vertices[i].y()) +
0190                         vertices[i].x();
0191
0192       if (((ver_i > p.y()) != (ver_j > p.y())) &&
0193           (p.x() < criteria || std::abs(p.x() - criteria) < tolerance))
0194         check = !check;
0195     }
0196   }
0197
0198   return check;
0199 }
0200
0201 bool isBoxTotalInsidePolygon(Point box[4], std::vector<Point> vertices) {
0202   bool pt_check = true;
0203   for (int i = 0; i < 4; i++)
0204     pt_check = pt_check && isPointInsidePolygon(box[i], vertices);
0205   return pt_check;
0206 }
0207
0208 bool isBoxPartialInsidePolygon(Point box[4], std::vector<Point> vertices) {
0209   bool pt_check = false;
0210   for (int i = 0; i < 4; i++)
0211     pt_check = pt_check || isPointInsidePolygon(box[i], vertices);
0212   return pt_check;
0213 }
0214
0215 std::vector<std::pair<double, double>>
0216 getPolygonVertices(std::pair<double, double> center, double radius, double angle_0, int numSides) {
0217   std::vector<std::pair<double, double>> vertices;
0218   double angle = 2 * M_PI / numSides; // calculate the angle between adjacent vertices
0219   for (int i = 0; i < numSides; i++) {
0220     double x = center.first + radius * cos(i * angle + angle_0);
0221     double y = center.second + radius * sin(i * angle + angle_0);
0222     vertices.emplace_back(x, y); // add the vertex to the vector
0223   }
0224   return vertices;
0225 }
0226
0227 } // namespace epic::geo