悠闲数学娱乐论坛(第2版) » Mathematica » 寻找图形中的共圆点

hbghlyj 发表于 2021-6-30 01:53

寻找图形中的共圆点

[i=s] 本帖最后由 hbghlyj 于 2021-6-30 19:40 编辑 [/i]

算法1.用内置的函数CircleThrough求三个点的外接圆,用判断第四个点是否在圆上.
算法2.
对三个点$(x_i, y_i),i=1, 2, 3$解线性方程组$2ax_i+2by_i+c=x_i^2+y_i^2, i=1,2,3$,如果无解就Continue[],如果有解,就取另外一个点(x,y),检验$2ax+2by+c=x^2+y^2$是否成立.
算法3.
计算行列式$\det\begin{bmatrix}
     x^2+y^2 & x & y & 1 \\
     x_0^2+y_0^2 & x_0 & y_0 & 1 \\
     x_1^2+y_1^2 & x_1 & y_1 & 1 \\
     x_2^2+y_2^2 & x_2 & y_2 & 1 \\
   \end{bmatrix} = 0
$
算法4.托勒密定理
算法5.对i循环,将$P_{i+1}-P_i,P_{i+2}-P_i,\cdots,P_n-P_i$关于单位圆反演,寻找共线的点.
算法6.有向角$P_1P_3P_2=P_1P_iP_2$
算法7.直线$l_{12}$与$l_{3i}$的倾斜角之和=$l_{13}l_{2i}$的倾斜角之和$(±π)$

hbghlyj 发表于 2021-6-30 03:18

[i=s] 本帖最后由 hbghlyj 于 2021-6-30 14:45 编辑 [/i]

以下面这组点为例:
p = {{0.07059703, 2.3311955963}, {-0.2312351881,
   1.4561916933}, {1.6323255814, 0.4529926561}, {1.8817163479,
   2.3923914504}, {0.274158443, 1.861522012}, {1.3099706592,
   1.5983182696}, {-0.633824898, 1.7300548894}, {2.6092295121,
   0.0762685712}}
----------
算法1.
AbsoluteTiming[Module[{}, S = {};
  For[i = 1, i <= Length[p] - 3, i++,
   For[j = i + 1, j <= Length[p] - 2, j++,
    For[k = j + 1, k <= Length[p] - 1,
     k++, {center, radius} = List @@ CircleThrough[p[[{i, j, k}]]];
     For[l = k + 1, l <= Length[p], l++,
      If[Abs[EuclideanDistance[p[[l]], center] - radius] < 1/10^3,
       AppendTo[S, {i, j, k, l}]]]]]];
  S /. Thread[{1, 2, 3, 4, 5, 6, 7, 8} -> {"a", "b", "c", "d", "e",
      "f", "g", "h"}]]]
输出为
{0.0021478, {{"a", "b", "c", "d"}, {"a", "b", "e", "g"}, {"a", "c",
   "e", "h"}, {"a", "d", "e", "f"}, {"b", "c", "g", "h"}, {"b", "d",
   "f", "g"}, {"c", "d", "f", "h"}, {"e", "f", "g", "h"}}}
测试多次,实际用时0.002~0.003秒
-----------
可以看出,虽然结果正确,但是计算过程中会产生很大的误差,p中的坐标都精确到$10^{-10}$,但是作为检验的点到圆心的距离与半径之差的阈值却只能设置为$10^{-3}$(我用$10^{-4}$试过一次,会漏掉两个圆)

hbghlyj 发表于 2021-6-30 03:45

算法4.
AbsoluteTiming[
Module[{},
  d = Association[
    Table[{s, t} -> EuclideanDistance[p[[s]], p[[t]]], {s,
      Length[p] - 1}, {t, s + 1, Length[p]}]]; S = {};
  For[i = 1, i <= Length[p] - 3, i++,
   For[j = i + 1, j <= Length[p] - 2, j++, ij = d[{i, j}];
    For[k = j + 1, k <= Length[p] - 1, k++, ik = d[{i, k}];
     jk = d[{j, k}];
     For[l = k + 1, l <= Length[p], l++, il = d[{i, l}];
      jl = d[{j, l}]; kl = d[{k, l}];
      If[Abs[(il jk - ik jl - ij kl) (il jk + ik jl - ij kl) (il jk -
            ik jl + ij kl) (il jk + ik jl + ij kl)] < 1/10^3,
       AppendTo[S, {i, j, k, l}]]]]]];
  S /. Thread[{1, 2, 3, 4, 5, 6, 7, 8} -> {"a", "b", "c", "d", "e",
      "f", "g", "h"}]]]
输出为
{0.000824, {{"a", "b", "c", "d"}, {"a", "b", "e", "g"}, {"a", "c",
   "e", "h"}, {"a", "d", "e", "f"}, {"b", "c", "g", "h"}, {"b", "d",
   "f", "g"}, {"c", "d", "f", "h"}, {"e", "f", "g", "h"}}}
测试多次,实际用时0.0008~0.0009秒

hbghlyj 发表于 2021-6-30 03:53

算法3.
AbsoluteTiming[
Module[{}, S = {};
  For[i = 1, i <= Length[p] - 3, i++, {x0, y0} = p[[i]];
   For[j = i + 1, j <= Length[p] - 2, j++, {x1, y1} = p[[j]];
    For[k = j + 1, k <= Length[p] - 1, k++, {x2, y2} = p[[k]];
     f[{x_, y_}] =
      Det[{{x^2 + y^2, x, y, 1}, {x0^2 + y0^2, x0, y0,
         1}, {x1^2 + y1^2, x1, y1, 1}, {x2^2 + y2^2, x2, y2, 1}}];
     For[l = k + 1, l <= Length[p], l++,
      If[Abs[f[p[[l]]]] < 1/10^3, AppendTo[S, {i, j, k, l}]]]]]];
  S /. Thread[{1, 2, 3, 4, 5, 6, 7, 8} -> {"a", "b", "c", "d", "e",
      "f", "g", "h"}]]]
输出为
{0.0042354, {{"a", "b", "c", "d"}, {"a", "b", "e", "g"}, {"a", "c",
   "e", "h"}, {"a", "d", "e", "f"}, {"b", "c", "g", "h"}, {"b", "d",
   "f", "g"}, {"c", "d", "f", "h"}, {"e", "f", "g", "h"}}}
测试多次,实际用时0.004~0.005秒

hbghlyj 发表于 2021-6-30 04:12

算法2.
AbsoluteTiming[
Module[{}, S = {};
  For[i = 1, i <= Length[p] - 3, i++, {x0, y0} = p[[i]];
   For[j = i + 1, j <= Length[p] - 2, j++, {x1, y1} = p[[j]];
    For[k = j + 1, k <= Length[p] - 1,
     k++, {x2, y2} = p[[k]]; {a, b, c} =
      LinearSolve[{{2 x0, 2 y0, 1}, {2 x1, 2 y1, 1}, {2 x2, 2 y2,
         1}}, {x0^2 + y0^2, x1^2 + y1^2, x2^2 + y2^2}];
     For[l = k + 1, l <= Length[p], l++, {x3, y3} = p[[l]];
      If[Abs[2 a x3 + 2 b y3 + c - x3^2 - y3^2] < 1/10^3,
       AppendTo[S, {i, j, k, l}]]]]]];
  S /. Thread[{1, 2, 3, 4, 5, 6, 7, 8} -> {"a", "b", "c", "d", "e",
      "f", "g", "h"}]]]
输出为
{0.0008629, {{"a", "b", "c", "d"}, {"a", "b", "e", "g"}, {"a", "c",
   "e", "h"}, {"a", "d", "e", "f"}, {"b", "c", "g", "h"}, {"b", "d",
   "f", "g"}, {"c", "d", "f", "h"}, {"e", "f", "g", "h"}}}
测试多次,实际用时0.0008~0.0009秒

hbghlyj 发表于 2021-6-30 14:41

[i=s] 本帖最后由 hbghlyj 于 2021-6-30 19:38 编辑 [/i]

算法6.
AbsoluteTiming[Module[{}, S = {};
  f[{x1_, y1_}, {x2_, y2_}] = (y1 x2 - y2 x1)/(x1 x2 + y1 y2);
  For[i = 1, i <= Length[p] - 3, i++,
   For[j = i + 1, j <= Length[p] - 2, j++, ij = p[[i]] - p[[j]];
    For[k = j + 1, k <= Length[p] - 1, k++,
     angle = f[ij, p[[k]] - p[[j]]];
     For[l = k + 1, l <= Length[p], l++,
      If[Abs[f[p[[i]] - p[[l]], p[[k]] - p[[l]]] - angle] < 1/10^3,
       AppendTo[S, {i, j, k, l}]]]]]];
  S /. Thread[{1, 2, 3, 4, 5, 6, 7, 8} -> {"a", "b", "c", "d", "e",
      "f", "g", "h"}]]]
输出为
{0.0008468, {{"a", "b", "c", "d"}, {"a", "b", "e", "g"}, {"a", "c",
   "e", "h"}, {"a", "d", "e", "f"}, {"b", "c", "g", "h"}, {"b", "d",
   "f", "g"}, {"c", "d", "f", "h"}, {"e", "f", "g", "h"}}}
测试多次,实际用时0.0008~0.0009秒

hbghlyj 发表于 2021-6-30 15:10

[i=s] 本帖最后由 hbghlyj 于 2021-6-30 15:46 编辑 [/i]

算法7.
AbsoluteTiming[
Module[{}, f[{x_, y_}] = ArcTan[y/x];
  d = Association[
    Table[{s, t} -> f[p[[s]] - p[[t]]], {s, Length[p] - 1}, {t, s + 1,
       Length[p]}]]; S = {};
  For[i = 1, i <= Length[p] - 3, i++,
   For[j = i + 1, j <= Length[p] - 2, j++, ij = d[{i, j}];
    For[k = j + 1, k <= Length[p] - 1, k++, ik = d[{i, k}];
     For[l = k + 1, l <= Length[p], l++,
      jl = d[{j, l}]; kl = d[{k, l}];
      If[Abs[ij + kl - ik - jl] < 1/10^3 ||
        Abs[ij + kl - ik - jl - Pi] < 1/10^3 ||
        Abs[ij + kl - ik - jl + Pi] < 1/10^3,
       AppendTo[S, {i, j, k, l}]]]]]];
  S /. Thread[{1, 2, 3, 4, 5, 6, 7, 8} -> {"a", "b", "c", "d", "e",
      "f", "g", "h"}]]]
输出为
{0.0009069, {{"a", "b", "c", "d"}, {"a", "b", "e", "g"}, {"a", "c",
   "e", "h"}, {"a", "d", "e", "f"}, {"b", "c", "g", "h"}, {"b", "d",
   "f", "g"}, {"c", "d", "f", "h"}, {"e", "f", "g", "h"}}}
测试多次,实际用时0.0008~0.0009秒

hbghlyj 发表于 2021-6-30 18:27

[i=s] 本帖最后由 hbghlyj 于 2021-6-30 19:34 编辑 [/i]

算法5.其中使用了[url=http://kuing.orzweb.net/viewthread.php?tid=8035&page=1&extra=#pid40235]这帖[/url]的算法3.[code]
FindDuplicates =
Function[x, Module[{counts}, counts = Select[Counts[x], # > 1 &];
   Flatten[Position[x, #]] & /@ Keys[counts]]];
FindCollinear = Function[q, Module[{S}, S = {}; g[{x_, y_}] = y/x;
   For[i = 1, i <= Length[q] - 2, i++, pi = q[[i]];
    dup =
     FindDuplicates[
      Table[Round[g[q[[j]] - pi], 10^(-3)], {j, i + 1, Length[q]}]];
    If[dup == {}, , S = S~Join~(Prepend[i] /@ (i + dup))]]; S]];
AbsoluteTiming[Module[{}, S = {}; f[{x_, y_}] = {x, y}/(x^2 + y^2);
  For[j = 1, j <= Length[p] - 3, j++,
   coll = FindCollinear[f[# - p[[j]]] & /@ p[[j + 1 ;;]]];
   If[coll != {}, S = S~Join~(Prepend[j] /@ (j + coll))]];
  S /. Thread[{1, 2, 3, 4, 5, 6, 7, 8} -> {"a", "b", "c", "d", "e",
      "f", "g", "h"}]]]
[/code]输出为
{0.000752, {{"a", "b", "c", "d"}, {"a", "b", "e", "g"}, {"a", "c",
   "e", "h"}, {"a", "d", "e", "f"}, {"b", "c", "g", "h"}, {"b", "d",
   "f", "g"}, {"c", "d", "f", "h"}, {"e", "f", "g", "h"}}}
测试多次,实际用时0.0007~0.0008秒

hbghlyj 发表于 2021-6-30 19:34

结论:算法5应该是最快的

查看完整版本: 寻找图形中的共圆点