// DistanceFromPointToLine computes the distance from a point p to a line segment lineStart/lineEnd
//
// Note: NON-ROBUST!
341313
func DistanceFromPointToLine(p, lineStart, lineEnd geom.Coord) float64 {
	// if start = end, then just compute distance to one of the endpoints
	if lineStart[0] == lineEnd[0] && lineStart[1] == lineEnd[1] {
		return internal.Distance2D(p, lineStart)
	}

	// otherwise use comp.graphics.algorithms Frequently Asked Questions method

	// (1) r = AC dot AB
	//         ---------
	//         ||AB||^2
	//
	// r has the following meaning:
	//   r=0 P = A
	//   r=1 P = B
	//   r<0 P is on the backward extension of AB
	//   r>1 P is on the forward extension of AB
	//   0<r<1 P is interior to AB

	len2 := (lineEnd[0]-lineStart[0])*(lineEnd[0]-lineStart[0]) + (lineEnd[1]-lineStart[1])*(lineEnd[1]-lineStart[1])
	r := ((p[0]-lineStart[0])*(lineEnd[0]-lineStart[0]) + (p[1]-lineStart[1])*(lineEnd[1]-lineStart[1])) / len2

	if r <= 0.0 {
		return internal.Distance2D(p, lineStart)
	}
	if r >= 1.0 {
		return internal.Distance2D(p, lineEnd)
	}

	// (2) s = (Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay)
	//         -----------------------------
	//                    L^2
	//
	// Then the distance from C to P = |s|*L.
	//
	// This is the same calculation as {@link #distancePointLinePerpendicular}.
	// Unrolled here for performance.
	s := ((lineStart[1]-p[1])*(lineEnd[0]-lineStart[0]) - (lineStart[0]-p[0])*(lineEnd[1]-lineStart[1])) / len2
	return math.Abs(s) * math.Sqrt(len2)
}
1013275
func (l *Line) Midpoint() *Point {
	return &Point{(l.a[0] + l.b[0]) / 2, (l.a[1] + l.b[1]) / 2}
}
96641
func (l *Loop) iteratorContainsPoint(it *ShapeIndexIterator, p Point) bool {
	// Test containment by drawing a line segment from the cell center to the
	// given point and counting edge crossings.
	aClipped := it.IndexCell().findByShapeID(0)
	inside := aClipped.containsCenter
	if len(aClipped.edges) > 0 {
		center := it.Center()
		crosser := NewEdgeCrosser(center, p)
		aiPrev := -2
		for _, ai := range aClipped.edges {
			if ai != aiPrev+1 {
				crosser.RestartAt(l.Vertex(ai))
			}
			aiPrev = ai
			inside = inside != crosser.EdgeOrVertexChainCrossing(l.Vertex(ai+1))
		}
	}
	return inside
}
22015
func (p *Point) PointDistanceSquared(p2 Point) float32 {
	return (p.X-p2.X)*(p.X-p2.X) + (p.Y-p2.Y)*(p.Y-p2.Y)
}
136507
func (r *RectBounder) AddPoint(b Point) {
	bLL := LatLngFromPoint(b)

	if r.bound.IsEmpty() {
		r.a = b
		r.aLL = bLL
		r.bound = r.bound.AddPoint(bLL)
		return
	}

	// First compute the cross product N = A x B robustly. This is the normal
	// to the great circle through A and B. We don't use RobustSign
	// since that method returns an arbitrary vector orthogonal to A if the two
	// vectors are proportional, and we want the zero vector in that case.
	n := r.a.Sub(b.Vector).Cross(r.a.Add(b.Vector)) // N = 2 * (A x B)

	// The relative error in N gets large as its norm gets very small (i.e.,
	// when the two points are nearly identical or antipodal). We handle this
	// by choosing a maximum allowable error, and if the error is greater than
	// this we fall back to a different technique. Since it turns out that
	// the other sources of error in converting the normal to a maximum
	// latitude add up to at most 1.16 * dblEpsilon, and it is desirable to
	// have the total error be a multiple of dblEpsilon, we have chosen to
	// limit the maximum error in the normal to be 3.84 * dblEpsilon.
	// It is possible to show that the error is less than this when
	//
	// n.Norm() >= 8 * sqrt(3) / (3.84 - 0.5 - sqrt(3)) * dblEpsilon
	//          = 1.91346e-15 (about 8.618 * dblEpsilon)
	nNorm := n.Norm()
	if nNorm < 1.91346e-15 {
		// A and B are either nearly identical or nearly antipodal (to within
		// 4.309 * dblEpsilon, or about 6 nanometers on the earth's surface).
		if r.a.Dot(b.Vector) < 0 {
			// The two points are nearly antipodal. The easiest solution is to
			// assume that the edge between A and B could go in any direction
			// around the sphere.
			r.bound = FullRect()
		} else {
			// The two points are nearly identical (to within 4.309 * dblEpsilon).
			// In this case we can just use the bounding rectangle of the points,
			// since after the expansion done by GetBound this Rect is
			// guaranteed to include the (lat,lng) values of all points along AB.
			r.bound = r.bound.Union(RectFromLatLng(r.aLL).AddPoint(bLL))
		}
		r.a = b
		r.aLL = bLL
		return
	}

	// Compute the longitude range spanned by AB.
	lngAB := s1.EmptyInterval().AddPoint(r.aLL.Lng.Radians()).AddPoint(bLL.Lng.Radians())
	if lngAB.Length() >= math.Pi-2*dblEpsilon {
		// The points lie on nearly opposite lines of longitude to within the
		// maximum error of the calculation. The easiest solution is to assume
		// that AB could go on either side of the pole.
		lngAB = s1.FullInterval()
	}

	// Next we compute the latitude range spanned by the edge AB. We start
	// with the range spanning the two endpoints of the edge:
	latAB := r1.IntervalFromPoint(r.aLL.Lat.Radians()).AddPoint(bLL.Lat.Radians())

	// This is the desired range unless the edge AB crosses the plane
	// through N and the Z-axis (which is where the great circle through A
	// and B attains its minimum and maximum latitudes). To test whether AB
	// crosses this plane, we compute a vector M perpendicular to this
	// plane and then project A and B onto it.
	m := n.Cross(r3.Vector{0, 0, 1})
	mA := m.Dot(r.a.Vector)
	mB := m.Dot(b.Vector)

	// We want to test the signs of "mA" and "mB", so we need to bound
	// the error in these calculations. It is possible to show that the
	// total error is bounded by
	//
	// (1 + sqrt(3)) * dblEpsilon * nNorm + 8 * sqrt(3) * (dblEpsilon**2)
	//   = 6.06638e-16 * nNorm + 6.83174e-31

	mError := 6.06638e-16*nNorm + 6.83174e-31
	if mA*mB < 0 || math.Abs(mA) <= mError || math.Abs(mB) <= mError {
		// Minimum/maximum latitude *may* occur in the edge interior.
		//
		// The maximum latitude is 90 degrees minus the latitude of N. We
		// compute this directly using atan2 in order to get maximum accuracy
		// near the poles.
		//
		// Our goal is compute a bound that contains the computed latitudes of
		// all S2Points P that pass the point-in-polygon containment test.
		// There are three sources of error we need to consider:
		// - the directional error in N (at most 3.84 * dblEpsilon)
		// - converting N to a maximum latitude
		// - computing the latitude of the test point P
		// The latter two sources of error are at most 0.955 * dblEpsilon
		// individually, but it is possible to show by a more complex analysis
		// that together they can add up to at most 1.16 * dblEpsilon, for a
		// total error of 5 * dblEpsilon.
		//
		// We add 3 * dblEpsilon to the bound here, and GetBound() will pad
		// the bound by another 2 * dblEpsilon.
		maxLat := math.Min(
			math.Atan2(math.Sqrt(n.X*n.X+n.Y*n.Y), math.Abs(n.Z))+3*dblEpsilon,
			math.Pi/2)

		// In order to get tight bounds when the two points are close together,
		// we also bound the min/max latitude relative to the latitudes of the
		// endpoints A and B. First we compute the distance between A and B,
		// and then we compute the maximum change in latitude between any two
		// points along the great circle that are separated by this distance.
		// This gives us a latitude change "budget". Some of this budget must
		// be spent getting from A to B; the remainder bounds the round-trip
		// distance (in latitude) from A or B to the min or max latitude
		// attained along the edge AB.
		latBudget := 2 * math.Asin(0.5*(r.a.Sub(b.Vector)).Norm()*math.Sin(maxLat))
		maxDelta := 0.5*(latBudget-latAB.Length()) + dblEpsilon

		// Test whether AB passes through the point of maximum latitude or
		// minimum latitude. If the dot product(s) are small enough then the
		// result may be ambiguous.
		if mA <= mError && mB >= -mError {
			latAB.Hi = math.Min(maxLat, latAB.Hi+maxDelta)
		}
		if mB <= mError && mA >= -mError {
			latAB.Lo = math.Max(-maxLat, latAB.Lo-maxDelta)
		}
	}
	r.a = b
	r.aLL = bLL
	r.bound = r.bound.Union(Rect{latAB, lngAB})
}
919260
func (p *SegmentedPath) LineTo(x, y float64) {
	p.Points = append(p.Points, x, y)
}
909092
func GeoPointFromLatLon(lat, lon float64) *GeoPoint {
	return &GeoPoint{Lat: lat, Lon: lon}
}
340424
func EuclideanDistance(dataPointX, dataPointY Float64Data) (distance float64, err error) {

	err = validateData(dataPointX, dataPointY)
	if err != nil {
		return math.NaN(), err
	}
	distance = 0
	for i := 0; i < len(dataPointX); i++ {
		distance = distance + ((dataPointX[i] - dataPointY[i]) * (dataPointX[i] - dataPointY[i]))
	}
	return math.Sqrt(distance), nil
}
10214
func (p Point) Distance(b Point) s1.Angle {
	return p.Vector.Angle(b.Vector)
}
1101075
func (p *Point) MidpointTo(p2 *Point) *Point {
	lat1 := p.lat * math.Pi / 180.0
	lat2 := p2.lat * math.Pi / 180.0

	lon1 := p.lng * math.Pi / 180.0
	dLon := (p2.lng - p.lng) * math.Pi / 180.0

	bx := math.Cos(lat2) * math.Cos(dLon)
	by := math.Cos(lat2) * math.Sin(dLon)

	lat3Rad := math.Atan2(
		math.Sin(lat1)+math.Sin(lat2),
		math.Sqrt(math.Pow(math.Cos(lat1)+bx, 2)+math.Pow(by, 2)),
	)
	lon3Rad := lon1 + math.Atan2(by, math.Cos(lat1)+bx)

	lat3 := lat3Rad * 180.0 / math.Pi
	lon3 := lon3Rad * 180.0 / math.Pi

	return NewPoint(lat3, lon3)
}
1357332
func Distance(p1, p2 orb.Point) float64 {
	d0 := (p1[0] - p2[0])
	d1 := (p1[1] - p2[1])
	return math.Sqrt(d0*d0 + d1*d1)
}
238422
func (s GeoDistanceSort) Point(lat, lon float64) GeoDistanceSort {
	s.points = append(s.points, GeoPointFromLatLon(lat, lon))
	return s
}
969429
func DistancePointToLine(point, lineStart, lineEnd geom.Coord) float64 {
	// if start = end, then just compute distance to one of the endpoints
	if Equals(lineStart, lineEnd) {
		return Distance(point, lineStart)
	}

	// otherwise use comp.graphics.algorithms Frequently Asked Questions method
	/*
	 * (1) r = AC dot AB
	 *         ---------
	 *         ||AB||^2
	 *
	 * r has the following meaning:
	 *   r=0 P = A
	 *   r=1 P = B
	 *   r<0 P is on the backward extension of AB
	 *   r>1 P is on the forward extension of AB
	 *   0<r<1 P is interior to AB
	 */

	len2 := (lineEnd[0]-lineStart[0])*(lineEnd[0]-lineStart[0]) + (lineEnd[1]-lineStart[1])*(lineEnd[1]-lineStart[1]) + (lineEnd[2]-lineStart[2])*(lineEnd[2]-lineStart[2])
	if math.IsNaN(len2) {
		panic("Ordinates must not be NaN")
	}
	r := ((point[0]-lineStart[0])*(lineEnd[0]-lineStart[0]) + (point[1]-lineStart[1])*(lineEnd[1]-lineStart[1]) + (point[2]-lineStart[2])*(lineEnd[2]-lineStart[2])) / len2

	if r <= 0.0 {
		return Distance(point, lineStart)
	}
	if r >= 1.0 {
		return Distance(point, lineEnd)
	}

	// compute closest point q on line segment
	qx := lineStart[0] + r*(lineEnd[0]-lineStart[0])
	qy := lineStart[1] + r*(lineEnd[1]-lineStart[1])
	qz := lineStart[2] + r*(lineEnd[2]-lineStart[2])
	// result is distance from p to q
	dx := point[0] - qx
	dy := point[1] - qy
	dz := point[2] - qz
	return math.Sqrt(dx*dx + dy*dy + dz*dz)
}
715303
func DistanceLineToLine(line1Start, line1End, line2Start, line2End geom.Coord) float64 {
	/**
	 * This calculation is susceptible to roundoff errors when
	 * passed large ordinate values.
	 * It may be possible to improve this by using {@link DD} arithmetic.
	 */
	if Equals(line1Start, line1End) {
		return DistancePointToLine(line1Start, line2Start, line2End)
	}
	if Equals(line2Start, line1End) {
		return DistancePointToLine(line2Start, line1Start, line1End)
	}

	/**
	 * Algorithm derived from http://softsurfer.com/Archive/algorithm_0106/algorithm_0106.htm
	 */
	a := VectorDot(line1Start, line1End, line1Start, line1End)
	b := VectorDot(line1Start, line1End, line2Start, line2End)
	c := VectorDot(line2Start, line2End, line2Start, line2End)
	d := VectorDot(line1Start, line1End, line2Start, line1Start)
	e := VectorDot(line2Start, line2End, line2Start, line1Start)

	denom := a*c - b*b
	if math.IsNaN(denom) {
		panic("Ordinates must not be NaN")
	}

	var s, t float64
	if denom <= 0.0 {
		/**
		 * The lines are parallel.
		 * In this case solve for the parameters s and t by assuming s is 0.
		 */
		s = 0
		// choose largest denominator for optimal numeric conditioning
		if b > c {
			t = d / b
		} else {
			t = e / c
		}
	} else {
		s = (b*e - c*d) / denom
		t = (a*e - b*d) / denom
	}
	switch {
	case s < 0:
		return DistancePointToLine(line1Start, line2Start, line2End)
	case s > 1:
		return DistancePointToLine(line1End, line2Start, line2End)
	case t < 0:
		return DistancePointToLine(line2Start, line1Start, line1End)
	case t > 1:
		return DistancePointToLine(line2End, line1Start, line1End)
	}
	/**
	 * The closest points are in interiors of segments,
	 * so compute them directly
	 */
	x1 := line1Start[0] + s*(line1End[0]-line1Start[0])
	y1 := line1Start[1] + s*(line1End[1]-line1Start[1])
	z1 := line1Start[2] + s*(line1End[2]-line1Start[2])

	x2 := line2Start[0] + t*(line2End[0]-line2Start[0])
	y2 := line2Start[1] + t*(line2End[1]-line2Start[1])
	z2 := line2Start[2] + t*(line2End[2]-line2Start[2])

	// length (p1-p2)
	return Distance(geom.Coord{x1, y1, z1}, geom.Coord{x2, y2, z2})
}
541260
func (c Cell) distanceInternal(targetXYZ Point, toInterior bool) s1.ChordAngle {
	// All calculations are done in the (u,v,w) coordinates of this cell's face.
	target := faceXYZtoUVW(int(c.face), targetXYZ)

	// Compute dot products with all four upward or rightward-facing edge
	// normals. dirIJ is the dot product for the edge corresponding to axis
	// I, endpoint J. For example, dir01 is the right edge of the Cell
	// (corresponding to the upper endpoint of the u-axis).
	dir00 := target.X - target.Z*c.uv.X.Lo
	dir01 := target.X - target.Z*c.uv.X.Hi
	dir10 := target.Y - target.Z*c.uv.Y.Lo
	dir11 := target.Y - target.Z*c.uv.Y.Hi
	inside := true
	if dir00 < 0 {
		inside = false // Target is to the left of the cell
		if c.vEdgeIsClosest(target, false) {
			return edgeDistance(-dir00, c.uv.X.Lo)
		}
	}
	if dir01 > 0 {
		inside = false // Target is to the right of the cell
		if c.vEdgeIsClosest(target, true) {
			return edgeDistance(dir01, c.uv.X.Hi)
		}
	}
	if dir10 < 0 {
		inside = false // Target is below the cell
		if c.uEdgeIsClosest(target, false) {
			return edgeDistance(-dir10, c.uv.Y.Lo)
		}
	}
	if dir11 > 0 {
		inside = false // Target is above the cell
		if c.uEdgeIsClosest(target, true) {
			return edgeDistance(dir11, c.uv.Y.Hi)
		}
	}
	if inside {
		if toInterior {
			return s1.ChordAngle(0)
		}
		// Although you might think of Cells as rectangles, they are actually
		// arbitrary quadrilaterals after they are projected onto the sphere.
		// Therefore the simplest approach is just to find the minimum distance to
		// any of the four edges.
		return minChordAngle(edgeDistance(-dir00, c.uv.X.Lo),
			edgeDistance(dir01, c.uv.X.Hi),
			edgeDistance(-dir10, c.uv.Y.Lo),
			edgeDistance(dir11, c.uv.Y.Hi))
	}

	// Otherwise, the closest point is one of the four cell vertices. Note that
	// it is *not* trivial to narrow down the candidates based on the edge sign
	// tests above, because (1) the edges don't meet at right angles and (2)
	// there are points on the far side of the sphere that are both above *and*
	// below the cell, etc.
	return minChordAngle(c.vertexChordDist2(target, false, false),
		c.vertexChordDist2(target, true, false),
		c.vertexChordDist2(target, false, true),
		c.vertexChordDist2(target, true, true))
}
1338489
func (calc *PointCentroidCalculator) AddPoint(point *geom.Point) {
	calc.AddCoord(geom.Coord(point.FlatCoords()))
}
443419
func distanceFromSegmentSquared(a, b, point []float64) float64 {
	x := a[0]
	y := a[1]
	dx := b[0] - x
	dy := b[1] - y

	if dx != 0 || dy != 0 {
		t := ((point[0]-x)*dx + (point[1]-y)*dy) / (dx*dx + dy*dy)

		if t > 1 {
			x = b[0]
			y = b[1]
		} else if t > 0 {
			x += dx * t
			y += dy * t
		}
	}

	dx = point[0] - x
	dy = point[1] - y

	return dx*dx + dy*dy
}