求二维曲线双切线的算法

Question

我正在寻找一种算法来计算曲线的双切线，即。曲线下方的线恰好与它相交于两点：

具体来说，我正在寻找双切线与曲线相交的两个 x 值。这些是我正在研究的实际曲线的例子，我已经手动绘制了我想要找到的切线。

从 Mathematica Stack Exchange 几乎一字不漏地引用 this post，我知道我正在寻找两个不同的 x 值，通用行描述这条曲线的函数有

相同的 y 值（即直线与曲线相交）
同导数（即直线与曲线相切）

如果我们知道描述曲线的函数 f(x)，这就有效，因为我们知道通用直线的方程是 y = mx + b。然后我们可以建立以下方程组并求解：

y₁ = f(x₁),
y₂ = f(x₂),
f'(x₁) = f'(x₂),
f'(x₁) = (y₂ - y₁)/(x₂ - x₁)

我遇到的问题有两个。我不知道这条线的功能，因为它来自测量，即使我知道，我也不知道如何使用 C#/MathNet 求解方程组。关于曲线，我唯一可以明确地说的是双切线的斜率为正，x 轴将从 -150 到 700 运行。

我尝试过的其他方法是使用修改后的凸包算法，并通过手动选择结点来拟合三次样条，但这两种尝试都不够准确。

所以我的问题是：

如何找到这条曲线的函数（最好使用 MathNet/C#）
如何求解上述方程组（最好使用 MathNet/C#）

感谢任何其他算法想法/建议！

谢谢

Other related post 在 Mathematica Stack Exchange 上

Answer 1

由于函数在分析上未知，唯一可行的方法是数值计算。使用凸包的想法似乎很有希望。根据您的示例数据，它似乎应该按如下方式工作：

计算下凸包，例如使用 Graham Scan 实际上分别计算上下一个（或最终 left/right 一个取决于实现）
取较长的船体段
也许较长的是您正在寻找的，或者最终您可能需要采用一些较长的（例如 2 到 5）并根据一些合适的想法创建它们的排名，例如：
- 线段上的点与线段本身之间的最大 y 距离
- 走最中心的那个
- 取坡度高的那个

您需要对您的数据进行试验（并了解其物理含义），以了解什么是最适合select您需要的那一段的排名标准。

Answer 2

输入数据

因为你没有以数字形式分享原始输入数据，所以我从图像中提取了它。 这部分是为其他人准备的...首先我在绘画中对它们进行了一些编辑：

用黑色填充背景

删除了许多双切点
内部裁剪

删除轴和 scales/grid
手动清除副切线留下的残余像素

然后我应用一些 DIP 将所有垂直线和水平线缩小到一个点，结果：

然后我提取任何非背景像素的 xy（第一个数字是 x，第二个 y 以像素为单位）：

int data0_xy[]=
    {
      33, 447,  41, 451,  49, 462,  56, 470,  64, 473,  72, 477,  80, 481,  87, 489,  95, 492, 103, 495, 111, 496, 118, 500, 126, 504, 134, 507, 142, 506, 149, 506,
     157, 510, 165, 514, 173, 512, 180, 512, 188, 512, 196, 515, 204, 515, 212, 513, 219, 512, 227, 514, 235, 514, 243, 513, 250, 512, 258, 511, 266, 512, 274, 511,
     281, 508, 289, 505, 297, 505, 305, 506, 312, 505, 320, 500, 328, 500, 336, 500, 343, 501, 351, 497, 359, 495, 367, 497, 374, 499, 382, 500, 390, 498, 398, 498,
     406, 501, 413, 502, 421, 501, 429, 501, 437, 502, 444, 504, 452, 506, 460, 504, 468, 504, 475, 504, 483, 505, 491, 505, 499, 503, 506, 500, 514, 501, 522, 503,
     530, 499, 537, 497, 545, 495, 553, 496, 561, 495, 569, 490, 576, 486, 584, 486, 592, 485, 600, 483, 607, 478, 615, 475, 623, 472, 631, 471, 638, 468, 646, 461,
     654, 457, 662, 454, 669, 452, 677, 450, 685, 443, 693, 436, 700, 431, 708, 429, 716, 425, 724, 417, 732, 410, 739, 405, 747, 400, 755, 397, 763, 387, 770, 380,
     778, 375, 786, 369, 794, 363, 801, 355, 809, 346, 817, 342, 825, 336, 832, 326, 840, 318, 848, 310, 856, 306, 863, 300, 871, 291, 879, 281, 887, 274, 894, 269,
     902, 263, 910, 252, 918, 244, 926, 237, 933, 230, 941, 224, 949, 215, 957, 206, 964, 202, 972, 197, 980, 187, 988, 179, 995, 173,1003, 169,1011, 162,1019, 155,
    1026, 148,1034, 141,1042, 136,1050, 133,1057, 126,1065, 119,1073, 115,1081, 113,1089, 108,1096, 102,1104,  97,1112,  97,1120,  94,1127,  87,1135,  84,1143,  82,
    1151,  84,1158,  80,1166,  77,1174,  74,1182,  72,1189,  73,1197,  72,1205,  69,1213,  67,1220,  68,1228,  68,1236,  64,1244,  62,1252,  61,1259,  63,1267,  60,
    1275,  57,1283,  53,1290,  49,1298,  49,1306,  47,1314,  40,1321,  35,1329,  31,1337,  27,1345,  20,1352,  10,
    };
int data1_xy[]=
    {
      33, 533,  41, 533,  49, 532,  56, 531,  64, 533,  72, 532,  80, 530,  87, 533,  95, 530, 103, 533, 111, 531, 118, 532, 126, 531, 134, 531, 142, 531, 149, 529,
     157, 530, 165, 530, 173, 529, 180, 526, 188, 529, 196, 527, 204, 525, 212, 526, 219, 526, 227, 525, 235, 523, 243, 520, 250, 522, 258, 520, 266, 520, 274, 517,
     281, 519, 289, 516, 297, 515, 305, 513, 312, 511, 320, 510, 328, 509, 336, 507, 343, 503, 351, 503, 359, 501, 367, 500, 374, 496, 382, 496, 390, 493, 398, 490,
     406, 488, 413, 485, 421, 482, 429, 480, 437, 476, 444, 469, 452, 467, 460, 466, 468, 461, 475, 456, 483, 451, 491, 449, 499, 442, 506, 438, 514, 432, 522, 426,
     530, 420, 537, 411, 545, 402, 553, 400, 561, 395, 569, 386, 576, 380, 584, 372, 592, 363, 600, 354, 607, 344, 615, 336, 623, 327, 631, 323, 638, 315, 646, 300,
     654, 299, 662, 291, 669, 284, 677, 273, 685, 266, 693, 256, 700, 250, 708, 243, 716, 237, 724, 227, 732, 217, 739, 210, 747, 197, 755, 185, 763, 184, 770, 151,
     778, 165, 786, 158, 794, 145, 801, 136, 809, 111, 817, 119, 825, 112, 832, 105, 840,  94, 848,  76, 856,  80, 863,  78, 871,  66, 879,  60, 887,  40, 894,  28,
     902,  38, 910,  42, 918,  41, 926,  46, 933,  47, 941,  41, 949,  42, 957,  39, 964,  48, 972,  52, 980,  49, 988,  57, 995,  51,1003,  67,1011,  78,1019,  85,
    1026,  90,1034,  96,1042, 106,1050, 114,1057, 123,1065, 132,1073, 126,1081, 137,1089, 157,1096, 169,1104, 175,1112, 169,1120, 189,1127, 191,1135, 204,1143, 209,
    1151, 221,1158, 229,1166, 235,1174, 226,1182, 241,1189, 244,1197, 257,1205, 264,1213, 260,1220, 273,1228, 275,1236, 280,1244, 285,1252, 288,1259, 290,1267, 292,
    1275, 290,1283, 297,1290, 298,1298, 298,1306, 301,1314, 300,1321, 304,1329, 304,1337, 298,1345, 297,1352, 297,
    };

所以现在我们拥有与您大致相同的数据...

双切线

因为你的数据不是那么大，你可以尝试暴力破解方法。

遍历所有点对

简单O(n^2)搜索。不要测试点对两次，这样第二个循环就不会测试在第一个循环中已经完成的点。
检测双切线是否有效

因此，如果我们计算测试双切线的法线向量（指向数据点的垂直向量）。然后所有向量 p(i)-p(bitangent) 必须与我们的法线具有非负（或者正，如果你想要恰好 2 个命中）点积。这导致另一个 O(n) 循环导致 O(n^3) 方法。如果任何点积过零，则丢弃该双切线并测试下一个。如果没有点积交叉，您会发现双切线，因此将前两个循环的 actaul 点作为新的双切线添加到列表中。

这将找到数据较低侧的所有双切线（如果翻转法向量或交叉条件，则在其上方）。所以现在你可以将你的启发式应用到你想要的select双切线。

我不会用 C# 编写代码，所以这里是简单的 C++ 示例：

const int n2=sizeof(data_xy)/sizeof(data_xy[0]);    // numbers in data
const int n=n2>>1;                                  // points in data
int bitangent[100],bitangents;                      // 2 poit indexes per bitangent, number of indexes in bitangent[]
// O(n^3) bruteforce bitangent search
int nx,ny,i2,j2,k2,dx,dy;
bitangents=0;
for (i2=0;i2<n2;i2+=2)      // loop all points (bitangent start)
 for (j2=i2+2;j2<n2;j2+=2)  // loop all yet untested points (bitangent end)
    {
    // normal
    ny=-(data_xy[j2+0]-data_xy[i2+0]);
    nx=+(data_xy[j2+1]-data_xy[i2+1]);
    // test overlap
    for (k2=0;k2<n2;k2+=2)
     if ((k2!=i2)&&(k2!=j2))
        {
        dx=data_xy[k2+0]-data_xy[i2+0];
        dy=data_xy[k2+1]-data_xy[i2+1];
        if ((dx*nx)+(dy*ny)<0) { k2=-1; break; } // if dot product is negative overlap occurs so throw solution away
        }
    if (k2>=0)
        {
        bitangent[bitangents]=i2; bitangents++;
        bitangent[bitangents]=j2; bitangents++;
        }
    }

我这样渲染（VCL）：

void draw()
    {
    int i2,j2,x,y,r=3;
    bmp->Canvas->Brush->Color=clWhite;
    bmp->Canvas->FillRect(TRect(0,0,xs,ys));

    // render data lines
    bmp->Canvas->Pen->Color=clBlack;
    for (i2=0;i2<n2;i2+=2)
        {
        x=data_xy[i2+0];
        y=data_xy[i2+1];
        if (!i2) bmp->Canvas->MoveTo(x,y);
        else     bmp->Canvas->LineTo(x,y);
        }

    // render bitangents lines
    bmp->Canvas->Pen->Color=clBlue;
    for (i2=0;i2<bitangents;i2+=2)
        {
        j2=bitangent[i2+0];
        x=data_xy[j2+0];
        y=data_xy[j2+1];
        bmp->Canvas->MoveTo(x,y);
        j2=bitangent[i2+1];
        x=data_xy[j2+0];
        y=data_xy[j2+1];
        bmp->Canvas->LineTo(x,y);
        }

    // render data points
    bmp->Canvas->Pen->Color=clBlack;
    bmp->Canvas->Brush->Color=clRed;
    for (i2=0;i2<n2;i2+=2)
        {
        x=data_xy[i2+0];
        y=data_xy[i2+1];
        bmp->Canvas->Ellipse(x-r,y-r,x+r,y+r);
        }
    // render bitangents points
    bmp->Canvas->Pen->Color=clBlack;
    bmp->Canvas->Brush->Color=clAqua;
    for (i2=0;i2<bitangents;i2++)
        {
        j2=bitangent[i2];
        x=data_xy[j2+0];
        y=data_xy[j2+1];
        bmp->Canvas->Ellipse(x-r,y-r,x+r,y+r);
        }

    Form1->Canvas->Draw(0,0,bmp);
    }

此处为您的数据输出：

如您所见，我不需要任何浮点运算或变量:) 为了尽可能简单，我使用了静态数组（没有动态分配）。如您所见，双切线比您最初想象的要多得多。

求二维曲线双切线的算法

Algorithm to find the bitangent of a 2D curve

c#

algorithm

mathnet-numerics