How can I check wether a point is inside the circumcircle of 3 points?

Question

Is there any easy solution? Or has anybody an example of an implementation?

Thanks, Jonas

Answer 1

The mathematical intuition behind TitouanT could be explained by the following:

1. The determinant calculates the signed volume of the parallelotop formed by the vectors.

2. (ax, ay, ax² + ay²) is the result of projecting the vector (ax, ay) onto the paraboloid f(v) = (vx, vy, ||v||²).

3. When you take the plane of the triangle in the paraboloid, every point within the circumcircle is in front of the plane, and every point outside the triangle is behind the plane.

Consider the matrix:

ax, ay, ax^2 + ay^2
bx, by, bx^2 + by^2
cx, cy, cx^2 + cy^2

This matrix can be interpreted as:

f(a)
f(b)
f(c)

and:

ax - px, ay - py, (ax - px)^2 + (ay - py)^2
bx - px, by - py, (bx - px)^2 + (by - py)^2
cx - px, cy - py, (cx - px)^2 + (cy - py)^2

as:

f(a) - f(p)
f(b) - f(p)
f(c) - f(p)

So what the matrix does is it project the points into the paraboloid and takes the vectors from p to the triangle vertices.

Remember, the determinant calculates the signed volume. If p is in front of the plane, the volume is positive. If p is behind the plane, the volume is negative.

To visualize this, here is a plot of the paraboloid, the triangle and the point p and it's projections. First picture the point is within the circumcircle and in the second one it's not.

Plot 1

Plot 2

Answer 2

(In case you are interested in a non-obvious/"crazy" kind of solution.)

One equivalent property of Delaunay triangulation is as follows: if you build a circumcircle of any triangle in the triangulation, it is guaranteed not to contain any other vertices of the triangulation.

Another equivalent property of Delaunay triangulation is: it maximizes the minimal triangle angle (i.e. maximizes it among all triangulations on the same set of points).

This suggests an algorithm for your test:

1. Consider triangle ABC built on the original 3 points.
2. If the test point P lies inside the triangle it is definitely inside the circle
3. If the test point P belongs to one of the "corner" regions (see the shaded regions in the picture below), it is definitely outside the circle

1. Otherwise (let's say P lies in region 1) consider two triangulations of quadrilateral ABCP: the original one contains the original triangle (red diagonal), and the alternate one with "flipped" diagonal (blue diagonal)

1. Determine which one if this triangulations is a Delaunay triangulation by testing the "flip" condition, e.g. by comparing α = min(∠1,∠4,∠5,∠8) vs. β = min(∠2,∠3,∠6,∠7).
2. If the original triangulation is a Delaunay triangulation (α > β), P lies outside the circle. If the alternate triangulation is a Delaunay triangulation (α < β), P lies inside the circle.

Done.

---

(Revisiting this answer after a while.)

This solution might not be as "crazy" as it might appear at the first sight. Note that in order to compare angles at steps 5 and 6 there's no need to calculate the actual angle values. It is sufficient to know their cosines (i.e. there's no need to involve trigonometric functions).

A C++ version of the code

#include <cmath>
#include <array>
#include <algorithm>

struct pnt_t
{
  int x, y;

  pnt_t ccw90() const
    { return { -y, x }; }

  double length() const
    { return std::hypot(x, y); }

  pnt_t &operator -=(const pnt_t &rhs)
  {
    x -= rhs.x;
    y -= rhs.y;
    return *this;
  }

  friend pnt_t operator -(const pnt_t &lhs, const pnt_t &rhs)
    { return pnt_t(lhs) -= rhs; }

  friend int operator *(const pnt_t &lhs, const pnt_t &rhs)
    { return lhs.x * rhs.x + lhs.y * rhs.y; }
};

int side(const pnt_t &a, const pnt_t &b, const pnt_t &p)
{
  int cp = (b - a).ccw90() * (p - a);
  return (cp > 0) - (cp < 0);
}

void make_ccw(std::array<pnt_t, 3> &t)
{
  if (side(t[0], t[1], t[2]) < 0)
    std::swap(t[0], t[1]);
}

double ncos(pnt_t a, const pnt_t &o, pnt_t b)
{
  a -= o;
  b -= o;
  return -(a * b) / (a.length() * b.length());
}

bool inside_circle(std::array<pnt_t, 3> t, const pnt_t &p)
{
  make_ccw(t);

  std::array<int, 3> s = 
    { side(t[0], t[1], p), side(t[1], t[2], p), side(t[2], t[0], p) };

  unsigned outside = std::count(std::begin(s), std::end(s), -1);
  if (outside != 1)
    return outside == 0;

  while (s[0] >= 0)
  {
    std::rotate(std::begin(t), std::begin(t) + 1, std::end(t));
    std::rotate(std::begin(s), std::begin(s) + 1, std::end(s));
  }

  double 
    min_org = std::min({
      ncos(t[0], t[1], t[2]), ncos(t[2], t[0], t[1]), 
      ncos(t[1], t[0], p), ncos(p, t[1], t[0]) }),
    min_alt = std::min({
      ncos(t[1], t[2], p), ncos(p, t[2], t[0]), 
      ncos(t[0], p, t[2]), ncos(t[2], p, t[1]) });

  return min_org <= min_alt;
}

and a couple of tests with arbitrarily chosen triangles and a large number of random points

---

Of course, the whole thing can be easily reformulated without even mentioning Delaunay triangulations. Starting from step 4 this solution is based in the property of the opposite angles of cyclic quadrilateral, which must sum to 180°.

Answer 3

Lets call

• a, b, c our three points,
• C the circumcircle of (a, b, c)
• and d an other point.

A fast way to determine if d is in C is to compute this determinant:

      | ax-dx, ay-dy, (ax-dx)² + (ay-dy)² |
det = | bx-dx, by-dy, (bx-dx)² + (by-dy)² |
      | cx-dx, cy-dy, (cx-dx)² + (cy-dy)² |

if a, b, c are in counter clockwise order then:

• if det equal 0 then d is on C
• if det > 0 then d is inside C
• if det < 0 then d is outside C

here is a javascript function that does just that:

function inCircle (ax, ay, bx, by, cx, cy, dx, dy) {
    let ax_ = ax-dx;
    let ay_ = ay-dy;
    let bx_ = bx-dx;
    let by_ = by-dy;
    let cx_ = cx-dx;
    let cy_ = cy-dy;
    return (
        (ax_*ax_ + ay_*ay_) * (bx_*cy_-cx_*by_) -
        (bx_*bx_ + by_*by_) * (ax_*cy_-cx_*ay_) +
        (cx_*cx_ + cy_*cy_) * (ax_*by_-bx_*ay_)
    ) > 0;
}

You might also need to check if your points are in counter clockwise order:

function ccw (ax, ay, bx, by, cx, cy) {
    return (bx - ax)*(cy - ay)-(cx - ax)*(by - ay) > 0;
}

I didn't place the ccw check inside the inCircle function because you shouldn't check it every time.

This process doesn't take any divisions or square root operation.

You can see the code in action there: https://titouant.github.io/testTriangle/

And the source there: https://github.com/TitouanT/testTriangle

Answer 4

Given 3 points (x1,y1),(x2,y2),(x3,y3) and the point you want to check is inside the circle defined by the above 3 points (x,y) you can do something like

/**
 *
 * @param x coordinate of point want to check if inside
 * @param y coordinate of point want to check if inside
 * @param cx center x
 * @param cy center y
 * @param r radius of circle
 * @return whether (x,y) is inside circle
 */
static boolean g(double x,double y,double cx,double cy,double r){
    return Math.sqrt((x-cx)*(x-cx)+(y-cy)*(y-cy))<r;
}

// check if (x,y) is inside circle defined by (x1,y1),(x2,y2),(x3,y3)
static boolean isInside(double x,double y,double x1,double y1,double x2,double y2,double x3,double y3){
    double m1 = (x1-x2)/(y2-y1);
    double m2 = (x1-x3)/(y3-y1);
    double b1 = ((y1+y2)/2) - m1*(x1+x2)/2;
    double b2 = ((y1+y3)/2) - m2*(x1+x3)/2;
    double xx = (b2-b1)/(m1-m2);
    double yy = m1*xx + b1;
    return g(x,y,xx,yy,Math.sqrt((xx-x1)*(xx-x1)+(yy-y1)*(yy-y1)));
}

public static void main(String[] args) {
    // if (0,1) is inside the circle defined by (0,0),(0,2),(1,1)
    System.out.println(isInside(0,1,0,0,0,2,1,1));
}

The method for getting an expression for the center of circle from 3 points goes from finding the intersection of the 2 perpendicular bisectors of 2 line segments, above I chose (x1,y1)-(x2,y2) and (x1,y1)-(x3,y3). Since you know a point on each perpendicular bisector, namely (x1+x2)/2 and (x1+x3)/2, and since you also know the slope of each perpendicular bisector, namely (x1-x2)/(y2-y1) and (x1-x3)/(y3-y1) from the above 2 line segments respectively, you can solve for the (x,y) where they intersect.

Answer 5

In this Math SE post of mine I included an equation which checks if four points are cocircular by computing a 4×4 determinant. By turning that equation into an inequality you can check for insideness.

If you want to know which direction the inequality has to go, conisder the case of a point very far away. In this case, the x²+y² term will dominate all other terms. So you can simply assume that for the point in question, this term is one while the three others are zero. Then pick the sign of your inequality so this value does not satisfy it, since this point is definitely outside but you want to characterize inside.

If numeric precision is an issue, this page by Prof. Shewchuk describes how to obtain consistent predicates for points expressed using regular double precision floating point numbers.