1package authoring23import (4 "fmt"5 "math"67 "github.com/FelineStateMachine/atlas/format/vnext"8)910func transformGeometry(geometry vnext.Geometry, mapping GeometryMap, target CoordinateSpace) (vnext.Geometry, error) {11 if familyOf(geometry.Kind) != mapping.Family {12 return vnext.Geometry{}, fmt.Errorf("geometry is %s, feature set requires %s", familyOf(geometry.Kind), mapping.Family)13 }14 transformed := vnext.Geometry{Kind: geometry.Kind, Parts: make([]vnext.GeometryPart, len(geometry.Parts))}15 for partIndex, part := range geometry.Parts {16 transformed.Parts[partIndex].Rings = make([][]vnext.Position, len(part.Rings))17 for ringIndex, ring := range part.Rings {18 if geometry.Kind == vnext.GeometryPolygon && (len(ring) < 4 || ring[0] != ring[len(ring)-1]) {19 return vnext.Geometry{}, fmt.Errorf("polygon ring is not closed")20 }21 positions := make([]vnext.Position, len(ring))22 for positionIndex, position := range ring {23 positions[positionIndex] = transformPosition(position, mapping.Transform)24 }25 transformed.Parts[partIndex].Rings[ringIndex] = positions26 }27 }28 clipped := clipGeometry(transformed, target.Extent)29 if len(clipped.Parts) == 0 {30 return vnext.Geometry{}, fmt.Errorf("transformed geometry falls outside target extent")31 }32 return clampGeometry(clipped, target.Extent), nil33}3435// clampGeometry removes floating-point dust left by affine projection and36// segment intersection at an exact target edge. Clipping has already decided37// topology; this only makes its closed-bound promise exact on disk.38func clampGeometry(geometry vnext.Geometry, extent [4]float64) vnext.Geometry {39 for partIndex := range geometry.Parts {40 for ringIndex := range geometry.Parts[partIndex].Rings {41 for positionIndex, position := range geometry.Parts[partIndex].Rings[ringIndex] {42 geometry.Parts[partIndex].Rings[ringIndex][positionIndex] = vnext.Position{43 math.Max(extent[0], math.Min(extent[2], position[0])),44 math.Max(extent[1], math.Min(extent[3], position[1])),45 }46 }47 }48 }49 return geometry50}5152func clipGeometry(geometry vnext.Geometry, extent [4]float64) vnext.Geometry {53 out := vnext.Geometry{Kind: geometry.Kind}54 switch geometry.Kind {55 case vnext.GeometryPoint:56 for _, part := range geometry.Parts {57 if len(part.Rings) == 1 && len(part.Rings[0]) == 1 && inside(part.Rings[0][0], extent) {58 out.Parts = append(out.Parts, part)59 }60 }61 case vnext.GeometryLineString:62 for _, part := range geometry.Parts {63 if len(part.Rings) != 1 {64 continue65 }66 for _, line := range clipLine(part.Rings[0], extent) {67 out.Parts = append(out.Parts, vnext.GeometryPart{Rings: [][]vnext.Position{line}})68 }69 }70 case vnext.GeometryPolygon:71 for _, part := range geometry.Parts {72 if len(part.Rings) == 0 {73 continue74 }75 exterior := clipRing(part.Rings[0], extent)76 if len(exterior) == 0 {77 continue78 }79 rings := [][]vnext.Position{exterior}80 for _, hole := range part.Rings[1:] {81 if clipped := clipRing(hole, extent); len(clipped) > 0 {82 rings = append(rings, clipped)83 }84 }85 out.Parts = append(out.Parts, vnext.GeometryPart{Rings: rings})86 }87 }88 return out89}9091func clipLine(line []vnext.Position, extent [4]float64) [][]vnext.Position {92 var out [][]vnext.Position93 var current []vnext.Position94 flush := func() {95 if len(current) >= 2 {96 out = append(out, current)97 }98 current = nil99 }100 for index := 1; index < len(line); index++ {101 start, end, held := clipSegment(line[index-1], line[index], extent)102 if !held {103 flush()104 continue105 }106 if len(current) == 0 || current[len(current)-1] != start {107 flush()108 current = []vnext.Position{start}109 }110 if current[len(current)-1] != end {111 current = append(current, end)112 }113 }114 flush()115 return out116}117118func clipRing(ring []vnext.Position, extent [4]float64) []vnext.Position {119 if len(ring) < 4 || ring[0] != ring[len(ring)-1] {120 return nil121 }122 current := append([]vnext.Position(nil), ring[:len(ring)-1]...)123 insideEdge := [4]func(vnext.Position) bool{124 func(point vnext.Position) bool { return point[0] >= extent[0] },125 func(point vnext.Position) bool { return point[0] <= extent[2] },126 func(point vnext.Position) bool { return point[1] >= extent[1] },127 func(point vnext.Position) bool { return point[1] <= extent[3] },128 }129 crossing := [4]func(vnext.Position, vnext.Position) vnext.Position{130 func(from, to vnext.Position) vnext.Position { return crossX(from, to, extent[0]) },131 func(from, to vnext.Position) vnext.Position { return crossX(from, to, extent[2]) },132 func(from, to vnext.Position) vnext.Position { return crossY(from, to, extent[1]) },133 func(from, to vnext.Position) vnext.Position { return crossY(from, to, extent[3]) },134 }135 for edge := range insideEdge {136 if len(current) == 0 {137 return nil138 }139 next := make([]vnext.Position, 0, len(current)+4)140 for index, to := range current {141 from := current[(index+len(current)-1)%len(current)]142 fromInside, toInside := insideEdge[edge](from), insideEdge[edge](to)143 switch {144 case fromInside && toInside:145 next = append(next, to)146 case fromInside:147 next = append(next, crossing[edge](from, to))148 case toInside:149 next = append(next, crossing[edge](from, to), to)150 }151 }152 current = next153 }154 if len(current) < 3 {155 return nil156 }157 return append(current, current[0])158}159160func clipSegment(from, to vnext.Position, extent [4]float64) (vnext.Position, vnext.Position, bool) {161 enter, leave := 0.0, 1.0162 dx, dy := to[0]-from[0], to[1]-from[1]163 narrow := func(p, q float64) bool {164 if p == 0 {165 return q >= 0166 }167 at := q / p168 if p < 0 {169 if at > leave {170 return false171 }172 enter = math.Max(enter, at)173 return true174 }175 if at < enter {176 return false177 }178 leave = math.Min(leave, at)179 return true180 }181 if !narrow(-dx, from[0]-extent[0]) || !narrow(dx, extent[2]-from[0]) ||182 !narrow(-dy, from[1]-extent[1]) || !narrow(dy, extent[3]-from[1]) {183 return from, to, false184 }185 return vnext.Position{from[0] + enter*dx, from[1] + enter*dy},186 vnext.Position{from[0] + leave*dx, from[1] + leave*dy}, true187}188189func crossX(from, to vnext.Position, x float64) vnext.Position {190 at := (x - from[0]) / (to[0] - from[0])191 return vnext.Position{x, from[1] + at*(to[1]-from[1])}192}193194func crossY(from, to vnext.Position, y float64) vnext.Position {195 at := (y - from[1]) / (to[1] - from[1])196 return vnext.Position{from[0] + at*(to[0]-from[0]), y}197}198199func sourceExtent(mapping GeometryMap, target CoordinateSpace) ([4]float64, error) {200 if mapping.Transform.Kind == "identity" {201 return target.Extent, nil202 }203 matrix := mapping.Transform.Matrix204 determinant := matrix[0]*matrix[4] - matrix[1]*matrix[3]205 if determinant == 0 {206 return [4]float64{}, fmt.Errorf("affine transform is singular")207 }208 inverse := CoordinateTransform{Kind: "affine", Matrix: [6]float64{209 matrix[4] / determinant,210 -matrix[1] / determinant,211 (matrix[1]*matrix[5] - matrix[4]*matrix[2]) / determinant,212 -matrix[3] / determinant,213 matrix[0] / determinant,214 (matrix[3]*matrix[2] - matrix[0]*matrix[5]) / determinant,215 }}216 corners := []vnext.Position{217 {target.Extent[0], target.Extent[1]}, {target.Extent[0], target.Extent[3]},218 {target.Extent[2], target.Extent[1]}, {target.Extent[2], target.Extent[3]},219 }220 extent := [4]float64{math.Inf(1), math.Inf(1), math.Inf(-1), math.Inf(-1)}221 for _, corner := range corners {222 position := transformPosition(corner, inverse)223 extent[0], extent[1] = math.Min(extent[0], position[0]), math.Min(extent[1], position[1])224 extent[2], extent[3] = math.Max(extent[2], position[0]), math.Max(extent[3], position[1])225 }226 return extent, nil227}228229func transformPosition(position vnext.Position, transform CoordinateTransform) vnext.Position {230 if transform.Kind == "identity" {231 return position232 }233 matrix := transform.Matrix234 return vnext.Position{235 matrix[0]*position[0] + matrix[1]*position[1] + matrix[2],236 matrix[3]*position[0] + matrix[4]*position[1] + matrix[5],237 }238}239240func inside(position vnext.Position, extent [4]float64) bool {241 return finite(position[0]) && finite(position[1]) &&242 position[0] >= extent[0] && position[0] <= extent[2] &&243 position[1] >= extent[1] && position[1] <= extent[3]244}245246func familyOf(kind vnext.GeometryKind) string {247 switch kind {248 case vnext.GeometryPoint:249 return "point"250 case vnext.GeometryLineString:251 return "path"252 case vnext.GeometryPolygon:253 return "area"254 default:255 return "unknown"256 }257}