你用极坐标分区的思路非常正确，但生成所有点组合再逐个判断的方式，本质还是暴力枚举——这对于N=12000来说，计算量是天文数字（哪怕分成三个部分，组合数也远超600亿），慢是必然的。我们可以利用质心在三角形内的几何判定规则，结合极角排序和双指针法，把复杂度从O(N³)直接降到O(N log N)，几秒就能出结果。

核心几何逻辑 #

质心 ( C_t ) 在三角形内部的充要条件是：三个点相对于 ( C_t ) 的极角，不能被一个180°的半平面完全覆盖。反过来想：所有不包含 ( C_t ) 的三角形，三个点的极角肯定都落在某个180°的范围内。所以我们只需要算出总三角形数，再减去“极角全在180°半平面内”的三角形数量，就能得到目标结果。

具体操作步骤 #

1. 转换为相对质心的极角 #

先把每个点的坐标转换为以 ( C_t ) 为原点的相对坐标，再计算极角（范围0到2π，避免负角度的麻烦）：

import numpy as np
import pandas as pd

# 假设你已经算出质心Ct的坐标(ct_x, ct_y)
PntsDF['rel_x'] = PntsDF['x'] - ct_x
PntsDF['rel_y'] = PntsDF['y'] - ct_y

# 计算极角，把负角度转换为0-2π范围内的值
PntsDF['theta'] = np.arctan2(PntsDF['rel_y'], PntsDF['rel_x'])
PntsDF['theta'] = np.where(PntsDF['theta'] < 0, PntsDF['theta'] + 2*np.pi, PntsDF['theta'])

2. 排序极角并扩展数组 #

把极角排序后，复制一份并加上2π——这样可以轻松处理跨0°（也就是2π）的半平面情况（比如从350°到10°的180°范围）：

theta_sorted = np.sort(PntsDF['theta'].values)
theta_extended = np.concatenate([theta_sorted, theta_sorted + 2*np.pi])
n = len(theta_sorted)

3. 双指针法统计不包含质心的三角形 #

对于每个点 ( i )（对应排序后的极角 ( \theta_i )），用双指针找到最远的点 ( j )，使得 ( \theta_j - \theta_i < 180°(\pi) )。这部分点的数量是 ( k = j - i - 1 )，从这 ( k ) 个点里选2个和 ( i ) 组成的三角形，肯定都不包含 ( C_t )：

count_exclude = 0
j = 0
for i in range(n):
    # 找到第一个不在θ_i的180°范围内的点
    while theta_extended[j] - theta_sorted[i] < np.pi:
        j += 1
    # 计算当前i对应的、不包含Ct的三角形数量
    k = j - i - 1
    if k >= 2:
        count_exclude += k * (k - 1) // 2  # 组合数C(k,2)

4. 计算最终结果 #

总三角形数是组合数 ( \binom{N}{3} )，减去不包含的数量就是包含 ( C_t ) 的三角形数：

total_triangles = n * (n - 1) * (n - 2) // 6
count_include = total_triangles - count_exclude

为什么这个方法快？ #

• 整个流程的时间瓶颈是排序（O(N log N)），对于N=12000来说，排序只需要几毫秒，后续的双指针遍历是O(N)的线性时间，全程用numpy数组运算，完全避开了pandas行级查找的低效操作。
• 不需要生成任何点组合，直接通过数学计算得到结果，彻底解决了组合数爆炸的问题。

小补充 #

如果你的数据量再大（比如N=1e5），这个方法依然适用——它的效率是线性对数级的，远优于任何暴力枚举的思路。

内容的提问来源于stack exchange，提问作者Valuex