github.com/richardwilkes/toolbox@v1.121.0/xmath/geom/line.go (about)

     1  // Copyright (c) 2016-2024 by Richard A. Wilkes. All rights reserved.
     2  //
     3  // This Source Code Form is subject to the terms of the Mozilla Public
     4  // License, version 2.0. If a copy of the MPL was not distributed with
     5  // this file, You can obtain one at http://mozilla.org/MPL/2.0/.
     6  //
     7  // This Source Code Form is "Incompatible With Secondary Licenses", as
     8  // defined by the Mozilla Public License, version 2.0.
     9  
    10  package geom
    11  
    12  import (
    13  	"math"
    14  
    15  	"github.com/richardwilkes/toolbox/xmath"
    16  )
    17  
    18  // LineIntersection determines the intersection of two lines, if any. A return of no points indicates no intersection.
    19  // One point indicates intersection at a single point. Two points indicates an overlapping line segment.
    20  func LineIntersection[T xmath.Numeric](a1, a2, b1, b2 Point[T]) []Point[T] {
    21  	aIsPt := a1.X == a2.X && a1.Y == a2.Y
    22  	bIsPt := b1.X == b2.X && b1.Y == b2.Y
    23  	switch {
    24  	case aIsPt && bIsPt:
    25  		if a1.X == b1.X && a1.Y == b1.Y {
    26  			return []Point[T]{a1}
    27  		}
    28  	case aIsPt:
    29  		if PointSegmentDistance(b1, b2, a1) == 0 {
    30  			return []Point[T]{a1}
    31  		}
    32  	case bIsPt:
    33  		if PointSegmentDistance(a1, a2, b1) == 0 {
    34  			return []Point[T]{b1}
    35  		}
    36  	default:
    37  		abdx := a1.X - b1.X
    38  		abdy := a1.Y - b1.Y
    39  		bdx := b2.X - b1.X
    40  		bdy := b2.Y - b1.Y
    41  		uat := bdx*abdy - bdy*abdx
    42  		adx := a2.X - a1.X
    43  		ady := a2.Y - a1.Y
    44  		ubt := adx*abdy - ady*abdx
    45  		ub := bdy*adx - bdx*ady
    46  		if ub != 0 {
    47  			// Not parallel, so find intersection point
    48  			a := uat / ub
    49  			if a >= 0 && a <= 1 {
    50  				b := ubt / ub
    51  				if b >= 0 && b <= 1 {
    52  					return []Point[T]{{X: a1.X + a*adx, Y: a1.Y + a*ady}}
    53  				}
    54  			}
    55  		} else if uat == 0 || ubt == 0 {
    56  			// Parallel, so check for overlap
    57  			var ub1, ub2 T
    58  			if xmath.Abs(adx) > xmath.Abs(ady) {
    59  				ub1 = (b1.X - a1.X) / adx
    60  				ub2 = (b2.X - a1.X) / adx
    61  			} else {
    62  				ub1 = (b1.Y - a1.Y) / ady
    63  				ub2 = (b2.Y - a1.Y) / ady
    64  			}
    65  			left := max(0, min(ub1, ub2))
    66  			right := min(1, max(ub1, ub2))
    67  			if left == right {
    68  				return []Point[T]{
    69  					{X: a2.X*left + a1.X*(1-left), Y: a2.Y*left + a1.Y*(1-left)},
    70  				}
    71  			}
    72  			return []Point[T]{
    73  				{X: a2.X*left + a1.X*(1-left), Y: a2.Y*left + a1.Y*(1-left)},
    74  				{X: a2.X*right + a1.X*(1-right), Y: a2.Y*right + a1.Y*(1-right)},
    75  			}
    76  		}
    77  	}
    78  	return nil
    79  }
    80  
    81  // PointSegmentDistance returns the distance from a point to a line segment. The distance measured is the distance
    82  // between the specified point and the closest point between the specified end points. If the specified point intersects
    83  // the line segment in between the end points, this function returns 0.
    84  func PointSegmentDistance[T xmath.Numeric](s1, s2, p Point[T]) T {
    85  	return T(math.Sqrt(float64(PointSegmentDistanceSquared(s1, s2, p))))
    86  }
    87  
    88  // PointSegmentDistanceSquared returns the square of the distance from a point to a line segment. The distance measured
    89  // is the distance between the specified point and the closest point between the specified end points. If the specified
    90  // point intersects the line segment in between the end points, this function returns 0.
    91  func PointSegmentDistanceSquared[T xmath.Numeric](s1, s2, p Point[T]) T {
    92  	vx := s2.X - s1.X
    93  	vy := s2.Y - s1.Y
    94  	px := p.X - s1.X
    95  	py := p.Y - s1.Y
    96  	dp := px*vx + py*vy
    97  	var projected T
    98  	if dp > 0 {
    99  		px = vx - px
   100  		py = vy - py
   101  		dp = px*vx + py*vy
   102  		if dp > 0 {
   103  			projected = dp * dp / (vx*vx + vy*vy)
   104  		}
   105  	}
   106  	lenSq := px*px + py*py - projected
   107  	if lenSq < 0 {
   108  		lenSq = 0
   109  	}
   110  	return lenSq
   111  }