计算几何：In Circle Test

1. 背景介绍

最近在看清华大学x的计算几何，《Delaunay Triangulation》章节，里面有提到一个非常重要的判别方法：

已知A，B和C三点，和它们形成的△ABC会有一个外接圆O，判断给定P点是否在圆O的内部

如下面教学PPT¹所示：

那这篇文章就讨论一下InCircle()函数的实现方法，有需要的童鞋可以直接取用呢~

特别感谢：@我们家的橙酱 （GitHub）对本文的帮助哒

2. 方法解析

其实PPT里面已经给到了详细的计算方法，就是计算A、B和C组成的行列式：

$$\left[\begin{array}{cccc}{a}_x& a_y& a_x^2+a_y^2& 1\\ b_x& b_y& b_x^2+b_y^2& 1\\ c_x& c_y& c_x^2+c_y^2& 1\\ p_x& p_y& p_x^2+p_y^2& 1\end{array}\right]$$

但是四阶行列式的计算非常麻烦，而且博主也没有学过线性代数，所以找找了这个计算等式的结果，因为最后一列全是1，所以应该是由简化的方法。因此，根据这篇文章《Adaptive Precision Floating-Point Arithmetic》的结论²，上面的四阶行列式可以转换成下面的三阶行列式：

$$\left[\begin{array}{ccc}{a}_x-p_x& a_y-p_y& (a_x-p_x{)}^2+(a_y-p_y{)}^2\\ b_x-p_x& b_y-p_y& (b_x-p_x{)}^2+(b_y-p_y{)}^2\\ c_x-p_x& c_y-p_y& (c_x-p_x{)}^2+(c_y-p_y{)}^2\end{array}\right]$$
对于三阶行列式，我们就可以直接使用对角线法则进行直接计算：

$$\begin{gathered} \left[\begin{array}{ccc}a& b& c\\ d& e& f\\ g& h& i\end{array}\right] \\ =a\ast e\ast i+b\ast f\ast g+d\ast h\ast c-c\ast e\ast g-b\ast d\ast i-f\ast h\ast a \\ =a\ast (e\ast i-f\ast h)+b\ast (f\ast g-d\ast i)+c\ast (d\ast h-e\ast g) \end{gathered}$$

对于行列式的结果R，可以有以下结论：

如果R > 0，则P点在圆O的内部；

如果R == 0，则P点在圆O上，即圆周上；

如果R < 0，则P点在圆O的外部；

这里拓展一下，这个In Circle Test不仅可以用在计算几何里面，如果有面试题问已知三点，问如何判断另一点在不在这三点形成三角形的外接圆内部，这类题目可以直接用这里的方法破解哦。

3. 代码解析

最后简单过一遍给到的代码哦，首先是InCircle()：

/**
 * determine if p lies in the circumcircle of triangle(a, b, c)
 *
 * Note that point a, b and c must be in counter-clock wise order
 *
 * Reference resource:
 * @see <a href=https://www.cs.cmu.edu/~quake/robust.html>Adaptive Precision Floating-Point</a>
 * @see <a href=https://www.cs.umd.edu/class/spring2020/cmsc754/Lects/lect13-delaun-alg.pdf>Delaunay Triangulations: Incremental Construction</a>
 *
 * */

public static
double inCircle( Vector a, Vector b, Vector c, Vector p ) {
    if ( MyMath.isEqualZero( Triangles.areaTwo( a, b, c ) ) )
        throw new ArithmeticException( "not a triangle in inCircleTest()" );

    final double DIFF_X_A_P = a.x - p.x;
    final double DIFF_Y_A_P = a.y - p.y;
    final double DIFF_X_B_P = b.x - p.x;
    final double DIFF_Y_B_P = b.y - p.y;
    final double DIFF_X_C_P = c.x - p.x;
    final double DIFF_Y_C_P = c.y - p.y;
    double[][] matrix = {
        { DIFF_X_A_P, DIFF_Y_A_P, DIFF_X_A_P * DIFF_X_A_P + DIFF_Y_A_P * DIFF_Y_A_P },
        { DIFF_X_B_P, DIFF_Y_B_P, DIFF_X_B_P * DIFF_X_B_P + DIFF_Y_B_P * DIFF_Y_B_P },
        { DIFF_X_C_P, DIFF_Y_C_P, DIFF_X_C_P * DIFF_X_C_P + DIFF_Y_C_P * DIFF_Y_C_P }
    };

    return new Determinant( matrix ).getResult();
}

代码基本就是上面思路的直接实现，但需要注意两点：

1）输入点a, b, c必须为逆时针顺序

根据教材³上面的假设：

Second, by relabeling a, b, and c, we may assume that the counterclockwise order of points is <abcd>

构成内接圆三角形的输入三点必须是逆时针顺序，否则行列式无法计算出正确的值。最新的代码已经使用Duliaty，使得任意三点输入都不受上述限制，具体实现思路大家可以参考Code。

2）三点共线的情况

if ( MyMath.isEqualZero( Triangles.areaTwo( a, b, c ) ) )
	throw new ArithmeticException( "not a triangle in inCircleTest()" );

这里是排除了A、B、C三点共线的情况，用的toLeft测试，算是一个特殊情况。至于三阶行列式的计算代码，大家可以在Determinant类里面找到，这里就不再赘述。

4. 拓展阅读

1. 计算几何课堂：几何寻路之旅
2. 向量复习（一）：定义、求解、四则运算、点积和叉积
3. 超详细！红黑树详解文章汇总（含代码）
4. 几何求交（一）：直线和直线的交点
5. 几何求交（二）：直线和圆的交点

5. 附录：项目代码

Description	Entry method
toLeft test	boolean toLeft( Vector base1, Vector base2, Vector point )
inCircle test	double inCircle( Vector a, Vector b, Vector c, Vector p )

6. 参考资料

1. 计算几何 | Computational Geometry
2. Adaptive Precision Floating-Point Arithmetic and Fast Robust Predicates for Computational Geometry
3. Delaunay Triangulations: Incremental Construction

7. 免责声明

※ 本文之中如有错误和不准确的地方，欢迎大家指正哒~
※ 此项目仅用于学习交流，请不要用于任何形式的商用用途，谢谢呢；

---

---

1. 计算几何 | Computational Geometry ↩︎

2. Adaptive Precision Floating-Point Arithmetic and Fast Robust Predicates for Computational Geometry ↩︎

3. Delaunay Triangulations: Incremental Construction ↩︎