seg.cpp

Go to the documentation of this file.

/*
 * This program source code file is part of KiCad, a free EDA CAD application.
 *
 * Copyright (C) 2013 CERN
 * @author Tomasz Wlostowski <[email protected]>
 * Copyright The KiCad Developers, see AUTHORS.txt for contributors.
 *
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License
 * as published by the Free Software Foundation; either version 2
 * of the License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <https://www.gnu.org/licenses/>.
 */
 
#include <algorithm>        // for min
#include <geometry/seg.h>
#include <math/util.h>      // for rescale
#include <math/vector2d.h>  // for VECTOR2I, VECTOR2
#include <trigo.h>          // for RAD2DEG
 
template <typename T>
int sgn( T aVal )
{
    return ( T( 0 ) < aVal ) - ( aVal < T( 0 ) );
}
 
template <typename T>
constexpr T sqrt_helper(T x, T lo, T hi)
{
  if (lo == hi)
    return lo;
 
  const T mid = (lo + hi + 1) / 2;
  if (x / mid < mid)
    return sqrt_helper<T>(x, lo, mid - 1);
  else
    return sqrt_helper(x, mid, hi);
}
 
template <typename T>
constexpr T ct_sqrt(T x)
{
  return sqrt_helper<T>(x, 0, x / 2 + 1);
}
 
template <typename T>
static constexpr T sqrt_max_typed = ct_sqrt( std::numeric_limits<T>::max() );
 
template <typename T>
T isqrt(T x)
{
  T sqrt_max = sqrt_max_typed<T>;
 
  if( x < 0 )
    return sqrt_max;
 
  T r = (T) std::sqrt( (double) x );
 
  while( r < sqrt_max && r * r < x )
    r++;
 
  while( r > sqrt_max || r * r > x )
    r--;
 
  return r;
}
 
 
SEG::ecoord SEG::SquaredDistance( const SEG& aSeg ) const
{
    // Handle zero-length segments (points) specially.
    // The Intersects() check below doesn't handle this case correctly because
    // the cross product with a zero vector is always zero, causing false positives.
    if( A == B )
        return aSeg.SquaredDistance( A );
 
    if( aSeg.A == aSeg.B )
        return SquaredDistance( aSeg.A );
 
    if( Intersects( aSeg ) )
        return 0;
 
    const VECTOR2I pts[4] =
    {
        aSeg.NearestPoint( A ) - A,
        aSeg.NearestPoint( B ) - B,
        NearestPoint( aSeg.A ) - aSeg.A,
        NearestPoint( aSeg.B ) - aSeg.B
    };
 
    ecoord m = VECTOR2I::ECOORD_MAX;
 
    for( int i = 0; i < 4; i++ )
        m = std::min( m, pts[i].SquaredEuclideanNorm() );
 
    return m;
}
 
 
EDA_ANGLE SEG::Angle( const SEG& aOther ) const
{
    EDA_ANGLE thisAngle = EDA_ANGLE( A - B ).Normalize180();
    EDA_ANGLE otherAngle = EDA_ANGLE( aOther.A - aOther.B ).Normalize180();
 
    return std::abs( ( thisAngle - otherAngle ).Normalize180() );
}
 
 
const VECTOR2I SEG::NearestPoint( const SEG& aSeg ) const
{
    if( OPT_VECTOR2I p = Intersect( aSeg ) )
        return *p;
 
    const VECTOR2I pts_origin[4] =
    {
        aSeg.NearestPoint( A ),
        aSeg.NearestPoint( B ),
        NearestPoint( aSeg.A ),
        NearestPoint( aSeg.B )
    };
 
 
    const VECTOR2I* pts_out[4] =
    {
        &A,
        &B,
        &pts_origin[2],
        &pts_origin[3]
    };
 
    const ecoord pts_dist[4] =
    {
        ( pts_origin[0] - A ).SquaredEuclideanNorm(),
        ( pts_origin[1] - B ).SquaredEuclideanNorm(),
        ( pts_origin[2] - aSeg.A ).SquaredEuclideanNorm(),
        ( pts_origin[3] - aSeg.B ).SquaredEuclideanNorm()
    };
 
    int min_i = 0;
 
    for( int i = 0; i < 4; i++ )
    {
        if( pts_dist[i] < pts_dist[min_i] )
            min_i = i;
    }
 
    return *pts_out[min_i];
}
 
 
bool SEG::NearestPoints( const SEG& aSeg, VECTOR2I& aPtA, VECTOR2I& aPtB, int64_t& aDistSq ) const
{
    if( OPT_VECTOR2I p = Intersect( aSeg ) )
    {
        aPtA = aPtB = *p;
        aDistSq = 0;
 
        return true;
    }
 
    const VECTOR2I pts_origin[4] =
    {
        aSeg.NearestPoint( A ),
        aSeg.NearestPoint( B ),
        NearestPoint( aSeg.A ),
        NearestPoint( aSeg.B )
    };
 
    const VECTOR2I* pts_a_out[4] =
    {
        &A,
        &B,
        &pts_origin[2],
        &pts_origin[3]
    };
 
    const VECTOR2I* pts_b_out[4] =
    {
        &pts_origin[0],
        &pts_origin[1],
        &aSeg.A,
        &aSeg.B
    };
 
    const ecoord pts_dist[4] =
    {
        ( pts_origin[0] - A ).SquaredEuclideanNorm(),
        ( pts_origin[1] - B ).SquaredEuclideanNorm(),
        ( pts_origin[2] - aSeg.A ).SquaredEuclideanNorm(),
        ( pts_origin[3] - aSeg.B ).SquaredEuclideanNorm()
    };
 
    int min_i = 0;
 
    for( int i = 0; i < 4; i++ )
    {
        if( pts_dist[i] < pts_dist[min_i] )
            min_i = i;
    }
 
    aPtA = *pts_a_out[min_i];
    aPtB = *pts_b_out[min_i];
    aDistSq = pts_dist[min_i];
 
    return true;
}
 
 
bool SEG::checkCollinearOverlap( const SEG& aSeg, bool useXAxis, bool aIgnoreEndpoints, VECTOR2I* aPt ) const
{
    // Extract coordinates based on the chosen axis
    int seg1_start, seg1_end, seg2_start, seg2_end;
    int coord1_start, coord1_end; // For calculating other axis coordinate
 
    if( useXAxis )
    {
        seg1_start = A.x;       seg1_end = B.x;
        seg2_start = aSeg.A.x;  seg2_end = aSeg.B.x;
        coord1_start = A.y;     coord1_end = B.y;
    }
    else
    {
        seg1_start = A.y;       seg1_end = B.y;
        seg2_start = aSeg.A.y;  seg2_end = aSeg.B.y;
        coord1_start = A.x;     coord1_end = B.x;
    }
 
    // Find segment ranges on the projection axis
    const int seg1_min = std::min( seg1_start, seg1_end );
    const int seg1_max = std::max( seg1_start, seg1_end );
    const int seg2_min = std::min( seg2_start, seg2_end );
    const int seg2_max = std::max( seg2_start, seg2_end );
 
    // Check for overlap
    const bool overlaps = seg1_max >= seg2_min && seg2_max >= seg1_min;
    if( !overlaps )
        return false;
 
    // Check if intersection is only at endpoints when aIgnoreEndpoints is true
    if( aIgnoreEndpoints )
    {
        // Calculate overlap region
        const int overlap_start = std::max( seg1_min, seg2_min );
        const int overlap_end = std::min( seg1_max, seg2_max );
 
        // If overlap region has zero length, segments only touch at endpoint
        if( overlap_start == overlap_end )
        {
            // Check if this endpoint touching involves actual segment endpoints
            // (not just projected endpoints due to min/max calculation)
            bool isEndpointTouch = false;
 
            // Check if the touch point corresponds to actual segment endpoints
            if( overlap_start == seg1_min || overlap_start == seg1_max )
            {
                // Touch point is at seg1's endpoint, check if it's also at seg2's endpoint
                if( overlap_start == seg2_min || overlap_start == seg2_max )
                {
                    isEndpointTouch = true;
                }
            }
 
            if( isEndpointTouch )
                return false;  // Ignore endpoint-only intersection
        }
    }
 
    // Calculate intersection point if requested
    if( aPt )
    {
        // Find midpoint of overlap region
        const int overlap_start = std::max( seg1_min, seg2_min );
        const int overlap_end = std::min( seg1_max, seg2_max );
        const int intersection_proj = ( overlap_start + overlap_end ) / 2;
 
        // Calculate corresponding coordinate on the other axis
        int intersection_other;
        if( seg1_end != seg1_start )
        {
            // Use this segment's line equation to find other coordinate
            intersection_other = coord1_start + static_cast<int>(
                rescale( intersection_proj - seg1_start, coord1_end - coord1_start, seg1_end - seg1_start ) );
        }
        else
        {
            // Degenerate segment (point) or perpendicular to projection axis
            intersection_other = coord1_start;
        }
 
        // Set result based on projection axis
        if( useXAxis )
            *aPt = VECTOR2I( intersection_proj, intersection_other );
        else
            *aPt = VECTOR2I( intersection_other, intersection_proj );
    }
 
    return true;
}
 
 
bool SEG::intersects( const SEG& aSeg, bool aIgnoreEndpoints, bool aLines, VECTOR2I* aPt ) const
{
    // Quick rejection: check if segment bounding boxes overlap
    // (Skip for line mode since infinite lines can intersect anywhere)
    if( !aLines )
    {
        const int this_min_x = std::min( A.x, B.x );
        const int this_max_x = std::max( A.x, B.x );
        const int this_min_y = std::min( A.y, B.y );
        const int this_max_y = std::max( A.y, B.y );
 
        const int other_min_x = std::min( aSeg.A.x, aSeg.B.x );
        const int other_max_x = std::max( aSeg.A.x, aSeg.B.x );
        const int other_min_y = std::min( aSeg.A.y, aSeg.B.y );
        const int other_max_y = std::max( aSeg.A.y, aSeg.B.y );
 
        if( this_max_x < other_min_x || other_max_x < this_min_x ||
            this_max_y < other_min_y || other_max_y < this_min_y )
        {
            return false;
        }
    }
 
    // Calculate direction vectors and offset vector using VECTOR2 operations
    // Using parametric form: P₁ = A + t*dir1, P₂ = aSeg.A + s*dir2
    const VECTOR2L dir1 = VECTOR2L( B ) - A;           // direction vector e
    const VECTOR2L dir2 = VECTOR2L( aSeg.B ) - aSeg.A; // direction vector f
    const VECTOR2L offset = VECTOR2L( aSeg.A ) - A;    // offset vector ac
    const ecoord determinant = dir2.Cross( dir1 );
 
    // Handle parallel/collinear case
    if( determinant == 0 )
    {
        // Check if lines are collinear (not just parallel) using cross product
        // Lines are collinear if offset vector is also parallel to direction vector
        const ecoord collinear_test = dir1.Cross( offset );
 
        if( collinear_test != 0 )
            return false;  // Parallel but not collinear
 
        // Lines are collinear - for infinite lines, they always intersect
        if( aLines )
        {
            // For infinite collinear lines, intersection point is ambiguous
            // Use the midpoint between the two segment start points as a reasonable choice
            if( aPt )
            {
                // If aSeg is degenerate (point), use its start point
                if( aSeg.A == aSeg.B )
                {
                    *aPt = aSeg.A;
                }
                else if( A == B )
                { // If this segment is degenerate (point), use its start point
                    *aPt = A;
                }
                else
                {
                    const VECTOR2I midpoint = ( A + aSeg.A ) / 2;
                    *aPt = midpoint;
                }
            }
            return true;
        }
 
        // For segments, check overlap using the axis with larger coordinate range
        const bool use_x_axis = std::abs( dir1.x ) >= std::abs( dir1.y );
        return checkCollinearOverlap( aSeg, use_x_axis, aIgnoreEndpoints, aPt );
    }
 
    // param2_num = f × ac (parameter for second segment: s = p/d)
    // param1_num = e × ac (parameter for first segment: t = q/d)
    const ecoord param2_num = dir2.Cross( offset );
    const ecoord param1_num = dir1.Cross( offset );
 
    // For segments (not infinite lines), check if intersection is within both segments
    if( !aLines )
    {
        // Parameters must be in [0,1] for intersection within segments
        // Since we're comparing t = q/d and s = p/d to [0,1], we need to handle sign of d
        if( determinant > 0 )
        {
            // d > 0: check 0 ≤ q ≤ d and 0 ≤ p ≤ d
            if( param1_num < 0 || param1_num > determinant ||
                param2_num < 0 || param2_num > determinant )
                return false;
        }
        else
        {
            // d < 0: check d ≤ q ≤ 0 and d ≤ p ≤ 0
            if( param1_num > 0 || param1_num < determinant ||
                param2_num > 0 || param2_num < determinant )
                return false;
        }
 
        // Optionally exclude endpoint intersections (when segments share vertices)
        if( aIgnoreEndpoints &&
            ( param1_num == 0 || param1_num == determinant ) &&
            ( param2_num == 0 || param2_num == determinant ) )
        {
            return false;
        }
    }
 
    if( aPt )
    {
        // Use parametric equation: intersection = aSeg.A + (q/d) * f
        const VECTOR2L scaled_dir2( rescale( param1_num, dir2.x, determinant ),
                                   rescale( param1_num, dir2.y, determinant ) );
        const VECTOR2L result = VECTOR2L( aSeg.A ) + scaled_dir2;
 
        // Verify result fits in coordinate type range
        constexpr ecoord max_coord = std::numeric_limits<VECTOR2I::coord_type>::max();
        constexpr ecoord min_coord = std::numeric_limits<VECTOR2I::coord_type>::min();
 
        if( result.x > max_coord || result.x < min_coord ||
            result.y > max_coord || result.y < min_coord )
        {
            return false;  // Intersection exists but coordinates overflow
        }
 
        *aPt = VECTOR2I( static_cast<int>( result.x ), static_cast<int>( result.y ) );
    }
 
    return true;
}
 
 
bool SEG::Intersects( const SEG& aSeg ) const
{
    return intersects( aSeg );
}
 
 
OPT_VECTOR2I SEG::Intersect( const SEG& aSeg, bool aIgnoreEndpoints, bool aLines ) const
{
    VECTOR2I ip;
 
    if( intersects( aSeg, aIgnoreEndpoints, aLines, &ip ) )
        return ip;
    else
        return OPT_VECTOR2I();
}
 
 
bool SEG::IntersectsLine( double aSlope, double aOffset, VECTOR2I& aIntersection ) const
{
    const VECTOR2L segA( A );
    const VECTOR2L segB( B );
    const VECTOR2L segDir = segB - segA;
 
    // Handle vertical segment case
    if( segDir.x == 0 )
    {
        // Vertical segment: x = A.x, find y on the line
        const double intersect_y = aSlope * A.x + aOffset;
 
        // Compare before rounding.  The line is unbounded, so a miss lands anywhere, and
        // narrowing it first would fault on a value this is about to discard
        const int seg_min_y = std::min( A.y, B.y );
        const int seg_max_y = std::max( A.y, B.y );
 
        if( intersect_y >= seg_min_y && intersect_y <= seg_max_y )
        {
            aIntersection = VECTOR2I( A.x, KiROUND( intersect_y ) );
            return true;
        }
        return false;
    }
 
    const VECTOR2L lineDir( 1000, static_cast<ecoord>( aSlope * 1000 ) );
    const ecoord cross_product = segDir.Cross( lineDir );
 
    if( cross_product == 0 )
    {
        // Parallel lines - check if segment lies on the line
        const double expected_y = aSlope * A.x + aOffset;
        const double diff = std::abs( A.y - expected_y );
 
        if( diff < 0.5 )
        {
            // Collinear: segment lies on the line, return midpoint
            aIntersection = ( A + B ) / 2;
            return true;
        }
 
        return false; // Parallel but not collinear
    }
 
    // Find intersection using parametric equations
    // Segment: P = segA + t * segDir
    // Line: y = aSlope * x + aOffset
    //
    // At intersection: segA.y + t * segDir.y = aSlope * (segA.x + t * segDir.x) + aOffset
    // Solving for t: t = (aSlope * segA.x + aOffset - segA.y) / (segDir.y - aSlope * segDir.x)
 
    const double numerator = aSlope * segA.x + aOffset - segA.y;
    const double denominator = segDir.y - aSlope * segDir.x;
 
    const double t = numerator / denominator;
 
    // Check if intersection is within segment bounds
    if( t >= 0.0 && t <= 1.0 )
    {
        aIntersection = KiROUND( segA.x + t * segDir.x, segA.y + t * segDir.y );
        return true;
    }
 
    return false;
}
 
 
SEG SEG::PerpendicularSeg( const VECTOR2I& aP ) const
{
    VECTOR2I slope( B - A );
    VECTOR2I endPoint = slope.Perpendicular() + aP;
 
    return SEG( aP, endPoint );
}
 
 
SEG SEG::ParallelSeg( const VECTOR2I& aP ) const
{
    VECTOR2I slope( B - A );
    VECTOR2I endPoint = slope + aP;
 
    return SEG( aP, endPoint );
}
 
 
bool SEG::Collide( const SEG& aSeg, int aClearance, int* aActual ) const
{
    // Handle negative clearance
    if( aClearance < 0 )
    {
        if( aActual )
            *aActual = 0;
 
        return false;
    }
 
    // Handle zero-length segments (points) specially.
    // The intersects() check below doesn't handle this case correctly because
    // the cross product with a zero vector is always zero, causing false positives.
    if( A == B )
    {
        int dist = aSeg.Distance( A );
 
        if( aActual )
            *aActual = dist;
 
        return dist == 0 || dist < aClearance;
    }
 
    if( aSeg.A == aSeg.B )
    {
        int dist = Distance( aSeg.A );
 
        if( aActual )
            *aActual = dist;
 
        return dist == 0 || dist < aClearance;
    }
 
    // Check for exact intersection first
    if( intersects( aSeg, false, false ) )
    {
        if( aActual )
            *aActual = 0;
 
        return true;
    }
 
    const ecoord clearance_sq = static_cast<ecoord>( aClearance ) * aClearance;
    ecoord min_dist_sq = VECTOR2I::ECOORD_MAX;
 
    auto checkDistance = [&]( ecoord dist, ecoord& min_dist ) -> bool
    {
        if( dist == 0 )
        {
            if( aActual )
                *aActual = 0;
 
            return true;
        }
 
        min_dist = std::min( min_dist, dist );
        return false; // Continue checking
    };
 
    // There are 4 points to check: start and end of this segment, and
    // start and end of the other segment.
    if( checkDistance( SquaredDistance( aSeg.A ), min_dist_sq ) ||
        checkDistance( SquaredDistance( aSeg.B ), min_dist_sq ) ||
        checkDistance( aSeg.SquaredDistance( A ), min_dist_sq ) ||
        checkDistance( aSeg.SquaredDistance( B ), min_dist_sq ) )
    {
        return true;
    }
 
    if( min_dist_sq < clearance_sq )
    {
        if( aActual )
            *aActual = static_cast<int>( isqrt( min_dist_sq ) );
 
        return true;
    }
 
    if( aActual )
        *aActual = static_cast<int>( isqrt( min_dist_sq ) );
 
    return false;
}
 
 
bool SEG::Contains( const VECTOR2I& aP, int aSqDistanceThreshold ) const
{
    /* Warning. This code does work in most - but not all of the cases.
       Take the segment (0,0) - (10,0) and the point (9, 1). With the distance threshold of 3,
       the point is assumed as contained in by the segment:
 
             (9,1)\  
       (0,0) ------ (10, 0)
 
       I've made the distance threshold configurable (using the default of value of 3),
       to not break any existing code that relies on the current rounding behaviour, but we need to
       think of more accurate (and degeneracy-free) solution. */
    
    return SquaredDistance( aP ) <= aSqDistanceThreshold;
}
 
 
const VECTOR2I SEG::NearestPoint( const VECTOR2I& aP ) const
{
    // Inlined for performance reasons
    VECTOR2L d;
    d.x = static_cast<int64_t>( B.x ) - A.x;
    d.y = static_cast<int64_t>( B.y ) - A.y;
 
    ecoord   l_squared( d.x * d.x + d.y * d.y );
 
    if( l_squared == 0 )
        return A;
 
    // Inlined for performance reasons
    VECTOR2L pa;
    pa.x = static_cast<int64_t>( aP.x ) - A.x;
    pa.y = static_cast<int64_t>( aP.y ) - A.y;
 
    ecoord t = d.Dot( pa );
 
    if( t < 0 )
        return A;
    else if( t > l_squared )
        return B;
 
    ecoord xp = rescale( t, (ecoord) d.x, l_squared );
    ecoord yp = rescale( t, (ecoord) d.y, l_squared );
 
    return VECTOR2<ecoord>( A.x + xp, A.y + yp );
}
 
 
const VECTOR2I SEG::ReflectPoint( const VECTOR2I& aP ) const
{
    VECTOR2I                d = B - A;
    VECTOR2I::extended_type l_squared = d.Dot( d );
    VECTOR2I::extended_type t = d.Dot( aP - A );
    VECTOR2<ecoord>         c;
 
    if( !l_squared )
    {
        c = aP;
    }
    else
    {
        c.x = A.x + rescale( t, static_cast<VECTOR2I::extended_type>( d.x ), l_squared );
        c.y = A.y + rescale( t, static_cast<VECTOR2I::extended_type>( d.y ), l_squared );
    }
 
    return VECTOR2<ecoord>( 2 * c.x - aP.x, 2 * c.y - aP.y );
}
 
 
VECTOR2I SEG::LineProject( const VECTOR2I& aP ) const
{
    VECTOR2I d = B - A;
    ecoord   l_squared = d.Dot( d );
 
    if( l_squared == 0 )
        return A;
 
    ecoord t = d.Dot( aP - A );
 
    ecoord xp = rescale( t, ecoord{ d.x }, l_squared );
    ecoord yp = rescale( t, ecoord{ d.y }, l_squared );
 
    return VECTOR2<ecoord>( A.x + xp, A.y + yp );
}
 
 
int SEG::Distance( const SEG& aSeg ) const
{
    return int( isqrt( SquaredDistance( aSeg ) ) );
}
 
 
int SEG::Distance( const VECTOR2I& aP ) const
{
    return int( isqrt( SquaredDistance( aP ) ) );
}
 
 
SEG::ecoord SEG::SquaredDistance( const VECTOR2I& aP ) const
{
    VECTOR2<ecoord> ab( ecoord( B.x ) - A.x, ecoord( B.y ) - A.y );
    VECTOR2<ecoord> ap( ecoord( aP.x ) - A.x, ecoord( aP.y ) - A.y );
 
    ecoord e = ap.Dot( ab );
 
    if( e <= 0 )
        return ap.SquaredEuclideanNorm();
 
    ecoord f = ab.SquaredEuclideanNorm();
 
    if( e >= f )
    {
        VECTOR2<ecoord> bp( ecoord( aP.x ) - B.x, ecoord( aP.y ) - B.y );
 
        return bp.Dot( bp );
    }
 
    const double g = ap.SquaredEuclideanNorm() - ( double( e ) * e ) / f;
 
    // The only way g can be negative is if there was a rounding error since
    // e is the projection of aP onto ab and therefore cannot be greater than
    // the length of ap and f is guaranteed to be greater than e, meaning
    // e * e / f cannot be greater than ap.SquaredEuclideanNorm()
    if( g < 0 || g > static_cast<double>( std::numeric_limits<ecoord>::max() ) )
        return 0;
 
    return KiROUND<double, ecoord>( g );
}
 
 
int SEG::LineDistance( const VECTOR2I& aP, bool aDetermineSide ) const
{
    ecoord p = ecoord{ A.y } - B.y;
    ecoord q = ecoord{ B.x } - A.x;
    ecoord r = -p * A.x - q * A.y;
    ecoord l = p * p + q * q;
    ecoord det = p * aP.x + q * aP.y + r;
    ecoord dist_sq = 0;
 
    if( l > 0 )
    {
        dist_sq  = rescale( det, det, l );
    }
 
    ecoord dist = isqrt( dist_sq );
 
    return static_cast<int>( aDetermineSide ? sgn( det ) * dist : std::abs( dist ) );
}
 
 
bool SEG::mutualDistanceSquared( const SEG& aSeg, ecoord& aD1, ecoord& aD2 ) const
{
    SEG a( *this );
    SEG b( aSeg );
 
    if( a.SquaredLength() < b.SquaredLength() )
        std::swap(a, b);
 
    ecoord p = ecoord{ a.A.y } - a.B.y;
    ecoord q = ecoord{ a.B.x } - a.A.x;
    ecoord r = -p * a.A.x - q * a.A.y;
 
    ecoord l = p * p + q * q;
 
    if( l == 0 )
        return false;
 
    ecoord det1 = p * b.A.x + q * b.A.y + r;
    ecoord det2 = p * b.B.x + q * b.B.y + r;
 
    ecoord dsq1 = rescale( det1, det1, l );
    ecoord dsq2 = rescale( det2, det2, l );
 
    aD1 = sgn( det1 ) * dsq1;
    aD2 = sgn( det2 ) * dsq2;
 
    return true;
}
 
bool SEG::ApproxCollinear( const SEG& aSeg, int aDistanceThreshold ) const
{
    ecoord thresholdSquared = Square( aDistanceThreshold );
    ecoord d1_sq, d2_sq;
 
    if( !mutualDistanceSquared( aSeg, d1_sq, d2_sq ) )
        return false;
 
    return std::abs( d1_sq ) <= thresholdSquared && std::abs( d2_sq ) <= thresholdSquared;
}
 
 
bool SEG::ApproxParallel( const SEG& aSeg, int aDistanceThreshold ) const
{
    ecoord thresholdSquared = Square( aDistanceThreshold );
    ecoord d1_sq, d2_sq;
 
    if( ! mutualDistanceSquared( aSeg, d1_sq, d2_sq ) )
        return false;
 
    return std::abs( d1_sq - d2_sq ) <= thresholdSquared;
}
 
 
bool SEG::ApproxPerpendicular( const SEG& aSeg ) const
{
    SEG perp = PerpendicularSeg( A );
 
    return aSeg.ApproxParallel( perp );
}