椭球方程包含许多点

Question

我有大量的像素颜色（96种不同的颜色）：

我希望得到某种数学定义的概率区域，如this question：

我现在看到的主要障碍 - Google上的所有方法主要是关于可视化和二维空间，但是没有找到方程系数的算法，如：

A 1 X + B 1 Y + X 1 Y +反射器+ Paxs +尺寸+反转+ Bai + Saiz = 0

this paper对我来说太难以在python中实现它。 :(

无论如何，我只想确定一些像素是否或多或少存在于我所拥有的辩证中。

我尝试使用scikit群集来制作它，但由于只有一组数据，我失败了。并且创建表示每个像素颜色的数组2563元素似乎是错误的方式。

我想知道是否有一种简单的方法来确定这个点群的边界？或者，也许我只是过度思考它，还有类似OpenCV cv2.inRange（）函数的东西？

Answer 1

这可以通过椭球多项式的优化和拟合来解决。但是，我会从更快的几何方法开始：

1. 找到平均点位置 这将是你的椭圆体的中心

   p0 = sum (p[i]) / n      // average
   i = { 0,1,2,3,...,n-1 }  // of all points

   如果您的点密度不均匀，则使用边界框中心更安全。所以找到xmin,ymin,zmin,xmax,ymax,zmax，他们之间的中间是你的中心。
2. 找到最远点到中心 这将给你主要的半轴

   pa = p[j];
   |p[j]-p0| >= |p[i]-p0|   // max
   i = { 0,1,2,3,...,n-1 }  // of all points

3. 找到第二个半轴 所以矢量pa-p0是平面的法线，其他半轴应该是平面。所以从那架飞机找到最远点p0：

   pb = p[j];  
   |p[j]-p0| >= |p[i]-p0|   // max
   dot(pa-p0,p[j]-p0) == 0  // but inly if inside plane
   i = { 0,1,2,3,...,n-1 }  // from all points

   请注意，dot产品的结果可能不会精确为零，因此最好对此类测试进行测试： |dot(pa-p0,p[j]-p0)| <= 1e-3 您可以使用任何您想要的阈值（应该基于椭圆体大小）。
4. 找到最后的半轴 所以我们知道最后的半轴应该垂直于两者 (pa-p0) AND (pb-p0) 所以找到这样的观点：

   pc = p[j];  
   |p[j]-p0| >= |p[i]-p0|   // max
   dot(pa-p0,p[j]-p0) == 0  // but inly if inside plane
   dot(pb-p0,p[j]-p0) == 0  // and perpendicular also to b semi-axis
   i = { 0,1,2,3,...,n-1 }  // from all points

5. 椭球 现在您拥有了形成椭球所需的所有参数。矢量 (pa-p0),(pb-p0),(pc-p0) 是椭圆体的基础向量（可以使用叉积使它们垂直）。它们的大小为您提供半径。而p0是中心。您还可以使用此参数方程：

   a=pa-p0;
   b=pb-p0;
   c=pc-p0;
   p(u,v) = p0 + a*cos(u)*cos(v)
               + b*cos(u)*sin(v)
               + c*sin(u);
   u = < -0.5*PI , +0.5*PI >
   v = < 0.0 , 2.0*PI >

整个过程只是O(n)，结果可以作为优化和拟合的起点，以加快它们的速度，而不会降低精度。如果您想进一步提高准确度请参阅：

• How approximation search works

子链接显示了适合的例子......

你也可以看看这个：

• Algorithms: Ellipse matching

这基本上类似于你的任务，但只在2D中仍然可以带给你一些想法。

Answer 2

有一种方法可以找到基于两个焦点的非完美对称椭圆体。可能它对于大量积分是有益的。至少，它非常简单（基于随机搜索）：

3DPoints - Array of some Amount of points
vecCenter - average of 3DPoints
AngleCosine - cos of angle between two vectors

RandomOrder(3DPoints)

vecFocus := (0, 0, 0)
for i := 0 to Amount:
    vecRadius := 3DPoints[i] - vecCenter

    // Change vecRadius direction to parallel if it's not
    if AngleCosine(vecFocus, vecRadius) < 0 then
        vecRadius *= -1
    vecFocus += (vecRadius - vecFocus) / Amount

通过向后传递第二次数组可以改善结果然后将焦点置于vecCenter +/- vecFocus。然后到焦点的距离可以通过某个系数提出或者用第二次通过来计算，如下所示：

FinalRadius := 0
for i := 0 to Amount:
    vecRadius := (3DPoints[i] - vecFoci1) + (3DPoints[i] - vecFoci2)
    FinalRadius := Max(FinalRadius, Length(vecRadius))