Алгоритм нахождения битангенса 2D кривой - PullRequest
2 голосов
/ 18 апреля 2020

Я ищу алгоритм для вычисления битангента кривой, ie. линия под кривой, которая пересекает ее ровно в двух точках: Example 1 Example 2

В частности, я ищу два значения x , где битангенс пересекает кривую. Это примеры реальных кривых, которые я изучаю, я вручную нарисовал касательную, которую я хотел бы найти.

Ссылка этот пост от Mathematica Stack Exchange почти дословно, я понимаю, что я ищу два различных значения x , для которых общая строка c и функция, описывающая эту кривую, имеют

  1. То же значение y (ie. линия касается кривой)
  2. Та же производная (ie. линия касается кривой)

Это работает, если мы знаем функцию f (x) , которая описывает кривую, потому что мы знаем, что уравнение общей линии c равно y = mx + b . Затем мы можем установить следующую систему уравнений и решить:

  1. y 1 = f (x 1 ) ,

  2. y 2 = f (x 2 ) ,

  3. f '(x 1 ) = f' (x 2 ) ,

  4. f '(x 1 ) = (y 2 - y 1 ) / (x 2 - x 1 )

У меня проблема в два раза. Я не знаю функцию линии, как это происходит при измерении, и даже если бы я знал, я не знаю, как решить систему уравнений, используя C# / Mat hNet. Единственное, что я могу сказать однозначно о кривой, это то, что наклон битангента будет положительным, а ось x будет работать в диапазоне от -150 до 700.

Другие вещи, которые я пробовал, - использование модифицированного выпуклого алгоритм корпуса и подгонка сплайна cubi c путем ручного выбора узловых точек, но обе эти попытки были недостаточно точными.

Итак, мои вопросы:

  1. Как мне найти функция этой кривой (в идеале, используя Mat hNet / C#)
  2. Как я могу решить вышеуказанную систему уравнений (также в идеале, используя Mat hNet / C#)

Любые другие идеи и советы алгоритма приветствуются!

Спасибо

Другое связанных сообщений на бирже Mathematica Stack Exchange

Ответы [ 2 ]

2 голосов
/ 19 апреля 2020

Поскольку функция не известна аналитически, единственный жизнеспособный подход работает численно. Идея использования выпуклой оболочки кажется многообещающей. Основываясь на данных вашего примера, кажется, что это должно работать следующим образом:

  • вычисляет нижнюю выпуклую оболочку, например, с помощью Graham Scan , которая фактически вычисляет отдельно верхнюю и нижнюю (или в конечном итоге) левый / правый в зависимости от реализации)
  • возьмите более длинные сегменты корпуса
  • возможно, чем длиннее тот, который вы ищете, или в конечном итоге вам может понадобиться некоторое количество более длинных ( например, от 2 до 5), и создайте их ранжирование на основе некоторой подходящей идеи, например:
    • максимальное расстояние по оси y между точками на отрезке и самим отрезком
    • выберите самое центральное
    • возьмите тот, у которого более высокий наклон

Вам необходимо поэкспериментировать с вашими данными (и знать их физический смысл), чтобы увидеть, какой рейтинг наиболее подходит критерии выбора одного сегмента, который вам нужен, из более длинных.

1 голос
/ 20 апреля 2020

Входные данные

Ну, так как вы не передали исходные входные данные в числовой форме, я извлек их из изображений. Эта часть предназначена для других ... Сначала я немного их отредактировал:

  1. закрасить фон черным цветом

    , который удаляет многие из точек касательных би

  2. кадрирование внутри

    , который удаляет оси и шкалы / сетку

  3. вручную очистить оставшиеся пиксели того, что осталось от битангента

Затем я применяю немного DIP к сожмите все вертикальные и горизонтальные линии в одну точку, вот результаты:

plot0 plot1

и затем я извлек 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,
    };

Так что теперь у нас более или менее те же данные, что и у вас ...

Битангенс

Ну, так как ваши данные не так велики, вы можете попробовать грубый подход.

  1. l oop через все пары точек

    простой O(n^2) поиск. Не проверяйте пары точек дважды, чтобы вторая l oop не проверяла точки, уже выполненные в первой l oop.

  2. , определите, является ли действительный битангенс

    , поэтому, если мы вычислим вектор нормали к проверенному битангенсу (перпендикулярный вектор к нему указывает на точки данных). Тогда все векторы p(i)-p(bitangent) должны иметь неотрицательное (или положительное, если вы хотите ровно 2 попадания) точечное произведение с нашим нормальным. Это привело к другому O(n) l oop, что привело к O(n^3) подходу. Если любое точечное произведение пересекает ноль, отбросьте этот битангенс и протестируйте следующий. Если ни одно точечное произведение не пересекает найденный битангенс, добавьте в список точки активации из первых двух петель как новый битангенс.

Это позволит найти все битангенты на нижней стороне данных ( или выше, если вы переверните нормальный вектор или условие пересечения). Так что теперь вы можете применить свою эвристику, чтобы выбрать нужный битангенс.

Я не кодирую в 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);
    }

А вот вывод для ваших данных:

output0 output1

Как видите, мне не нужны никакие операции или переменные с плавающей запятой :) Чтобы сделать это как можно более простым, я использовал массивы stati c (без динамического выделения c). Как видите, битангенс гораздо больше, чем вы думаете на первый взгляд.

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...