github.com/agnivade/pgm@v0.0.0-20210528073050-e2df0d9cb72d/hull.go (about) 1 package pgm 2 3 import ( 4 "math" 5 ) 6 7 type hull []point 8 9 type corner int 10 11 const ( 12 upperRight corner = iota + 1 13 upperLeft 14 lowerLeft 15 lowerRight 16 ) 17 18 func cross(a, b, o point) float64 { 19 return (a.x-o.x)*(b.y-o.y) - (a.y-o.y)*(b.x-o.x) 20 } 21 22 // Uses the monotone chain algorithm. 23 // https://en.wikibooks.org/wiki/Algorithm_Implementation/Geometry/Convex_hull/Monotone_chain 24 func buildHull(pts []point) hull { 25 // points will already be sorted. 26 var lower, upper []point 27 28 // Build lower hull 29 for i := 0; i < len(pts); i++ { 30 for len(lower) >= 2 && cross(lower[len(lower)-2], lower[len(lower)-1], pts[i]) <= 0 { 31 lower = lower[:len(lower)-1] 32 } 33 lower = append(lower, pts[i]) 34 } 35 36 // Build upper hull 37 for i := len(pts) - 1; i >= 0; i-- { 38 for len(upper) >= 2 && cross(upper[len(upper)-2], upper[len(upper)-1], pts[i]) <= 0 { 39 upper = upper[:len(upper)-1] 40 } 41 upper = append(upper, pts[i]) 42 } 43 44 lower = lower[:len(lower)-1] 45 upper = upper[:len(upper)-1] 46 return append(lower, upper...) 47 } 48 49 func (h hull) getSmallestRectangle() rectangle { 50 rectangles := []rectangle{} 51 i := newCaliper(h, h.getIndex(upperRight), 90) 52 j := newCaliper(h, h.getIndex(upperLeft), 180) 53 k := newCaliper(h, h.getIndex(lowerLeft), 270) 54 l := newCaliper(h, h.getIndex(lowerRight), 0) 55 56 for l.currentAngle < 90 { 57 rectangles = append(rectangles, newRectangle( 58 l.getIntersection(i), 59 i.getIntersection(j), 60 j.getIntersection(k), 61 k.getIntersection(l), 62 )) 63 64 smallestTheta := getSmallestTheta(i, j, k, l) 65 66 i.rotateBy(smallestTheta) 67 j.rotateBy(smallestTheta) 68 k.rotateBy(smallestTheta) 69 l.rotateBy(smallestTheta) 70 } 71 72 index := 0 73 area := math.MaxFloat64 74 75 for i, r := range rectangles { 76 tmp := r.getArea() 77 if tmp < area { 78 area = tmp 79 index = i 80 } 81 } 82 83 return rectangles[index] 84 } 85 86 func (h hull) getIndex(c corner) int { 87 index := 0 88 pt := h[index] 89 90 for i := 1; i < len(h); i++ { 91 tmp := h[i] 92 change := false 93 94 switch c { 95 case upperRight: 96 change = (tmp.x > pt.x || (tmp.x == pt.x && tmp.y > pt.y)) 97 case upperLeft: 98 change = (tmp.y > pt.y || (tmp.y == pt.y && tmp.x < pt.x)) 99 case lowerLeft: 100 change = (tmp.x < pt.x || (tmp.x == pt.x && tmp.y < pt.y)) 101 case lowerRight: 102 change = (tmp.y < pt.y || (tmp.y == pt.y && tmp.x > pt.x)) 103 } 104 105 if change { 106 index = i 107 pt = tmp 108 } 109 } 110 111 return index 112 }