Calculate minimum bounding circle radius for geolocation coordinates using Welzl's algorithm

min_bounding_circle_radius.go

package helpers

import (
	"math"
	"math/rand"
)

type Point struct {
	Lat, Lng float64
}

type Circle struct {
	Center Point
	Radius float64
}

// Find the radius of the circle that contains all the given geo coordinates
// using the Min Bounding Circle algorithm
func MinBoundingCircleRadius(points *[]types.Point) float64 {
	return welzl(points).Radius
}

func welzl(points []Point) Circle {
	// Shuffle the points to ensure a random order
	rand.Shuffle(len(points), func(i, j int) { points[i], points[j] = points[j], points[i] })
	return welzlRecursive(points, []Point{})
}

func welzlRecursive(points, boundary []Point) Circle {
	if len(boundary) == 3 {
		return makeCircle(boundary)
	}
	if len(points) == 0 && len(boundary) == 2 {
		return makeCircle(boundary)
	}
	if len(points) == 0 || len(boundary) == 3 {
		if len(boundary) == 0 {
			return Circle{}
		} else if len(boundary) == 1 {
			return Circle{Center: boundary[0], Radius: 0}
		} else if len(boundary) == 2 {
			return makeCircle(boundary)
		}
	}

	// Try without the last point
	circle := welzlRecursive(points[:len(points)-1], boundary)

	if distance(points[len(points)-1], circle.Center) <= circle.Radius {
		return circle
	}

	// Try with the last point
	return welzlRecursive(points[:len(points)-1], append(boundary, points[len(points)-1]))
}

func makeCircle(boundary []Point) Circle {
	if len(boundary) == 0 {
		return Circle{}
	} else if len(boundary) == 1 {
		return Circle{Center: boundary[0], Radius: 0}
	} else if len(boundary) == 2 {
		center := Point{(boundary[0].Lat + boundary[1].Lat) / 2, (boundary[0].Lng + boundary[1].Lng) / 2}
		radius := distance(boundary[0], center)
		return Circle{Center: center, Radius: radius}
	}

	// Three points: Compute the circumcenter and circumradius
	a, b, c := boundary[0], boundary[1], boundary[2]
	d := 2 * (a.Lat*(b.Lng-c.Lng) + b.Lat*(c.Lng-a.Lng) + c.Lat*(a.Lng-b.Lng))
	ux := ((a.Lat*a.Lat+a.Lng*a.Lng)*(b.Lng-c.Lng) + (b.Lat*b.Lat+b.Lng*b.Lng)*(c.Lng-a.Lng) + (c.Lat*c.Lat+c.Lng*c.Lng)*(a.Lng-b.Lng)) / d
	uy := ((a.Lat*a.Lat+a.Lng*a.Lng)*(c.Lat-b.Lat) + (b.Lat*b.Lat+b.Lng*b.Lng)*(a.Lat-c.Lat) + (c.Lat*c.Lat+c.Lng*c.Lng)*(b.Lat-a.Lat)) / d
	center := Point{ux, uy}
	radius := distance(center, a)
	return Circle{Center: center, Radius: radius}
}

func distance(p1, p2 Point) float64 {
	lat1, lng1 := p1.Lat, p1.Lng
	lat2, lng2 := p2.Lat, p2.Lng
	radlat1 := math.Pi * lat1 / 180.0
	radlng1 := math.Pi * lng1 / 180.0
	radlat2 := math.Pi * lat2 / 180.0
	radlng2 := math.Pi * lng2 / 180.0
	dlat := radlat2 - radlat1
	dlng := radlng2 - radlng1
	a := math.Pow(math.Sin(dlat/2), 2) + math.Cos(radlat1)*math.Cos(radlat2)*math.Pow(math.Sin(dlng/2), 2)
	c := 2 * math.Atan2(math.Sqrt(a), math.Sqrt(1-a))
	return 6371000 * c // Earth radius in meters
}