Search code examples
c#algorithmmathnet-numerics

Algorithm to find the bitangent of a 2D curve

I'm looking for an algorithm to calculate the bitangent of a curve, ie. the line below a curve that intersects it at exactly two points: Example 1 Example 2

Specifically, I am looking for the two x-values where the bitangent intersects the curve. These are examples of actual curves that I am studying, I have manually drawn the tangent I would like to find.

Referencing this post from the Mathematica Stack Exchange almost verbatim, I understand that I am looking for two distinct x values for which a generic line and the function describing this curve have

  1. The same y value (ie. the line touches the curve)
  2. The same derivative (ie. the line is tangent to the curve)

This works if we know the function f(x) that describes the curve, because we know the equation of a generic line is y = mx + b. We can then set up the following system of equations and solve:

  1. y1 = f(x1),

  2. y2 = f(x2),

  3. f'(x1) = f'(x2),

  4. f'(x1) = (y2 - y1)/(x2 - x1)

The issue I am having is two fold. I don't know the function of the line, as it comes from a measurement, and even if I did, I don't know how to solve a system of equations using C# / MathNet. The only thing I can say definitively about the curve is the slope of the bitangent will be positive, and the x-axis will run from -150 to 700.

Other things I've tried is using a modified convex hull algorithm, and fitting a cubic spline by manually selecting the knot points, but both these attempts were insufficiently accurate.

So my questions are:

  1. How can I find the function of this curve (ideally using MathNet / C#)
  2. How can I solve the above system of equations (also ideally using MathNet / C#)

Any other algorithm ideas / advice is appreciated!

Thanks

Other related posts on the Mathematica Stack Exchange

Solution

  • Input data

    Well as you did not shared the original input data in numerical form I extracted it from the images. This part is meant for others... First I edited them in paint a bit:

    1. fill the background by Black color

      that removes many of the bi tangent points

    2. crop inside

      that removes axises and scales/grid

    3. manually clear the remnant pixels of what has been left from the bitangent

    Then I apply some DIP to shrink all vertical and horizontal lines to single point here the results:

    plot0 plot1

    and then I extracted xy of any non background pixels (first number is x, second y in pixels):

    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,
    };

    So now we have more or less the same data as you do...

    Bitangent

    Well as your data is not that big you can try bruteforce approach.

    1. loop through all point pairs

      simple O(n^2) search. Do not test point pairs twice so the second loop will not test points already done in first loop.

    2. detect if valid bitangent

      so if we compute normal vector to tested bitangent (perpendicular vector to it pointing towards data points). Then all the vectors p(i)-p(bitangent) must have non negative (or positive if you want exactly 2 hits) dot product with our normal. This lead to another O(n) loop resulting in O(n^3) approach. If any dot product crosses zero throw away this bitangent and test the next one. If no dot product crosses you found bitangent so add the actaul points from first two loops as a new bitangent to the list.

    This will find all the bitangents on the lowwer side of data (or above it if you flip the normal vector or crossing condition). So now you can apply your heuristics to select bitangent you want.

    I do not code in C# so here simple C++ example:

    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++;
        }
    }

    I render like this (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);
    }

    And here output for your data:

    output0 output1

    As you can see I did not need any floating point operation or variable :) To keep this as simple as possible I used static arrays (no dynamic allocation). As you can see there are quite more bitangents than you think at first.