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

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 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°.