Rewrite CircleLine to avoid LineProjectionReflection

Author: simonlindholm

content/geometry/CircleLine.h

@@ -3,23 +3,22 @@
  * Date: 2019-10-29
  * License: CC0
  * Source: https://vlecomte.github.io/cp-geo.pdf
- * Description: Finds the intersection between a circle and a line. Returns a
- * vector of either 0, 1, or 2 intersection points. P is intended to be
- * Point<double>
+ * Description: Finds the intersection between a circle and a line.
+ * Returns a vector of either 0, 1, or 2 intersection points.
+ * P is intended to be Point<double>.
  * Status: unit tested
  */
 
 #pragma once
 
 #include "Point.h"
-#include "lineDistance.h"
-#include "LineProjectionReflection.h"
 
 template<class P>
 vector<P> circleLine(P c, double r, P a, P b) {
-	double h2 = r*r - a.cross(b,c)*a.cross(b,c)/(b-a).dist2();
+	P ab = b - a, p = a + ab * (c-a).dot(ab) / ab.dist2();
+	double s = a.cross(b, c), h2 = r*r - s*s / ab.dist2();
 	if (h2 < 0) return {};
-	P p = lineProj(a, b, c), h = (b-a).unit() * sqrt(h2);
 	if (h2 == 0) return {p};
+	P h = ab.unit() * sqrt(h2);
 	return {p - h, p + h};
 }

stress-tests/geometry/CircleLine.cpp

@@ -1,12 +1,11 @@
 #include "../utilities/template.h"
+#include "../utilities/randGeo.h"
 
+#include "../../content/geometry/lineDistance.h"
 #include "../../content/geometry/CircleLine.h"
 
 typedef Point<double> P;
 int main() {
-    cin.sync_with_stdio(0);
-    cin.tie(0);
-    cin.exceptions(cin.failbit);
     {
         auto res = circleLine(P(0, 0), 1, P(-1, -1), P(1, 1));
         assert(res.size() == 2);
@@ -21,5 +20,31 @@ int main() {
         auto res = circleLine(P(4, 4), 1, P(0,  0), P(5, 0));
         assert(res.size() == 0);
     }
+	rep(it,0,100000) {
+		P a = randIntPt(5);
+		P b = randIntPt(5);
+		P c = randIntPt(5);
+		if (a == b) {
+			// Not a well defined line
+			continue;
+		}
+		double r = sqrt(rand() % 49);
+		vector<P> points = circleLine(c, r, a, b);
+
+		// Soundness
+		assert(sz(points) <= 2);
+		for (P p : points) {
+			// Point is on circle
+			assert(abs((p - c).dist() - r) < 1e-6);
+			// Point is on line
+			assert(lineDist(a, b, p) < 1e-6);
+		}
+
+		// Best-effort completeness check:
+		// in some easy cases we must have points in the intersection.
+		if ((a - c).dist() < r - 1e-6 || (b - c).dist() < r - 1e-6 || ((a + b) / 2 - c).dist() < r - 1e-6) {
+			assert(!points.empty());
+		}
+	}
     cout<<"Tests passed!"<<endl;
 }