kicad/polygon/math_for_graphics.cpp

// math for graphics utility routines and RC, from FreePCB

#include <vector>

#include <math.h>
#include <float.h>
#include <limits.h>

#include "fctsys.h"

#include "PolyLine.h"

using namespace std;


/**
 * Function TestLineHit
 * test for hit on line segment i.e. a point within a given distance from segment
 * @param x, y = cursor coords
 * @param xi,yi and xf,yf = the end-points of the line segment
 * @param dist = maximum distance for hit
 * return true if dist < distance between the point and the segment
 */
bool TestLineHit( int xi, int yi, int xf, int yf, int x, int y, double dist )
{
    double dd;

    // test for vertical or horizontal segment
    if( xf==xi )
    {
        // vertical segment
        dd = fabs( (double)(x-xi) );
        if( dd<dist && ( (yf>yi && y<yf && y>yi) || (yf<yi && y>yf && y<yi) ) )
            return true;
    }
    else if( yf==yi )
    {
        // horizontal segment
        dd = fabs( (double)(y-yi) );
        if( dd<dist && ( (xf>xi && x<xf && x>xi) || (xf<xi && x>xf && x<xi) ) )
            return true;
    }
    else
    {
        // oblique segment
        // find a,b such that (xi,yi) and (xf,yf) lie on y = a + bx
        double b = (double)(yf-yi)/(xf-xi);
        double a = (double)yi-b*xi;
        // find c,d such that (x,y) lies on y = c + dx where d=(-1/b)
        double d = -1.0/b;
        double c = (double)y-d*x;
        // find nearest point to (x,y) on line segment (xi,yi) to (xf,yf)
        double xp = (a-c)/(d-b);
        double yp = a + b*xp;
        // find distance
        dd = sqrt((x-xp)*(x-xp)+(y-yp)*(y-yp));
        if( fabs(b)>0.7 )
        {
            // line segment more vertical than horizontal
            if( dd<dist && ( (yf>yi && yp<yf && yp>yi) || (yf<yi && yp>yf && yp<yi) ) )
                return 1;
        }
        else
        {
            // line segment more horizontal than vertical
            if( dd<dist && ( (xf>xi && xp<xf && xp>xi) || (xf<xi && xp>xf && xp<xi) ) )
                return true;
        }
    }
    return false;	// no hit
}


// set EllipseKH struct to describe the ellipse for an arc
//
int MakeEllipseFromArc( int xi, int yi, int xf, int yf, int style, EllipseKH * el )
{
    // arc (quadrant of ellipse)
    // convert to clockwise arc
    int xxi, xxf, yyi, yyf;
    if( style == CPolyLine::ARC_CCW )
    {
        xxi = xf;
        xxf = xi;
        yyi = yf;
        yyf = yi;
    }
    else
    {
        xxi = xi;
        xxf = xf;
        yyi = yi;
        yyf = yf;
    }
    // find center and radii of ellipse
    double xo=0, yo=0;
    if( xxf > xxi && yyf > yyi )
    {
        xo = xxf;
        yo = yyi;
        el->theta1 = M_PI;
        el->theta2 = M_PI/2.0;
    }
    else if( xxf < xxi && yyf > yyi )
    {
        xo = xxi;
        yo = yyf;
        el->theta1 = -M_PI/2.0;
        el->theta2 = -M_PI;
    }
    else if( xxf < xxi && yyf < yyi )
    {
        xo = xxf;
        yo = yyi;
        el->theta1 = 0.0;
        el->theta2 = -M_PI/2.0;
    }
    else if( xxf > xxi && yyf < yyi )
    {
        xo = xxi;
        yo = yyf;
        el->theta1 = M_PI/2.0;
        el->theta2 = 0.0;
    }
    el->Center.X = xo;
    el->Center.Y = yo;
    el->xrad = abs(xf-xi);
    el->yrad = abs(yf-yi);
#if 0
    el->Phi = 0.0;
    el->MaxRad = el->xrad;
    el->MinRad = el->yrad;
    if( el->MaxRad < el->MinRad )
    {
        el->MaxRad = el->yrad;
        el->MinRad = el->xrad;
        el->Phi = M_PI/2.0;
    }
#endif
    return 0;
}

// find intersections between line segment (xi,yi) to (xf,yf)
// and line segment (xi2,yi2) to (xf2,yf2)
// the line segments may be arcs (i.e. quadrant of an ellipse) or straight
// returns number of intersections found (max of 2)
// returns coords of intersections in arrays x[2], y[2]
//
int FindSegmentIntersections( int xi, int yi, int xf, int yf, int style,
                                 int xi2, int yi2, int xf2, int yf2, int style2,
                                 double x[], double y[] )
{
    double xr[12], yr[12];
    int iret = 0;

    if( max(xi,xf) < min(xi2,xf2)
        || min(xi,xf) > max(xi2,xf2)
        || max(yi,yf) < min(yi2,yf2)
        || min(yi,yf) > max(yi2,yf2) )
        return 0;

    if( style != CPolyLine::STRAIGHT && style2 != CPolyLine::STRAIGHT )
    {
        // two identical arcs intersect
        if( style == style2 && xi == xi2 && yi == yi2 && xf == xf2 && yf == yf2 )
        {
            if( x && y )
            {
                x[0] = xi;
                y[0] = yi;
            }
            return 1;
        }
        else if( style != style2 && xi == xf2 && yi == yf2 && xf == xi2 && yf == yi2 )
        {
            if( x && y )
            {
                x[0] = xi;
                y[0] = yi;
            }
            return 1;
        }
    }

    if( style == CPolyLine::STRAIGHT && style2 == CPolyLine::STRAIGHT )
    {
        // both straight-line segments
        int x, y;
        bool bYes = TestForIntersectionOfStraightLineSegments( xi, yi, xf, yf, xi2, yi2, xf2, yf2, &x, &y );
        if( !bYes )
            return 0;
        xr[0] = x;
        yr[0] = y;
        iret = 1;
    }
    else if( style == CPolyLine::STRAIGHT )
    {
        // first segment is straight, second segment is an arc
        int ret;
        double x1r, y1r, x2r, y2r;
        if( xf == xi )
        {
            // vertical first segment
            double a = xi;
            double b = DBL_MAX/2.0;
            ret = FindLineSegmentIntersection( a, b, xi2, yi2, xf2, yf2, style2,
                &x1r, &y1r, &x2r, &y2r );
        }
        else
        {
            double b = (double)(yf-yi)/(double)(xf-xi);
            double a = yf - b*xf;
            ret = FindLineSegmentIntersection( a, b, xi2, yi2, xf2, yf2, style2,
                &x1r, &y1r, &x2r, &y2r );
        }
        if( ret == 0 )
            return 0;
        if( InRange( x1r, xi, xf ) && InRange( y1r, yi, yf ) )
        {
            xr[iret] = x1r;
            yr[iret] = y1r;
            iret++;
        }
        if( ret == 2 )
        {
            if( InRange( x2r, xi, xf ) && InRange( y2r, yi, yf ) )
            {
                xr[iret] = x2r;
                yr[iret] = y2r;
                iret++;
            }
        }
    }
    else if( style2 == CPolyLine::STRAIGHT )
    {
        // first segment is an arc, second segment is straight
        int ret;
        double x1r, y1r, x2r, y2r;
        if( xf2 == xi2 )
        {
            // vertical second segment
            double a = xi2;
            double b = DBL_MAX/2.0;
            ret = FindLineSegmentIntersection( a, b, xi, yi, xf, yf, style,
                &x1r, &y1r, &x2r, &y2r );
        }
        else
        {
            double b = (double)(yf2-yi2)/(double)(xf2-xi2);
            double a = yf2 - b*xf2;
            ret = FindLineSegmentIntersection( a, b, xi, yi, xf, yf, style,
                &x1r, &y1r, &x2r, &y2r );
        }
        if( ret == 0 )
            return 0;
        if( InRange( x1r, xi2, xf2 ) && InRange( y1r, yi2, yf2 ) )
        {
            xr[iret] = x1r;
            yr[iret] = y1r;
            iret++;
        }
        if( ret == 2 )
        {
            if( InRange( x2r, xi2, xf2 ) && InRange( y2r, yi2, yf2 ) )
            {
                xr[iret] = x2r;
                yr[iret] = y2r;
                iret++;
            }
        }
    }
    else
    {
        // both segments are arcs
        EllipseKH el1;
        EllipseKH el2;
        MakeEllipseFromArc( xi, yi, xf, yf, style, &el1 );
        MakeEllipseFromArc( xi2, yi2, xf2, yf2, style2, &el2 );
        int n;
        if( el1.xrad+el1.yrad > el2.xrad+el2.yrad )
            n = GetArcIntersections( &el1, &el2 );
        else
            n = GetArcIntersections( &el2, &el1 );
        iret = n;
    }
    if( x && y )
    {
        for( int i=0; i<iret; i++ )
        {
            x[i] = xr[i];
            y[i] = yr[i];
        }
    }
    return iret;
}

// find intersection between line y = a + bx and line segment (xi,yi) to (xf,yf)
// if b > DBL_MAX/10, assume vertical line at x = a
// the line segment may be an arc (i.e. quadrant of an ellipse)
// return 0 if no intersection
// returns 1 or 2 if intersections found
// sets coords of intersections in *x1, *y1, *x2, *y2
// if no intersection, returns min distance in dist
//
int FindLineSegmentIntersection( double a, double b, int xi, int yi, int xf, int yf, int style,
                                double * x1, double * y1, double * x2, double * y2,
                                double * dist )
{
    double xx = 0, yy = 0;	//Init made to avoid C compil "uninitialized" warning
    bool bVert = false;
    if( b > DBL_MAX/10.0 )
        bVert = true;

    if( xf != xi )
    {
        // non-vertical segment, get intersection
        if( style == CPolyLine::STRAIGHT || yf == yi )
        {
            // horizontal or oblique straight segment
            // put into form y = c + dx;
            double d = (double)(yf-yi)/(double)(xf-xi);
            double c = yf - d*xf;
            if( bVert )
            {
                // if vertical line, easy
                if( InRange( a, xi, xf ) )
                {
                    *x1 = a;
                    *y1 = c + d*a;
                    return 1;
                }
                else
                {
                    if( dist )
                        *dist = min( abs(a-xi), abs(a-xf) );
                    return 0;
                }
            }
            if( fabs(b-d) < 1E-12 )
            {
                // parallel lines
                if( dist )
                {
                    *dist = GetPointToLineDistance( a, b, xi, xf );
                }
                return 0;	// lines parallel
            }
            // calculate intersection
            xx = (c-a)/(b-d);
            yy = a + b*(xx);
            // see if intersection is within the line segment
            if( yf == yi )
            {
                // horizontal line
                if( (xx>=xi && xx>xf) || (xx<=xi && xx<xf) )
                    return 0;
            }
            else
            {
                // oblique line
                if( (xx>=xi && xx>xf) || (xx<=xi && xx<xf)
                    || (yy>yi && yy>yf) || (yy<yi && yy<yf) )
                    return 0;
            }
        }
        else if( style == CPolyLine::ARC_CW || style == CPolyLine::ARC_CCW )
        {
            // arc (quadrant of ellipse)
            // convert to clockwise arc
            int xxi, xxf, yyi, yyf;
            if( style == CPolyLine::ARC_CCW )
            {
                xxi = xf;
                xxf = xi;
                yyi = yf;
                yyf = yi;
            }
            else
            {
                xxi = xi;
                xxf = xf;
                yyi = yi;
                yyf = yf;
            }
            // find center and radii of ellipse
            double xo = xxf, yo = yyi, rx, ry;		// Init made to avoid C compil warnings
            if( xxf > xxi && yyf > yyi )
            {
                xo = xxf;
                yo = yyi;
            }
            else if( xxf < xxi && yyf > yyi )
            {
                xo = xxi;
                yo = yyf;
            }
            else if( xxf < xxi && yyf < yyi )
            {
                xo = xxf;
                yo = yyi;
            }
            else if( xxf > xxi && yyf < yyi )
            {
                xo = xxi;
                yo = yyf;
            }
            rx = fabs( (double)(xxi-xxf) );
            ry = fabs( (double)(yyi-yyf) );
            bool test;
            double xx1, xx2, yy1, yy2, aa;
            if( bVert )
            {
                // shift vertical line to coordinate system of ellipse
                aa = a - xo;
                test = FindVerticalLineEllipseIntersections( rx, ry, aa, &yy1, &yy2 );
                if( !test )
                    return 0;
                // shift back to PCB coordinates
                yy1 += yo;
                yy2 += yo;
                xx1 = a;
                xx2 = a;
            }
            else
            {
                // shift line to coordinate system of ellipse
                aa = a + b*xo - yo;
                test = FindLineEllipseIntersections( rx, ry, aa, b, &xx1, &xx2 );
                if( !test )
                    return 0;
                // shift back to PCB coordinates
                yy1 = aa + b*xx1;
                xx1 += xo;
                yy1 += yo;
                yy2 = aa + b*xx2;
                xx2 += xo;
                yy2 += yo;
            }
            int npts = 0;
            if( (xxf>xxi && xx1<xxf && xx1>xxi) || (xxf<xxi && xx1<xxi && xx1>xxf) )
            {
                if( (yyf>yyi && yy1<yyf && yy1>yyi) || (yyf<yyi && yy1<yyi && yy1>yyf) )
                {
                    *x1 = xx1;
                    *y1 = yy1;
                    npts = 1;
                }
            }
            if( (xxf>xxi && xx2<xxf && xx2>xxi) || (xxf<xxi && xx2<xxi && xx2>xxf) )
            {
                if( (yyf>yyi && yy2<yyf && yy2>yyi) || (yyf<yyi && yy2<yyi && yy2>yyf) )
                {
                    if( npts == 0 )
                    {
                        *x1 = xx2;
                        *y1 = yy2;
                        npts = 1;
                    }
                    else
                    {
                        *x2 = xx2;
                        *y2 = yy2;
                        npts = 2;
                    }
                }
            }
            return npts;
        }
        else
            wxASSERT(0);
    }
    else
    {
        // vertical line segment
        if( bVert )
            return 0;
        xx = xi;
        yy = a + b*xx;
        if( (yy>=yi && yy>yf) || (yy<=yi && yy<yf) )
            return 0;
    }
    *x1 = xx;
    *y1 = yy;
    return 1;
}

/**
 * Function TestForIntersectionOfStraightLineSegments
 * Test for intersection of line segments
 * If lines are parallel, returns false
 * If true, returns also intersection coords in x, y
 * if false, returns min. distance in dist (may be 0.0 if parallel)
 */
bool TestForIntersectionOfStraightLineSegments( int x1i, int y1i, int x1f, int y1f,
                                       int x2i, int y2i, int x2f, int y2f,
                                       int * x, int * y, double * d )
{
    double a, b, dist;
    // first, test for intersection
    if( x1i == x1f && x2i == x2f )
    {
        // both segments are vertical, can't intersect
    }
    else if( y1i == y1f && y2i == y2f )
    {
        // both segments are horizontal, can't intersect
    }
    else if( x1i == x1f && y2i == y2f )
    {
        // first seg. vertical, second horizontal, see if they cross
        if( InRange( x1i, x2i, x2f )
            && InRange( y2i, y1i, y1f ) )
        {
            if( x )
                *x = x1i;
            if( y )
                *y = y2i;
            if( d )
                *d = 0.0;
            return true;
        }
    }
    else if( y1i == y1f && x2i == x2f )
    {
        // first seg. horizontal, second vertical, see if they cross
        if( InRange( y1i, y2i, y2f )
            && InRange( x2i, x1i, x1f ) )
        {
            if( x )
                *x = x2i;
            if( y )
                *y = y1i;
            if( d )
                *d = 0.0;
            return true;
        }
    }
    else if( x1i == x1f )
    {
        // first segment vertical, second oblique
        // get a and b for second line segment, so that y = a + bx;
        b = (double)(y2f-y2i)/(x2f-x2i);
        a = (double)y2i - b*x2i;
        double x1, y1, x2, y2;
        int test = FindLineSegmentIntersection( a, b, x1i, y1i, x1f, y1f, CPolyLine::STRAIGHT,
            &x1, &y1, &x2, &y2 );
        if( test )
        {
            if( InRange( y1, y1i, y1f ) && InRange( x1, x2i, x2f ) && InRange( y1, y2i, y2f ) )
            {
                if( x )
                    *x = (int) x1;
                if( y )
                    *y = (int) y1;
                if( d )
                    *d = 0.0;
                return true;
            }
        }
    }
    else if( y1i == y1f )
    {
        // first segment horizontal, second oblique
        // get a and b for second line segment, so that y = a + bx;
        b = (double)(y2f-y2i)/(x2f-x2i);
        a = (double)y2i - b*x2i;
        double x1, y1, x2, y2;
        int test = FindLineSegmentIntersection( a, b, x1i, y1i, x1f, y1f, CPolyLine::STRAIGHT,
            &x1, &y1, &x2, &y2 );
        if( test )
        {
            if( InRange( x1, x1i, x1f ) && InRange( x1, x2i, x2f ) && InRange( y1, y2i, y2f ) )
            {
                if( x )
                    *x = (int) x1;
                if( y )
                    *y = (int) y1;
                if( d )
                    *d = 0.0;
                return true;
            }
        }
    }
    else if( x2i == x2f )
    {
        // second segment vertical, first oblique
        // get a and b for first line segment, so that y = a + bx;
        b = (double)(y1f-y1i)/(x1f-x1i);
        a = (double)y1i - b*x1i;
        double x1, y1, x2, y2;
        int test = FindLineSegmentIntersection( a, b, x2i, y2i, x2f, y2f, CPolyLine::STRAIGHT,
            &x1, &y1, &x2, &y2 );
        if( test )
        {
            if( InRange( x1, x1i, x1f ) &&  InRange( y1, y1i, y1f ) && InRange( y1, y2i, y2f ) )
            {
                if( x )
                    *x = (int) x1;
                if( y )
                    *y = (int) y1;
                if( d )
                    *d = 0.0;
                return true;
            }
        }
    }
    else if( y2i == y2f )
    {
        // second segment horizontal, first oblique
        // get a and b for second line segment, so that y = a + bx;
        b = (double)(y1f-y1i)/(x1f-x1i);
        a = (double)y1i - b*x1i;
        double x1, y1, x2, y2;
        int test = FindLineSegmentIntersection( a, b, x2i, y2i, x2f, y2f, CPolyLine::STRAIGHT,
            &x1, &y1, &x2, &y2 );
        if( test )
        {
            if( InRange( x1, x1i, x1f ) && InRange( y1, y1i, y1f ) )
            {
                if( x )
                    *x = (int) x1;
                if( y )
                    *y = (int) y1;
                if( d )
                    *d = 0.0;
                return true;
            }
        }
    }
    else
    {
        // both segments oblique
        if( (long)(y1f-y1i)*(x2f-x2i) != (long)(y2f-y2i)*(x1f-x1i) )
        {
            // not parallel, get a and b for first line segment, so that y = a + bx;
            b = (double)(y1f-y1i)/(x1f-x1i);
            a = (double)y1i - b*x1i;
            double x1, y1, x2, y2;
            int test = FindLineSegmentIntersection( a, b, x2i, y2i, x2f, y2f, CPolyLine::STRAIGHT,
                &x1, &y1, &x2, &y2 );
            // both segments oblique
            if( test )
            {
                if( InRange( x1, x1i, x1f ) && InRange( y1, y1i, y1f ) )
                {
                    if( x )
                        *x = (int) x1;
                    if( y )
                        *y = (int) y1;
                    if( d )
                        *d = 0.0;
                    return true;
                }
            }
        }
    }
    // don't intersect, get shortest distance between each endpoint and the other line segment
    dist = GetPointToLineSegmentDistance( x1i, y1i, x2i, y2i, x2f, y2f );
    double xx = x1i;
    double yy = y1i;
    double dd = GetPointToLineSegmentDistance( x1f, y1f, x2i, y2i, x2f, y2f );
    if( dd < dist )
    {
        dist = dd;
        xx = x1f;
        yy = y1f;
    }
    dd = GetPointToLineSegmentDistance( x2i, y2i, x1i, y1i, x1f, y1f );
    if( dd < dist )
    {
        dist = dd;
        xx = x2i;
        yy = y2i;
    }
    dd = GetPointToLineSegmentDistance( x2f, y2f, x1i, y1i, x1f, y1f );
    if( dd < dist )
    {
        dist = dd;
        xx = x2f;
        yy = y2f;
    }
    if( x )
        *x = (int) xx;
    if( y )
        *y = (int) yy;
    if( d )
        *d = dist;
    return false;
}


/*solves the Quadratic equation = a*x*x + b*x + c
*/
bool Quadratic( double a, double b, double c, double *x1, double *x2 )
{
    double root = b*b - 4.0*a*c;
    if( root < 0.0 )
        return false;
    root = sqrt( root );
    *x1 = (-b+root)/(2.0*a);
    *x2 = (-b-root)/(2.0*a);
    return true;
}

// finds intersections of vertical line at x
// with ellipse defined by (x^2)/(a^2) + (y^2)/(b^2) = 1;
// returns true if solution exist, with solutions in y1 and y2
// else returns false
//
bool FindVerticalLineEllipseIntersections( double a, double b, double x, double *y1, double *y2 )
{
    double y_sqr = (1.0-(x*x)/(a*a))*b*b;
    if( y_sqr < 0.0 )
        return false;
    *y1 = sqrt(y_sqr);
    *y2 = -*y1;
    return true;
}

// finds intersections of straight line y = c + dx
// with ellipse defined by (x^2)/(a^2) + (y^2)/(b^2) = 1;
// returns true if solution exist, with solutions in x1 and x2
// else returns false
//
bool FindLineEllipseIntersections( double a, double b, double c, double d, double *x1, double *x2 )
{
    // quadratic terms
    double A = d*d+b*b/(a*a);
    double B = 2.0*c*d;
    double C = c*c-b*b;
    return Quadratic( A, B, C, x1, x2 );
}



// Get clearance between 2 segments
// Returns point in segment closest to other segment in x, y
// in clearance > max_cl, just returns max_cl and doesn't return x,y
//
int GetClearanceBetweenSegments( int x1i, int y1i, int x1f, int y1f, int style1, int w1,
                                   int x2i, int y2i, int x2f, int y2f, int style2, int w2,
                                   int max_cl, int * x, int * y )
{
    // check clearance between bounding rectangles
    int test = max_cl + w1/2 + w2/2;
    if( min(x1i,x1f)-max(x2i,x2f) > test )
        return max_cl;
    if( min(x2i,x2f)-max(x1i,x1f) > test )
        return max_cl;
    if( min(y1i,y1f)-max(y2i,y2f) > test )
        return max_cl;
    if( min(y2i,y2f)-max(y1i,y1f) > test )
        return max_cl;

    if( style1 == CPolyLine::STRAIGHT && style1 == CPolyLine::STRAIGHT )
    {
        // both segments are straight lines
        int xx, yy;
        double dd;
        TestForIntersectionOfStraightLineSegments( x1i, y1i, x1f, y1f,
            x2i, y2i, x2f, y2f, &xx, &yy, &dd );
        int d = max( 0, (int)dd - w1/2 - w2/2 );
        if( x )
            *x = xx;
        if( y )
            *y = yy;
        return d;
    }

    // not both straight-line segments
    // see if segments intersect
    double xr[2];
    double yr[2];
    test = FindSegmentIntersections( x1i, y1i, x1f, y1f, style1, x2i, y2i, x2f, y2f, style2, xr, yr );
    if( test )
    {
        if( x )
            *x = (int) xr[0];
        if( y )
            *y = (int) yr[0];
        return 0;
    }

    // at least one segment is an arc
    EllipseKH el1;
    EllipseKH el2;
    bool bArcs;
    int xi=0, yi=0, xf=0, yf=0;
    if( style2 == CPolyLine::STRAIGHT )
    {
        // style1 = arc, style2 = straight
        MakeEllipseFromArc( x1i, y1i, x1f, y1f, style1, &el1 );
        xi = x2i;
        yi = y2i;
        xf = x2f;
        yf = y2f;
        bArcs = false;
    }
    else if( style1 == CPolyLine::STRAIGHT )
    {
        // style2 = arc, style1 = straight
        xi = x1i;
        yi = y1i;
        xf = x1f;
        yf = y1f;
        MakeEllipseFromArc( x2i, y2i, x2f, y2f, style2, &el1 );
        bArcs = false;
    }
    else
    {
        // style1 = arc, style2 = arc
        MakeEllipseFromArc( x1i, y1i, x1f, y1f, style1, &el1 );
        MakeEllipseFromArc( x2i, y2i, x2f, y2f, style2, &el2 );
        bArcs = true;
    }
    const int NSTEPS = 32;

    if( el1.theta2 > el1.theta1 )
    {
        wxASSERT(0);
    }
    if( bArcs && el2.theta2 > el2.theta1 )
    {
        wxASSERT(0);
    }

    // test multiple points in both segments
    double th1;
    double th2;
    double len2;
    if( bArcs )
    {
        th1 = el2.theta1;
        th2 = el2.theta2;
        len2 = max(el2.xrad, el2.yrad);
    }
    else
    {
        th1 = 1.0;
        th2 = 0.0;
        len2 = abs(xf-xi)+abs(yf-yi);
    }
    double s_start = el1.theta1;
    double s_end = el1.theta2;
    double s_start2 = th1;
    double s_end2 = th2;
    double dmin = DBL_MAX;
    double xmin = 0, ymin = 0, smin = 0, smin2 = 0;		// Init made to avoid C compil warnings

    int nsteps = NSTEPS;
    int nsteps2 = NSTEPS;
    double step = (s_start-s_end)/(nsteps-1);
    double step2 = (s_start2-s_end2)/(nsteps2-1);
    while( (step * max(el1.xrad, el1.yrad)) > 0.1*NM_PER_MIL
        && (step2 * len2) > 0.1*NM_PER_MIL )
    {
        step = (s_start-s_end)/(nsteps-1);
        for( int i=0; i<nsteps; i++ )
        {
            double s;
            if( i < nsteps-1 )
                s = s_start - i*step;
            else
                s = s_end;
            double x = el1.Center.X + el1.xrad*cos(s);
            double y = el1.Center.Y + el1.yrad*sin(s);
            // if not an arc, use s2 as fractional distance along line
            step2 = (s_start2-s_end2)/(nsteps2-1);
            for( int i2=0; i2<nsteps2; i2++ )
            {
                double s2;
                if( i2 < nsteps2-1 )
                    s2 = s_start2 - i2*step2;
                else
                    s2 = s_end2;
                double x2, y2;
                if( !bArcs )
                {
                    x2 = xi + (xf-xi)*s2;
                    y2 = yi + (yf-yi)*s2;
                }
                else
                {
                    x2 = el2.Center.X + el2.xrad*cos(s2);
                    y2 = el2.Center.Y + el2.yrad*sin(s2);
                }
                double d = Distance( (int) x, (int) y, (int) x2, (int) y2 );
                if( d < dmin )
                {
                    dmin = d;
                    xmin = x;
                    ymin = y;
                    smin = s;
                    smin2 = s2;
                }
            }
        }
        if( step > step2 )
        {
            s_start = min(el1.theta1, smin + step);
            s_end = max(el1.theta2, smin - step);
            step = (s_start - s_end)/nsteps;
        }
        else
        {
            s_start2 = min(th1, smin2 + step2);
            s_end2 = max(th2, smin2 - step2);
            step2 = (s_start2 - s_end2)/nsteps2;
        }
    }
    if( x )
        *x = (int) xmin;
    if( y )
        *y = (int) ymin;
    return max(0, (int)dmin-w1/2-w2/2);	// allow for widths
}



// Get min. distance from (x,y) to line y = a + bx
// if b > DBL_MAX/10, assume vertical line at x = a
// returns closest point on line in xp, yp
//
double GetPointToLineDistance( double a, double b, int x, int y, double * xpp, double * ypp )
{
    if( b > DBL_MAX/10 )
    {
        // vertical line
        if( xpp && ypp )
        {
            *xpp = a;
            *ypp = y;
        }
        return abs(a-x);
    }
    // find c,d such that (x,y) lies on y = c + dx where d=(-1/b)
    double d = -1.0/b;
    double c = (double)y-d*x;
    // find nearest point to (x,y) on line through (xi,yi) to (xf,yf)
    double xp = (a-c)/(d-b);
    double yp = a + b*xp;
    if( xpp && ypp )
    {
        *xpp = xp;
        *ypp = yp;
    }
    // find distance
    return Distance( x, y, (int) xp, (int) yp );
}

/***********************************************************************************/
double GetPointToLineSegmentDistance( int x, int y, int xi, int yi, int xf, int yf )
/***********************************************************************************/
/**
 * Function GetPointToLineSegmentDistance
 * Get distance between line segment and point
 * @param x,y = point
 * @param	xi,yi and xf,yf = the end-points of the line segment
 * @return the distance
 */
{
    // test for vertical or horizontal segment
    if( xf==xi )
    {
        // vertical line segment
        if( InRange( y, yi, yf ) )
            return abs( x - xi );
        else
            return min( Distance( x, y, xi, yi ), Distance( x, y, xf, yf ) );
    }
    else if( yf==yi )
    {
        // horizontal line segment
        if( InRange( x, xi, xf ) )
            return abs( y - yi );
        else
            return min( Distance( x, y, xi, yi ), Distance( x, y, xf, yf ) );
    }
    else
    {
        // oblique segment
        // find a,b such that (xi,yi) and (xf,yf) lie on y = a + bx
        double b = (double)(yf-yi)/(xf-xi);
        double a = (double)yi-b*xi;
        // find c,d such that (x,y) lies on y = c + dx where d=(-1/b)
        double d = -1.0/b;
        double c = (double)y-d*x;
        // find nearest point to (x,y) on line through (xi,yi) to (xf,yf)
        double xp = (a-c)/(d-b);
        double yp = a + b*xp;
        // find distance
        if( InRange( xp, xi, xf ) && InRange( yp, yi, yf ) )
            return Distance( x, y, (int) xp, (int) yp );
        else
            return min( Distance( x, y, xi, yi ), Distance( x, y, xf, yf ) );
    }
}

// test for value within range
//
bool InRange( double x, double xi, double xf )
{
    if( xf>xi )
    {
        if( x >= xi && x <= xf )
            return true;
    }
    else
    {
        if( x >= xf && x <= xi )
            return true;
    }
    return false;
}

// Get distance between 2 points
//
double Distance( int x1, int y1, int x2, int y2 )
{
    double d;
    d = sqrt( (double)(x1-x2)*(x1-x2) + (double)(y1-y2)*(y1-y2) );
    if( d > INT_MAX || d < INT_MIN )
    {
        wxASSERT(0);
    }
    return (int)d;
}

// this finds approximate solutions
// note: this works best if el2 is smaller than el1
//
int GetArcIntersections( EllipseKH * el1, EllipseKH * el2,
                        double * x1, double * y1, double * x2, double * y2 )
{
    if( el1->theta2 > el1->theta1 )
    {
        wxASSERT(0);
    }
    if( el2->theta2 > el2->theta1 )
    {
        wxASSERT(0);
    }

    const int NSTEPS = 32;
    double xret[2], yret[2];

    double xscale = 1.0/el1->xrad;
    double yscale = 1.0/el1->yrad;
    // now transform params of second ellipse into reference frame
    // with origin at center if first ellipse,
    // scaled so the first ellipse is a circle of radius = 1.0
    double xo = (el2->Center.X - el1->Center.X)*xscale;
    double yo = (el2->Center.Y - el1->Center.Y)*yscale;
    double xr = el2->xrad*xscale;
    double yr = el2->yrad*yscale;
    // now test NSTEPS positions in arc, moving clockwise (ie. decreasing theta)
    double step = M_PI/((NSTEPS-1)*2.0);
    double d_prev=0, th_prev;
    double th_interp;
    double th1;
    int n = 0;
    for( int i=0; i<NSTEPS; i++ )
    {
        double theta;
        if( i < NSTEPS-1 )
            theta = el2->theta1 - i*step;
        else
            theta = el2->theta2;
        double x = xo + xr*cos(theta);
        double y = yo + yr*sin(theta);
        double d = 1.0 - sqrt(x*x + y*y);
        if( i>0 )
        {
            bool bInt = false;
            if( d >= 0.0 && d_prev <= 0.0 )
            {
                th_interp = theta + (step*(-d_prev))/(d-d_prev);
                bInt = true;
            }
            else if( d <= 0.0 && d_prev >= 0.0 )
            {
                th_interp = theta + (step*d_prev)/(d_prev-d);
                bInt = true;
            }
            if( bInt )
            {
                x = xo + xr*cos(th_interp);
                y = yo + yr*sin(th_interp);
                th1 = atan2( y, x );
                if( th1 <= el1->theta1 && th1 >= el1->theta2 )
                {
                    xret[n] = x*el1->xrad + el1->Center.X;
                    yret[n] = y*el1->yrad + el1->Center.Y;
                    n++;
                    if( n > 2 )
                    {
                        wxASSERT(0);
                    }
                }
            }
        }
        d_prev = d;
        th_prev = theta;
    }
    if( x1 )
        *x1 = xret[0];
    if( y1 )
        *y1 = yret[0];
    if( x2 )
        *x2 = xret[1];
    if( y2 )
        *y2 = yret[1];
    return n;
}

// this finds approximate solution
//
//double GetSegmentClearance( EllipseKH * el1, EllipseKH * el2,
double GetArcClearance( EllipseKH * el1, EllipseKH * el2,
                     double * x1, double * y1 )
{
    const int NSTEPS = 32;

    if( el1->theta2 > el1->theta1 )
    {
        wxASSERT(0);
    }
    if( el2->theta2 > el2->theta1 )
    {
        wxASSERT(0);
    }

    // test multiple positions in both arcs, moving clockwise (ie. decreasing theta)
    double th_start = el1->theta1;
    double th_end = el1->theta2;
    double th_start2 = el2->theta1;
    double th_end2 = el2->theta2;
    double dmin = DBL_MAX;
    double xmin=0, ymin=0, thmin=0, thmin2=0;

    int nsteps = NSTEPS;
    int nsteps2 = NSTEPS;
    double step = (th_start-th_end)/(nsteps-1);
    double step2 = (th_start2-th_end2)/(nsteps2-1);
    while( (step * max(el1->xrad, el1->yrad)) > 1.0*NM_PER_MIL
        && (step2 * max(el2->xrad, el2->yrad)) > 1.0*NM_PER_MIL )
    {
        step = (th_start-th_end)/(nsteps-1);
        for( int i=0; i<nsteps; i++ )
        {
            double theta;
            if( i < nsteps-1 )
                theta = th_start - i*step;
            else
                theta = th_end;
            double x = el1->Center.X + el1->xrad*cos(theta);
            double y = el1->Center.Y + el1->yrad*sin(theta);
            step2 = (th_start2-th_end2)/(nsteps2-1);
            for( int i2=0; i2<nsteps2; i2++ )
            {
                double theta2;
                if( i2 < nsteps2-1 )
                    theta2 = th_start2 - i2*step2;
                else
                    theta2 = th_end2;
                double x2 = el2->Center.X + el2->xrad*cos(theta2);
                double y2 = el2->Center.Y + el2->yrad*sin(theta2);
                double d = Distance( (int) x, (int) y, (int) x2, (int) y2 );
                if( d < dmin )
                {
                    dmin = d;
                    xmin = x;
                    ymin = y;
                    thmin = theta;
                    thmin2 = theta2;
                }
            }
        }
        if( step > step2 )
        {
            th_start = min(el1->theta1, thmin + step);
            th_end = max(el1->theta2, thmin - step);
            step = (th_start - th_end)/nsteps;
        }
        else
        {
            th_start2 = min(el2->theta1, thmin2 + step2);
            th_end2 = max(el2->theta2, thmin2 - step2);
            step2 = (th_start2 - th_end2)/nsteps2;
        }
    }
    if( x1 )
        *x1 = xmin;
    if( y1 )
        *y1 = ymin;
    return dmin;
}