...
 
Commits (3)
......@@ -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};
}
......@@ -16,7 +16,7 @@
\kactlimport{CircleIntersection.h}
\kactlimport{CircleTangents.h}
% \kactlimport{CircleLine.h}
% \kactlimport{CirclePolygonIntersection.h}
\kactlimport{CirclePolygonIntersection.h}
\kactlimport{circumcircle.h}
\kactlimport{MinimumEnclosingCircle.h}
......
......@@ -12,7 +12,6 @@
\section{Dynamic programming}
\kactlimport{KnuthDP.h}
\kactlimport{DivideAndConquerDP.h}
\newpage % the rest fits nicely on one page
\section{Debugging tricks}
\begin{itemize}
......@@ -23,6 +22,7 @@
\end{itemize}
\section{Optimization tricks}
\verb@__builtin_ia32_ldmxcsr(40896);@ disables denormals (which make floats 20x slower near their minimum value).
\subsection{Bit hacks}
\begin{itemize}
\item \verb@x & -x@ is the least bit in \texttt{x}.
......
#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;
}