numerical instability in collinear segments

README.md

Results in a massive dx2 value which breaks return values for dx and dy:

Output is:

$ ./calc-intersection-test
m_width * (v1.y - v0.y): 0.200000000023408
m_width * (v2.y - v1.y): 0.300000000000018
dx2: 9771986970684.640625000000000
1 9771986970692.2246093750000000 8934388088183.3027343750000000

calc_intersection_test.cpp

//#include "agg_math.h"
#include "agg_vertex_sequence.h"
#include <iostream>
#include <iomanip>
/*
clang++ -o calc-intersection-test calc_intersection_test.cpp -Ideps/agg/include 
*/
namespace agg {
    bool calc_intersection_(double ax, double ay, double bx, double by,
                            double cx, double cy, double dx, double dy,
                            double* x, double* y)
    {
        double num = (ay-cy) * (dx-cx) - (ax-cx) * (dy-cy);
        double den = (bx-ax) * (dy-cy) - (by-ay) * (dx-cx);
        if(fabs(den) < intersection_epsilon) return false;
        double r = num / den;
        *x = ax + r * (bx-ax);
        *y = ay + r * (by-ay);
        return true;
    }
}

int main() {
  
    agg::vertex_dist v0(7.146342662755522,82.90081882764909);
    agg::vertex_dist v1(7.583842662776501,83.30081882769591);
    agg::vertex_dist v2(7.583842662776513,83.90081882769594);
    double len1 = 1.0915155748371188;
    double len2 = 0.0000000000000307;
    double m_width = 0.5000000000000000;
    double dx1 = m_width * (v1.y - v0.y) / len1;
    double dy1 = m_width * (v1.x - v0.x) / len1;
    double dx2 = m_width * (v2.y - v1.y) / len2;
    double dy2 = m_width * (v2.x - v1.x) / len2;
    std::clog << std::setprecision(15) << std::fixed << "m_width * (v1.y - v0.y): " << (double)(m_width * (v1.y - v0.y)) << "\n";
    std::clog << std::setprecision(15) << std::fixed << "m_width * (v2.y - v1.y): " << (double)(m_width * (v2.y - v1.y)) << "\n";
    std::clog << std::setprecision(15) << std::fixed << "dx2: " << dx2 << "\n"; // way to big!!!!
    double dx = (dx1 + dx2) / 2;
    double dy = (dy1 + dy2) / 2;
    bool ret = agg::calc_intersection_(v0.x + dx1, v0.y - dy1,
                                       v1.x + dx1, v1.y - dy1,
                                       /*here*/v1.x + dx2/*here*/, v1.y - dy2,
                                       v2.x + dx2, v2.y - dy2,
                                       &dx, &dy);
    std::clog << ret << " " << std::setprecision(16) << std::fixed << dx << " " << dy << "\n";
    return 0;
}