给定任意平面上的噪声数据点优化圆心

Question

我似乎没有找到一种算法来计算 3D 中任意平面上给定一组噪声数据（在圆周上）的圆心。 我确实尝试了不同的方法来计算中心，一次只取 3 个点，然后对中心进行平均，但我相信必须有更好的解决方案来通过最小化到所有点的距离来优化中心。 有人有想法吗？

任何想法或方法将不胜感激！

Answer 1

最后我找到了几个解决方案调整我的搜索词：

https://de.mathworks.com/matlabcentral/answers/475212-circle-least-squares-fit-for-3d-data

https://jekel.me/2015/Least-Squares-Sphere-Fit/

两种方法都非常有效。在第二个链接中，您必须先计算平面，这可以使用 SVD 来完成。基本上我认为这两种方法非常相似并且产生基本相同的结果。

Answer 2

我是这样做的：

我有一个辅助结构（称为 Aprox），它包含中心、半径（我用它的正方形来优化速度）和边界点（两个或三个）。

我有以下功能：

CircumscribedCircle_2Points(P1, P2):Aprox（3D 中线段的中心）

CircumscribedCircle_3Points（P1，P2，P3）：Aprox（圆心追踪），由https://en.wikipedia.org/wiki/Circumscribed_circle#Cartesian_coordinates_from_cross-_and_dot-products

编程

CircumscribedCircle_IsIn($Aprox, $Pnt):boolean (It is in the given Approx)

我从前两点计算 Aprox，然后遍历其他点。如果它在圈内，我什么都不做，如果它在圈外，我计算一个新的 Aprox 并再次运行（所有点）。 规则是：如果平面上有 4 个点，其中一个点总是在其他三个点给出的圆圈内。

function CircumscribedCircle_2Points($P1, $P2,
  $Aprox=[])
{
$Aprox = as_New('Borders', [$P1, $P2]);
$Aprox['Center'] = v3d_AvgV3($P1, $P2);
$Aprox['RadSq'] = v3d_mV3AbsKv($Aprox['Center'], $P1);
$Aprox['RadSq'] = $Aprox['RadSq'] + 0.0001; //Float rounding errors
return $Aprox;
}

function CircumscribedCircle_3Points($P1, $P2, $P3,
  $Aprox=[],$JmP=0,$Jm=0,$u1=0,$u2=0,$u3=0)
//https://en.wikipedia.org/wiki/Circumscribed_circle#Cartesian_coordinates_from_cross-_and_dot-products
{
$Aprox = as_New('Borders', [$P1, $P2, $P3]);
$JmP = v3d_Abs(v3d_vxV3(v3d_mV3($P1,$P2),v3d_mV3($P2,$P3)));
$Jm = 2 * $JmP**2;
$u1 = v3d_mV3AbsKv($P2,$P3) * v3d_sxV3(v3d_mV3($P1,$P2),v3d_mV3($P1,$P3)) / $Jm;
$u2 = v3d_mV3AbsKv($P1,$P3) * v3d_sxV3(v3d_mV3($P2,$P1),v3d_mV3($P2,$P3)) / $Jm;
$u3 = v3d_mV3AbsKv($P1,$P2) * v3d_sxV3(v3d_mV3($P3,$P1),v3d_mV3($P3,$P2)) / $Jm;
$Aprox['Center'] = v3d_pV3pV3(v3d_xSkl($P1,$u1), v3d_xSkl($P2,$u2), v3d_xSkl($P3,$u3));
$Aprox['RadSq'] = v3d_mV3AbsKv($P1,$P2) * v3d_mV3AbsKv($P2,$P3) * v3d_mV3AbsKv($P3,$P1) / (2*$JmP)**2;
$Aprox['RadSq'] = $Aprox['RadSq'] + 0.0001; //Float rounding errors
return $Aprox;
}

function CircumscribedCircle_IsIn($Aprox, $Pnt)
{
return v3d_mV3AbsKv($Aprox['Center'], $Pnt) <= $Aprox['RadSq'];
}


function CircumscribedCircle_SetAprox($Aprox, $Pnt,
  $HrLen=0,$i=0,$TestAprox=[],$Ok=true,$j=0, $TestRet=[])
{
$HrLen=ar_Count($Aprox['Borders']); //2 or 3
for ($i=0; $i<$HrLen; $i++)
  {
  $TestAprox = CircumscribedCircle_2Points($Aprox['Borders'][$i], $Pnt);
  $Ok = true;
  for ($j=0; $j<$HrLen; $j++)
    {
    if ($i==$j) continue;
    if (CircumscribedCircle_IsIn($TestAprox, $Aprox['Borders'][$j])) continue;
    $Ok = false;
    break;
    }
  if ($Ok) return $TestAprox;
  }
if ($HrLen==2)
  {
  return CircumscribedCircle_3Points($Aprox['Borders'][0], $Aprox['Borders'][1], $Pnt);
  }
$TestRet = ar_New();
$TestAprox = CircumscribedCircle_3Points($Aprox['Borders'][0], $Aprox['Borders'][1], $Pnt);
if (CircumscribedCircle_IsIn($TestAprox, $Aprox['Borders'][2])) ar_Push($TestRet, $TestAprox);
$TestAprox = CircumscribedCircle_3Points($Aprox['Borders'][0], $Aprox['Borders'][2], $Pnt);
if (CircumscribedCircle_IsIn($TestAprox, $Aprox['Borders'][1])) ar_Push($TestRet, $TestAprox);
$TestAprox = CircumscribedCircle_3Points($Aprox['Borders'][1], $Aprox['Borders'][2], $Pnt);
if (CircumscribedCircle_IsIn($TestAprox, $Aprox['Borders'][0])) ar_Push($TestRet, $TestAprox);
ar_Sort($TestRet, function($a, $b){return $a['RadSq']-$b['RadSq'];});
return $TestRet[0];
}

function CircumscribedCircle($Arr,
  $Aprox=[], $i=0, $l=0, $max=0)
{
$Aprox = CircumscribedCircle_2Points($Arr[0], $Arr[1]);
for ($i=0, $l=ar_Count($Arr); $i<$l; $i++)
  {
  if (CircumscribedCircle_IsIn($Aprox, $Arr[$i])) continue;
  $Aprox = CircumscribedCircle_SetAprox($Aprox, $Arr[$i]);
  if (!$Aprox) return false;
  $i=-1;
  }
return as_New('Center',$Aprox['Center'], 'Radius', nm_Sqrt($Aprox['RadSq']));
}

//-------------------------
in JS
const x = 0;
const y = 1;
const z = 2;
in PHP
define('x', 0);
define('y', 1);
define('z', 2);

function v3d_AvgV3($v, $w)
//Diameter of two vectors (center of the line segment)
{
return ar_New(($v[x]+$w[x])/2, ($v[y]+$w[y])/2, ($v[z]+$w[z])/2);
}

function v3d_Abs($v)
//Absolute vector size
{
return nm_Sqrt($v[x]**2 + $v[y]**2 + $v[z]**2);
}

function v3d_vxV3($v,$w)
//Vector product
{
return ar_New($v[y]*$w[z]-$w[y]*$v[z], $v[z]*$w[x]-$w[z]*$v[x], $v[x]*$w[y]-$w[x]*$v[y]);
}

function v3d_mV3Abs($v, $w)
//Vector difference
{
return nm_Sqrt(($v[x]-$w[x])**2 + ($v[y]-$w[y])**2 + ($v[z]-$w[z])**2);
}

function v3d_mV3AbsKv($v, $w)
//Square of vector difference
{
return ($v[x]-$w[x])**2 + ($v[y]-$w[y])**2 + ($v[z]-$w[z])**2;
}

function v3d_sxV3($v,$w)
//Scalar product
{
return $v[x]*$w[x] + $v[y]*$w[y] + $v[z]*$w[z];
}

function v3d_pV3pV3($v,$w,$u)
//Sum 3 3D vect
{
return ar_New($v[x]+$w[x]+$u[x], $v[y]+$w[y]+$u[y], $v[z]+$w[z]+$u[z]);
}

function v3d_xSkl($v,$k)
//3D vect multiple scalar
{
return ar_New($v[x]*$k, $v[y]*$k, $v[z]*$k);
}

代码是通用的，可以在 PHP 和 Javascript 中使用 :-) 我将辅助变量声明为可选参数（以避免使用 var 关键字）

as_* 是与关联数组相关的函数，ar_* 是与数值数组相关的函数，nm_* 是数值型的，从名字就可以看出它们的功能。 v3d_* 是与 3d 向量相关的函数，列在脚注中。

用www.DeepL.com/Translator翻译（免费版）