func directedHausdorffDistance()

in s2/rect.go [519:601]
43 lines of code
12 McCabe index (conditional complexity)

func directedHausdorffDistance(lngDiff s1.Angle, a, b r1.Interval) s1.Angle {
	// By symmetry, we can assume a's longitude is 0 and b's longitude is
	// lngDiff. Call b's two endpoints bLo and bHi. Let H be the hemisphere
	// containing a and delimited by the longitude line of b. The Voronoi diagram
	// of b on H has three edges (portions of great circles) all orthogonal to b
	// and meeting at bLo cross bHi.
	// E1: (bLo, bLo cross bHi)
	// E2: (bHi, bLo cross bHi)
	// E3: (-bMid, bLo cross bHi), where bMid is the midpoint of b
	//
	// They subdivide H into three Voronoi regions. Depending on how longitude 0
	// (which contains edge a) intersects these regions, we distinguish two cases:
	// Case 1: it intersects three regions. This occurs when lngDiff <= π/2.
	// Case 2: it intersects only two regions. This occurs when lngDiff > π/2.
	//
	// In the first case, the directed Hausdorff distance to edge b can only be
	// realized by the following points on a:
	// A1: two endpoints of a.
	// A2: intersection of a with the equator, if b also intersects the equator.
	//
	// In the second case, the directed Hausdorff distance to edge b can only be
	// realized by the following points on a:
	// B1: two endpoints of a.
	// B2: intersection of a with E3
	// B3: farthest point from bLo to the interior of D, and farthest point from
	//     bHi to the interior of U, if any, where D (resp. U) is the portion
	//     of edge a below (resp. above) the intersection point from B2.

	if lngDiff < 0 {
		panic("impossible: negative lngDiff")
	}
	if lngDiff > math.Pi {
		panic("impossible: lngDiff > Pi")
	}

	if lngDiff == 0 {
		return s1.Angle(a.DirectedHausdorffDistance(b))
	}

	// Assumed longitude of b.
	bLng := lngDiff
	// Two endpoints of b.
	bLo := PointFromLatLng(LatLng{s1.Angle(b.Lo), bLng})
	bHi := PointFromLatLng(LatLng{s1.Angle(b.Hi), bLng})

	// Cases A1 and B1.
	aLo := PointFromLatLng(LatLng{s1.Angle(a.Lo), 0})
	aHi := PointFromLatLng(LatLng{s1.Angle(a.Hi), 0})
	maxDistance := maxAngle(
		DistanceFromSegment(aLo, bLo, bHi),
		DistanceFromSegment(aHi, bLo, bHi))

	if lngDiff <= math.Pi/2 {
		// Case A2.
		if a.Contains(0) && b.Contains(0) {
			maxDistance = maxAngle(maxDistance, lngDiff)
		}
		return maxDistance
	}

	// Case B2.
	p := bisectorIntersection(b, bLng)
	pLat := LatLngFromPoint(p).Lat
	if a.Contains(float64(pLat)) {
		maxDistance = maxAngle(maxDistance, p.Angle(bLo.Vector))
	}

	// Case B3.
	if pLat > s1.Angle(a.Lo) {
		intDist, ok := interiorMaxDistance(r1.Interval{a.Lo, math.Min(float64(pLat), a.Hi)}, bLo)
		if ok {
			maxDistance = maxAngle(maxDistance, intDist)
		}
	}
	if pLat < s1.Angle(a.Hi) {
		intDist, ok := interiorMaxDistance(r1.Interval{math.Max(float64(pLat), a.Lo), a.Hi}, bHi)
		if ok {
			maxDistance = maxAngle(maxDistance, intDist)
		}
	}

	return maxDistance
}