ellipse.go

Documentation: oss.terrastruct.com/d2/lib/geo

     1  package geo
     2  
     3  import (
     4  	"math"
     5  	"sort"
     6  )
     7  
     8  type Ellipse struct {
     9  	Center *Point
    10  	Rx     float64
    11  	Ry     float64
    12  }
    13  
    14  func NewEllipse(center *Point, rx, ry float64) *Ellipse {
    15  	return &Ellipse{
    16  		Center: center,
    17  		Rx:     rx,
    18  		Ry:     ry,
    19  	}
    20  }
    21  
    22  func (e Ellipse) Intersections(segment Segment) []*Point {
    23  	// we check for intersections between the ellipse and line segment in the following manner:
    24  	// 0. we compute ignoring the ellipse's position, as if it were centered at 0,0 for a simpler computation
    25  	// 1. translate the line segment variables to match the ellipse's "translation" to 0,0
    26  	// 2. get the (infinite) line equation for the given line segment
    27  	// 3. compute the intersections between the line and ellipse
    28  	// 4. filter out intersections that are on the line but not on the line segment
    29  	// 5. apply the inverse translation to the intersection to get the result relative to where the ellipse and line actually are
    30  	intersections := []*Point{}
    31  
    32  	// ellipse equation: (x-cx)^2/rx^2 + (y-cy)^2/ry^2 = 1
    33  	// in the form: x^2/a^2 + y^2/b^2 = 1
    34  	a := e.Rx
    35  	b := e.Ry
    36  	a2 := a * a
    37  	b2 := b * b
    38  	if a <= 0 || b <= 0 {
    39  		return nil
    40  	}
    41  
    42  	// line for a line segment between (x1,y1) and (x2, y2): (see https://en.wikipedia.org/wiki/Linear_equation#Two-point_form)
    43  	//   y - y1 = ((y2 - y1)/(x2 - x1))(x - x1)
    44  	x1 := segment.Start.X - e.Center.X
    45  	y1 := segment.Start.Y - e.Center.Y
    46  	x2 := segment.End.X - e.Center.X
    47  	y2 := segment.End.Y - e.Center.Y
    48  
    49  	// Handle the edge case of a vertical line (avoiding dividing by 0)
    50  	if x1 == x2 {
    51  		// Ellipse solutions for a given x (from wolfram "(x^2)/(a^2)+(y^2)/(b^2)=1"):
    52  		// 1. y = +b*root(a^2 - x^2)/a
    53  		intersection1 := NewPoint(x1, +b*math.Sqrt(a2-x1*x1)/a)
    54  		// 2. y = -b*root(a^2 - x^2)/a
    55  		intersection2 := intersection1.Copy()
    56  		intersection2.Y *= -1
    57  
    58  		isPointOnLine := func(p *Point) bool {
    59  			ps := []float64{p.Y, y1, y2}
    60  			sort.Slice(ps, func(i, j int) bool {
    61  				return ps[i] < ps[j]
    62  			})
    63  			return ps[1] == p.Y
    64  		}
    65  
    66  		if isPointOnLine(intersection1) {
    67  			intersections = append(intersections, intersection1)
    68  		}
    69  		// when y = 0, intersection2 will be a duplicate of intersection1
    70  		if intersection2.Y != 0.0 && isPointOnLine(intersection2) {
    71  			intersections = append(intersections, intersection2)
    72  		}
    73  
    74  		for _, p := range intersections {
    75  			p.X += e.Center.X
    76  			p.Y += e.Center.Y
    77  		}
    78  		return intersections
    79  	}
    80  
    81  	// converting line to form: y = mx + c
    82  	// from: y - y1 = ((y2-y1)/(x2-x1))(x - x1)
    83  	// m = (y2-y1)/(x2-x1), c = y1 - m*x1
    84  	m := (y2 - y1) / (x2 - x1)
    85  	c := y1 - m*x1
    86  	isPointOnLine := func(p *Point) bool {
    87  		return PrecisionCompare(p.DistanceToLine(NewPoint(x1, y1), NewPoint(x2, y2)), 0, PRECISION) == 0
    88  	}
    89  
    90  	// "(x^2)/(a^2)+(y^2)/(b^2) =1, y = mx + c" solutions from wolfram:
    91  	// if (a^2)(m^2) + b^2 != 0, and ab != 0
    92  	// 2 solutions 1 with +, 1 with -
    93  	// x = ( -c m a^2 {+/-} root(a^2 b^2 (a^2 m^2 + b^2 - c^2)) ) / (a^2 m^2 + b^2)
    94  	// y = ( c b^2 {+/-} m root(a^2 b^2 (a^2 m^2 + b^2 - c^2)) ) / (a^2 m^2 + b^2)
    95  	denom := a2*m*m + b2
    96  	// Note: we already checked a and b != 0 so denom == 0 is impossible (assuming no imaginary numbers)
    97  	root := math.Sqrt(a2 * b2 * (denom - c*c))
    98  
    99  	intersection1 := NewPoint((-m*c*a2+root)/denom, (c*b2+m*root)/denom)
   100  	intersection2 := NewPoint((-m*c*a2-root)/denom, (c*b2-m*root)/denom)
   101  
   102  	if isPointOnLine(intersection1) {
   103  		intersections = append(intersections, intersection1)
   104  	}
   105  	if !intersection1.Equals(intersection2) && isPointOnLine(intersection2) {
   106  		intersections = append(intersections, intersection2)
   107  	}
   108  
   109  	for _, p := range intersections {
   110  		p.X += e.Center.X
   111  		p.Y += e.Center.Y
   112  	}
   113  	return intersections
   114  }
   115
View as plain text