Split 2022 vision math code into library

This allows it to be used in target estimator and camera reader in a
addition to blob detector. Will be used for comparing size of projected
blobs to detected blobs, as well as for interpolating between localizer
outputs.

Signed-off-by: milind-u <milind.upadhyay@gmail.com>
Change-Id: I80a726e93d638106431614459837e8671bcb7251
diff --git a/y2022/vision/blob_detector.cc b/y2022/vision/blob_detector.cc
index 2c68a27..68832fd 100644
--- a/y2022/vision/blob_detector.cc
+++ b/y2022/vision/blob_detector.cc
@@ -8,6 +8,7 @@
 #include "aos/time/time.h"
 #include "opencv2/features2d.hpp"
 #include "opencv2/imgproc.hpp"
+#include "y2022/vision/geometry.h"
 
 DEFINE_uint64(red_delta, 100,
               "Required difference between green pixels vs. red");
@@ -84,105 +85,6 @@
   return blob_stats;
 }
 
-namespace {
-
-// Linear equation in the form ax + by = c
-struct Line {
- public:
-  double a, b, c;
-
-  std::optional<cv::Point2d> Intersection(const Line &l) const {
-    // Use Cramer's rule to solve for the intersection
-    const double denominator = Determinant(a, b, l.a, l.b);
-    const double numerator_x = Determinant(c, b, l.c, l.b);
-    const double numerator_y = Determinant(a, c, l.a, l.c);
-
-    std::optional<cv::Point2d> intersection = std::nullopt;
-    // Return nullopt if the denominator is 0, meaning the same slopes
-    if (denominator != 0) {
-      intersection =
-          cv::Point2d(numerator_x / denominator, numerator_y / denominator);
-    }
-
-    return intersection;
-  }
-
- private:  // Determinant of [[a, b], [c, d]]
-  static double Determinant(double a, double b, double c, double d) {
-    return (a * d) - (b * c);
-  }
-};
-
-struct Circle {
- public:
-  cv::Point2d center;
-  double radius;
-
-  static std::optional<Circle> Fit(std::vector<cv::Point2d> centroids) {
-    CHECK_EQ(centroids.size(), 3ul);
-    // For the 3 points, we have 3 equations in the form
-    // (x - h)^2 + (y - k)^2 = r^2
-    // Manipulate them to solve for the center and radius
-    // (x1 - h)^2 + (y1 - k)^2 = r^2 ->
-    // x1^2 + h^2 - 2x1h + y1^2 + k^2 - 2y1k = r^2
-    // Also, (x2 - h)^2 + (y2 - k)^2 = r^2
-    // Subtracting these two, we get
-    // x1^2 - x2^2 - 2h(x1 - x2) + y1^2 - y2^2 - 2k(y1 - y2) = 0 ->
-    // h(x1 - x2) + k(y1 - y2) = (-x1^2 + x2^2 - y1^2 + y2^2) / -2
-    // Doing the same with equations 1 and 3, we get the second linear equation
-    // h(x1 - x3) + k(y1 - y3) = (-x1^2 + x3^2 - y1^2 + y3^2) / -2
-    // Now, we can solve for their intersection and find the center
-    const auto l =
-        Line{centroids[0].x - centroids[1].x, centroids[0].y - centroids[1].y,
-             (-std::pow(centroids[0].x, 2) + std::pow(centroids[1].x, 2) -
-              std::pow(centroids[0].y, 2) + std::pow(centroids[1].y, 2)) /
-                 -2.0};
-    const auto m =
-        Line{centroids[0].x - centroids[2].x, centroids[0].y - centroids[2].y,
-             (-std::pow(centroids[0].x, 2) + std::pow(centroids[2].x, 2) -
-              std::pow(centroids[0].y, 2) + std::pow(centroids[2].y, 2)) /
-                 -2.0};
-    const auto center = l.Intersection(m);
-
-    std::optional<Circle> circle = std::nullopt;
-    if (center) {
-      // Now find the radius
-      const double radius = cv::norm(centroids[0] - *center);
-      circle = Circle{*center, radius};
-    }
-    return circle;
-  }
-
-  double DistanceTo(cv::Point2d p) const {
-    const auto p_prime = TranslateToOrigin(p);
-    // Now, the distance is simply the difference between distance from the
-    // origin to p' and the radius.
-    return std::abs(cv::norm(p_prime) - radius);
-  }
-
-  double AngleOf(cv::Point2d p) const {
-    auto p_prime = TranslateToOrigin(p);
-    // Flip the y because y values go downwards.
-    p_prime.y *= -1;
-    return std::atan2(p_prime.y, p_prime.x);
-  }
-
-  // TODO(milind): handle wrapping around 2pi
-  bool InAngleRange(cv::Point2d p, double theta_min, double theta_max) const {
-    const double theta = AngleOf(p);
-    return (theta >= theta_min && theta <= theta_max);
-  }
-
- private:
-  // Translate the point on the circle
-  // as if the circle's center is the origin (0,0)
-  cv::Point2d TranslateToOrigin(cv::Point2d p) const {
-    return cv::Point2d(p.x - center.x, p.y - center.y);
-  }
-};
-
-}  // namespace
-
 void BlobDetector::FilterBlobs(BlobResult *blob_result) {
   std::vector<std::vector<cv::Point>> filtered_blobs;
   std::vector<BlobStats> filtered_stats;