PixelSeeding/test/CircleEq_t.cpp

0001 #include "RecoTracker/PixelSeeding/interface/CircleEq.h"
0002 #include <cassert>
0003
0004 struct OriCircle {
0005   using T = float;
0006
0007   float radius = 0;
0008   float x_center = 0;
0009   float y_center = 0;
0010
0011   constexpr OriCircle(T x1, T y1, T x2, T y2, T x3, T y3) { compute(x1, y1, x2, y2, x3, y3); }
0012
0013   // dca to origin
0014   constexpr T dca0() const { return std::sqrt(x_center * x_center + y_center * y_center) - radius; }
0015
0016   // dca to given point
0017   constexpr T dca(T x, T y) const {
0018     x -= x_center;
0019     y -= y_center;
0020     return std::sqrt(x * x + y * y) - radius;
0021   }
0022
0023   constexpr void compute(T x1, T y1, T x2, T y2, T x3, T y3) {
0024     auto det = (x1 - x2) * (y2 - y3) - (x2 - x3) * (y1 - y2);
0025
0026     auto offset = x2 * x2 + y2 * y2;
0027
0028     auto bc = (x1 * x1 + y1 * y1 - offset) * 0.5f;
0029
0030     auto cd = (offset - x3 * x3 - y3 * y3) * 0.5f;
0031
0032     auto idet = 1.f / det;
0033
0034     x_center = (bc * (y2 - y3) - cd * (y1 - y2)) * idet;
0035     y_center = (cd * (x1 - x2) - bc * (x2 - x3)) * idet;
0036
0037     radius = std::sqrt((x2 - x_center) * (x2 - x_center) + (y2 - y_center) * (y2 - y_center));
0038   }
0039 };
0040
0041 #include <iostream>
0042
0043 template <typename T>
0044 bool equal(T a, T b) {
0045   //  return float(a-b)==0;
0046   return std::abs(float(a - b)) < std::abs(0.01f * a);
0047 }
0048
0049 int main() {
0050   float r1 = 4, r2 = 8, r3 = 15;
0051   for (float phi = -3; phi < 3.1; phi += 0.5) {
0052     float x1 = r1 * cos(phi);
0053     float x2 = r2 * cos(phi);
0054     float y1 = r1 * sin(phi);
0055     float y2 = r2 * sin(phi);
0056     for (float phi3 = phi - 0.31; phi3 < phi + 0.31; phi3 += 0.05) {
0057       float x3 = r3 * cos(phi3);
0058       float y3 = r3 * sin(phi3);
0059
0060       OriCircle ori(x1, y1, x2, y2, x3, y3);
0061       CircleEq<float> eq(x1, y1, x2, y2, x3, y3);
0062       // std::cout << "r " << ori.radius <<' '<< eq.radius() << std::endl;
0063       assert(equal(ori.radius, std::abs(eq.radius())));
0064       auto c = eq.center();
0065       auto dir = eq.cosdir();
0066       assert(equal(1.f, dir.first * dir.first + dir.second * dir.second));
0067       assert(equal(ori.x_center, c.first));
0068       assert(equal(ori.y_center, c.second));
0069       // std::cout << "dca " << ori.dca0() <<' '<< eq.radius()*eq.dca0() << std::endl;
0070       assert(equal(std::abs(ori.dca0()), std::abs(eq.radius() * eq.dca0())));
0071       // std::cout << "dca " << ori.dca(1.,1.) <<' '<< eq.radius()*eq.dca(1.,1.) << std::endl;
0072       assert(equal(std::abs(ori.dca(1., 1.)), std::abs(eq.radius() * eq.dca(1., 1.))));
0073     }
0074   }
0075
0076   return 0;
0077 }