package border

import (
        "image"
)

const (
        Outer = 0
        Hole  = 1
)

// Contour represents a single contour/border extracted from an image.
// It also tracks its parents and children.
type Contour struct {

        // Points making up the contour
        Points []image.Point

        Id int

        // Outer or Hole.
        BorderType int

        // Id of parent
        ParentId int

        // ParentCollision indicates if colliding with parent. Just an optimisation for quick removal later on.
        ParentCollision bool

        // Parent links to contours parent
        Parent *Contour

        // Children links to contours children
        Children []*Contour

        // ConflictingContours is a map of contours that we KNOW we conflict with. This may be the parent or other
        // siblings
        ConflictingContours map[int]bool // hate to use maps here... but want uniqueness

        // usable or not. Not filtering out but marking that we may not use it. (say if we're conflicting with another contour)
        Usable bool
}

// NewContour create new contour
func NewContour(id int) *Contour {
        c := Contour{}
        c.Id = id
        c.BorderType = Hole
        c.ConflictingContours = make(map[int]bool)
        c.Usable = true
        return &c
}

// GetAllPoints returns all points in the contour and all children.
func (c *Contour) GetAllPoints() []image.Point {

        var allPoints []image.Point

        for _, p := range c.Points {
                allPoints = append(allPoints, p)
        }

        for _, ch := range c.Children {
                points := ch.GetAllPoints()
                allPoints = append(allPoints, points...)
        }

        return allPoints
}

package border

import (
        "errors"
        "image"

        "github.com/kpfaulkner/borders/common"
)

// dirDelta determines which direction we move based on the direction (0-7) index.
// Fixed-size array so the type itself documents the 8-direction invariant and the
// compiler can prove dirDelta[dir] is in bounds for any `dir` it can range-check.
var dirDelta = [8]image.Point{{0, -1}, {1, -1}, {1, 0}, {1, 1}, {0, 1}, {-1, 1}, {-1, 0}, {-1, -1}}

// FindContours takes a SuzukiImage and determines the Contours that are present.
// It returns the single parent contour which in turn has all other contours as children or further
// generations.
func FindContours(img *common.SuzukiImage) (*Contour, error) {
        nbd := 1
        lnbd := 1

        // contours[id] holds the contour with that id. Contour ids are assigned
        // sequentially starting at 1, so direct slice indexing replaces a map.
        // Index 0 is a nil sentinel — no contour ever has id 0.
        contours := []*Contour{nil, NewContour(1)}

        height := img.Height
        width := img.Width

        for i := 0; i < height; i++ {
                lnbd = 1
                for j := 0; j < width; j++ {
                        fji := img.GetXY(j, i)
                        isOuter := fji == 1 && (j == 0 || img.GetXY(j-1, i) == 0)
                        isHole := fji >= 1 && (j == width-1 || img.GetXY(j+1, i) == 0)
                        if isOuter || isHole {

                                var contourPrime *Contour
                                contour := NewContour(1)
                                from := image.Point{j, i}
                                parentId := 0
                                if isOuter {
                                        nbd += 1
                                        from = from.Sub(image.Point{1, 0})
                                        contour.BorderType = Outer
                                        contourPrime = contours[lnbd]
                                        if contourPrime.BorderType == Outer {
                                                parentId = contourPrime.ParentId
                                        } else {
                                                parentId = contourPrime.Id
                                        }
                                } else {
                                        nbd += 1
                                        if fji > 1 {
                                                lnbd = fji
                                        }
                                        contourPrime = contours[lnbd]
                                        from = from.Add(image.Point{1, 0})
                                        contour.BorderType = Hole
                                        if contourPrime.BorderType == Outer {
                                                parentId = contourPrime.Id
                                        } else {
                                                parentId = contourPrime.ParentId
                                        }
                                }

                                p0 := image.Point{j, i}
                                border, collectionIndices, err := createBorder(img, p0, from, nbd)
                                if err != nil {
                                        return nil, err
                                }

                                if len(border) == 0 {
                                        border = append(border, p0)
                                        img.Set(p0, -1*nbd)
                                }

                                if parentId != 0 {
                                        parent := contours[parentId]
                                        parent.Children = append(parent.Children, contour)
                                        contour.Parent = contours[parentId]
                                }
                                contour.ParentId = parentId
                                contour.Points = border
                                contour.Id = nbd
                                // nbd was incremented to len(contours), so append lands the new contour at index nbd.
                                contours = append(contours, contour)
                                addCollisionFlag(contour, parentId, contours, collectionIndices)
                        }
                        if fji != 0 && fji != 1 {
                                lnbd = fji
                                if lnbd < 0 {
                                        lnbd *= -1
                                }
                        }
                }
        }

        finalContour := contours[1]

        // image was padded... so now shift every co-ord by -1,-1
        if img.HasPadding() {
                shiftContour(finalContour)
        }
        return finalContour, nil
}

// clockwise determines direction if we have 'dir' and turn clockwise
func clockwise(dir int) int {
        return (dir + 1) % 8
}

// counterClockwise determines direction if we have 'dir' and turn counterclockwise
func counterClockwise(dir int) int {
        return (dir + 7) % 8
}

// move moves the current point (pixel) in the direction 'dir'.
// Returns ok=false when the neighbour is out of bounds or background.
func move(pixel image.Point, img *common.SuzukiImage, dir int) (image.Point, bool) {
        newP := pixel.Add(dirDelta[dir])
        if newP.X < 0 || newP.X >= img.Width || newP.Y < 0 || newP.Y >= img.Height {
                return image.Point{}, false
        }
        if img.Get(newP) == 0 {
                return image.Point{}, false
        }
        return newP, true
}

// calcDir returns index of dirDelta that matches direction taken.
func calcDir(from image.Point, to image.Point) (int, error) {
        delta := to.Sub(from)
        for i, d := range dirDelta {
                if d.X == delta.X && d.Y == delta.Y {
                        return i, nil
                }
        }

        return 0, errors.New("unable to determine direction")
}

// createBorder returns the slice of Points making up the border/contour
// Also returns list of nbd's that are colliding with this. Can use to help create
// tree with collision info later.
func createBorder(img *common.SuzukiImage, p0 image.Point, p2 image.Point, nbd int) ([]image.Point, map[int]bool, error) {

        // track which borders have conflicts
        collisionIndices := make(map[int]bool)

        border := []image.Point{}
        dir, err := calcDir(p0, p2)
        if err != nil {
                return nil, nil, err
        }

        moved := clockwise(dir)
        var p1 image.Point
        foundP1 := false
        for moved != dir {
                if newP, ok := move(p0, img, moved); ok {
                        p1 = newP
                        foundP1 = true
                        break
                }
                moved = clockwise(moved)
        }

        if !foundP1 {
                return []image.Point{}, collisionIndices, nil
        }
        p2 = p1
        p3 := p0

        for {
                dir, err = calcDir(p3, p2)
                if err != nil {
                        return nil, nil, err
                }
                moved = counterClockwise(dir)
                var p4 image.Point
                var done [8]bool
                for {
                        var ok bool
                        p4, ok = move(p3, img, moved)
                        if ok {
                                break
                        }
                        done[moved] = true
                        moved = counterClockwise(moved)
                }

                // detect if colliding with something else (ie not 0 nor 1)
                curP3 := img.Get(p3)
                if curP3 != 1 {
                        if curP3 < 0 {
                                curP3 *= -1
                        }

                        absNbd := nbd
                        if absNbd < 0 {
                                absNbd *= -1
                        }
                        collisionIndices[curP3] = true
                        collisionIndices[absNbd] = true
                }

                border = append(border, p3)
                if p3.Y == img.Height-1 || done[2] {
                        img.Set(p3, -1*nbd)
                } else if img.Get(p3) == 1 {
                        img.Set(p3, nbd)
                }

                if p4.X == p0.X && p4.Y == p0.Y && p3.X == p1.X && p3.Y == p1.Y {
                        break
                }

                p2 = p3
                p3 = p4
        }

        return border, collisionIndices, nil
}

// addCollisionFlag mark contours with collisions with other contours.
func addCollisionFlag(contour *Contour, parentId int, contours []*Contour, collisionIndices map[int]bool) {
        for contour1 := range collisionIndices {

                // quick indicator to say colliding with parent.
                if contour1 == parentId {
                        contour.ParentCollision = true
                }

                for contour2 := range collisionIndices {
                        if contour1 != contour2 {
                                c1 := contours[contour1]
                                c1.ConflictingContours[contour2] = true
                        }
                }
        }
}

func shiftContour(contour *Contour) {
        for i, _ := range contour.Points {
                contour.Points[i].X = contour.Points[i].X - 1
                contour.Points[i].Y = contour.Points[i].Y - 1
        }

        for _, child := range contour.Children {
                shiftContour(child)
        }
}

package border

import (
        "fmt"
        "image"
        "image/color"
        "image/png"
        _ "image/png"
        "os"

        "github.com/kpfaulkner/borders/common"
        image2 "github.com/kpfaulkner/borders/image"
)

// LoadImage loads a PNG and returns a SuzukiImage. Currently restricted to PNG but will eventually expand
// to include other formats.
//
// Erode parameter forces the eroding of the image before converting to a SuzukiImage.
// See https://en.wikipedia.org/wiki/Erosion_(morphology) for explanation
//
// Dilate parameter forces the dilating of the image before converting to a SuzukiImage. Likewise, see
// https://en.wikipedia.org/wiki/Dilation_(morphology) for explanation.
//
// The combination of Erode and Dilate helps remove any "spikes" that may appear in the generated boundary.
// erode and dilate will usually be 0 (none) or 1 (single pixel spikes)
func LoadImage(filename string, erode int, dilate int) (*common.SuzukiImage, error) {

        f, err := os.Open(filename)
        if err != nil {
                return nil, err
        }
        defer f.Close()

        img, _, err := image.Decode(f)
        if err != nil {
                return nil, err
        }

        // If any pixels on the edges are populated, then we need to pad this out by 1 pixel on each side.
        // This will be reversed later.
        requirePadding := doesImageRequirePadding(img)

        // need border to be black. Pad edges with 1 black pixel
        si := common.NewSuzukiImage(img.Bounds().Dx(), img.Bounds().Dy(), requirePadding)

        paddingOffset := 0
        if requirePadding {
                paddingOffset = 1
        }

        // dumb... but convert to own image format for now.
        for y := 0; y < img.Bounds().Dy(); y++ {
                for x := 0; x < img.Bounds().Dx(); x++ {
                        cc := 0
                        c := img.At(x, y)
                        r, g, b, _ := c.RGBA()
                        if !(r == 0 && g == 0 && b == 0) {
                                cc = 1
                        }
                        si.SetXY(x+paddingOffset, y+paddingOffset, cc)
                }

        }

        if erode != 0 {
                si, err = image2.Erode(si, erode)
                if err != nil {
                        return nil, err
                }
        }

        if dilate != 0 {
                si, err = image2.Dilate(si, dilate)
                if err != nil {
                        return nil, err
                }
        }

        return si, nil
}

// check down each edge to see if populated, if so, it will require padding
func doesImageRequirePadding(img image.Image) bool {

        // down left/right edge
        for y := 0; y < img.Bounds().Dy(); y++ {
                leftEdgeX := 0
                c := img.At(leftEdgeX, y)
                r, g, b, _ := c.RGBA()
                if !(r == 0 && g == 0 && b == 0) {
                        return true
                }

                rightEdgeX := img.Bounds().Dx() - 1
                c = img.At(rightEdgeX, y)
                r, g, b, _ = c.RGBA()
                if !(r == 0 && g == 0 && b == 0) {
                        return true
                }
        }

        // across top and bottom
        for x := 0; x < img.Bounds().Dx(); x++ {
                topEdgeY := 0
                c := img.At(x, topEdgeY)
                r, g, b, _ := c.RGBA()
                if !(r == 0 && g == 0 && b == 0) {
                        return true
                }

                bottomEdgeY := img.Bounds().Dy() - 1
                c = img.At(x, bottomEdgeY)
                r, g, b, _ = c.RGBA()
                if !(r == 0 && g == 0 && b == 0) {
                        return true
                }
        }

        return false
}

// SaveImage saves a SuzukiImage as a PNG to given filename.
// Although currently only PNG, will make this more generic in the future.
func SaveImage(filename string, si *common.SuzukiImage) error {

        upLeft := image.Point{0, 0}
        lowRight := image.Point{si.Width, si.Height}

        img := image.NewRGBA(image.Rectangle{upLeft, lowRight})

        for x := 0; x < si.Width; x++ {
                for y := 0; y < si.Height; y++ {
                        p := si.GetXY(x, y)
                        if p == 1 {
                                img.Set(x, y, color.White)
                        } else {
                                img.Set(x, y, color.Black)
                        }
                }
        }

        f, err := os.Create(filename)
        if err != nil {
                return err
        }
        defer f.Close()
        return png.Encode(f, img)
}

// SaveContourSliceImage saves a contour (and all child contours) as a PNG.
// Width and height are the dimensions of the image to save.
//
// flipBook is a bool to indicate that when each child contour is added to the image, the image should be
// saved to a new file with the filename suffixed with the count of contours added so far. This is useful
// for debugging and visualising the contours as they are added.
//
// minContourSize indicates if minimum number of points that make up a contour. If contour contains fewer, then
// do NOT save.
func SaveContourSliceImage(filename string, c *Contour, width int, height int, flipBook bool, minContourSize int) (err error) {
        upLeft := image.Point{0, 0}
        lowRight := image.Point{width, height}

        img := image.NewRGBA(image.Rectangle{upLeft, lowRight})

        // naive fill
        for x := 0; x < width; x++ {
                for y := 0; y < height; y++ {
                        img.Set(x, y, color.Black)
                }
        }

        colour := 0
        count := 0

        if err := drawContour(img, c, flipBook, minContourSize, colour, &count, filename); err != nil {
                return err
        }
        f, err := os.Create(filename)
        if err != nil {
                return err
        }
        defer func() {
                err = f.Close()
        }()

        return png.Encode(f, img)
}

// drawContour saves a contour to the provided image and then recursively calls to save children to same image.
func drawContour(img *image.RGBA, c *Contour, flipBook bool, minContourSize int, colour int, count *int, filename string) error {

        colours := []color.RGBA{
                {255, 0, 0, 255},
                {255, 106, 0, 255},
                {255, 216, 0, 255},
                {0, 255, 0, 255},
                {127, 255, 197, 255},
                {72, 0, 255, 255},
                {255, 127, 182, 255},
        }

        max := len(colours)
        if c.BorderType == Outer {
                colour = 0
        }

        // draw contour itself.
        if len(c.Points) > 0 && len(c.Points) > minContourSize {

                colourToUse := colours[colour]
                for _, p := range c.Points {
                        img.Set(p.X, p.Y, colourToUse)
                }
                colour++
                if colour >= max {
                        colour = 0
                }

                // save new image per contour added...  crazy
                if flipBook {
                        fn := fmt.Sprintf("%s-%d.png", filename, *count)
                        f, err := os.Create(fn)
                        if err != nil {
                                return err
                        }
                        encErr := png.Encode(f, img)
                        if cErr := f.Close(); encErr == nil {
                                encErr = cErr
                        }
                        if encErr != nil {
                                return encErr
                        }
                }
                *count = *count + 1
        }

        for _, child := range c.Children {
                colour++
                if colour >= max {
                        colour = 0
                }
                *count = *count + 1
                if err := drawContour(img, child, flipBook, minContourSize, colour, count, filename); err != nil {
                        return err
                }
        }

        return nil
}

package common

import (
        "fmt"
        "image"
        "strings"
)

// SuzukiImage is the basic structure we use to define an image when trying to find contours.
//
// Storage is int32 per cell: each cell holds either 0/1 (background/foreground)
// or a signed contour id assigned during border following. Contour ids are
// bounded by the pixel count, so int32 is sufficient for any image up to ~2^31
// pixels and halves memory vs. int on 64-bit platforms.
type SuzukiImage struct {
        Width   int
        Height  int
        data    []int32
        dataLen int

        // Indicates if a 1 pixel padding has been applied to around the image.
        // This helps with imagery where it goes RIGHT up to the edge.
        hasPadding bool
}

// NewSuzukiImage creates a new SuzukiImage of specific dimensions.
func NewSuzukiImage(width int, height int, hasPadding bool) *SuzukiImage {
        si := SuzukiImage{}
        padding := 0
        if hasPadding {
                padding = 2
        }
        si.Width = width + padding
        si.Height = height + padding
        si.data = make([]int32, si.Width*si.Height)
        si.dataLen = si.Width * si.Height // just saves us calculating a lot
        si.hasPadding = hasPadding
        return &si
}

func NewSuzukiImageFromData(width int, height int, hasPadding bool, data []int) *SuzukiImage {
        si := NewSuzukiImage(width, height, hasPadding)
        if hasPadding {
                // Copy each row of the unpadded source into the interior of the padded
                // buffer, leaving the 1-pixel zero border intact.
                for y := 0; y < height; y++ {
                        dst := (y+1)*si.Width + 1
                        src := y * width
                        for i := 0; i < width; i++ {
                                si.data[dst+i] = int32(data[src+i])
                        }
                }
        } else {
                for i, v := range data {
                        si.data[i] = int32(v)
                }
        }
        return si
}

// GetAllData returns a copy of the underlying data as []int.
func (si *SuzukiImage) GetAllData() []int {
        out := make([]int, len(si.data))
        for i, v := range si.data {
                out[i] = int(v)
        }
        return out
}

// Get returns the value of a given point
func (si *SuzukiImage) Get(p image.Point) int {
        return int(si.data[p.Y*si.Width+p.X])
}

// GetXY returns the value of a given x/y
func (si *SuzukiImage) GetXY(x int, y int) int {
        return int(si.data[y*si.Width+x])
}

// Set sets the value at a given point
func (si *SuzukiImage) Set(p image.Point, val int) {
        si.data[p.Y*si.Width+p.X] = int32(val)
}

// SetXY sets the value at a given x/y
func (si *SuzukiImage) SetXY(x int, y int, val int) {
        si.data[y*si.Width+x] = int32(val)
}

func (si *SuzukiImage) HasPadding() bool {
        return si.hasPadding
}

// DisplayAsText generates a string of a given image. This is purely used for debugging SMALL images
func (si *SuzukiImage) DisplayAsText() []string {
        s := []string{}
        for y := 0; y < si.Height; y++ {
                ss := si.data[y*si.Width : (y*si.Width + si.Width)]
                t := []string{}
                for _, i := range ss {
                        t = append(t, fmt.Sprintf("%d", i))
                }
                s = append(s, strings.Join(t, " ")+"\n")
        }

        return s
}

// Equals checks if two SuzukiImages are equal.
func (si *SuzukiImage) Equals(other *SuzukiImage) bool {
        if si.Width != other.Width || si.Height != other.Height {
                return false
        }

        for i := 0; i < si.dataLen; i++ {
                if si.data[i] != other.data[i] {
                        return false
                }
        }
        return true
}

package converters

import (
        "errors"
        "image"
        "math"

        "github.com/kpfaulkner/borders/border"
        "github.com/peterstace/simplefeatures/geom"
)

const (
        EarthRadius       = 6378137.0
        toleranceInMetres = 2

        degreesToRadiansRatio = math.Pi / 180.0
        radiansToDegreesRatio = 180.0 / math.Pi
)

type PointConverter func(x float64, y float64) (float64, float64)

// NewSlippyToLatLongConverter returns a function that converts slippy tile coordinates to lat/long.
func NewSlippyToLatLongConverter(slippyXOffset float64, slippyYOffset float64, scale int) func(X float64, Y float64) (float64, float64) {
        latLongN := math.Pow(2, float64(scale))
        f := func(x float64, y float64) (float64, float64) {
                long, lat := slippyCoordsToLongLat(slippyXOffset, slippyYOffset, x, y, latLongN)
                return long, lat
        }
        return f
}

// LatLongToSlippy converts lat/long to slippy tile coordinates.
func LatLongToSlippy(latDegrees float64, longDegrees float64, scale int) (float64, float64) {
        n := math.Exp2(float64(scale))
        x := int(math.Floor((longDegrees + 180.0) / 360.0 * n))
        if float64(x) >= n {
                x = int(n - 1)
        }
        y := int(math.Floor((1.0 - math.Log(math.Tan(latDegrees*math.Pi/180.0)+1.0/math.Cos(latDegrees*math.Pi/180.0))/math.Pi) / 2.0 * n))
        return float64(x), float64(y)
}

// ConvertContourToPolygon converts the contours (set of x/y coords) to geometries commonly used in the GIS space
// Convert to polygons, then simplify (if required) while still in "pixel space"
// Only then apply conversions which may be to lat/long (or any other conversions).
// Simplifying while in "pixel space" simplifies the simplification degTolerance calculation.
// params:
//
//                simplify: Simplify the resulting polygons
//            minPoints: Minimum number of vertices for a polygon to be considered valid. If less than this, then will be discarded. 0 means no minimum
//                degTolerance: Tolerance in pixels when simplifying. If set to 0, then will use defaults.
//                multiPolygonOnly: If the geometry is results in a GeometryCollection, then extract out the multipolygon part and return that.
//                pointConverters: Used to convert point co-ord systems. eg. slippy to lat/long.
func ConvertContourToPolygon(c *border.Contour, scale int, simplify bool, minPoints int, tolerance float64, multiPolygonOnly bool, pointConverters ...PointConverter) (*geom.Geometry, error) {
        polygons := []geom.Polygon{}

        err := convertContourToPolygons(c, minPoints, &polygons)
        if err != nil {
                return nil, err
        }

        mp := geom.NewMultiPolygon(polygons)

        if simplify {
                if tolerance == 0 {
                        tolerance = generateSimplifyTolerance(scale)
                }
                gg := mp.AsGeometry()
                simplifiedGeom, err := gg.Simplify(tolerance, geom.NoValidate{})
                if err != nil {
                        return nil, err
                }

                if multiPolygonOnly {
                        if simplifiedGeom.Type() == geom.TypeMultiPolygon {
                                mp, _ = simplifiedGeom.AsMultiPolygon()
                                return returnConvertedGeometry(&mp, pointConverters...), nil
                        }
                        return nil, errors.New("unable to filter multipolygon from geometry collection")
                }

                mp, ok := simplifiedGeom.AsMultiPolygon()
                if ok {
                        return returnConvertedGeometry(&mp, pointConverters...), nil
                } else {
                        return nil, errors.New("unable to convert simplified geom to multipolygon")
                }
        }
        return returnConvertedGeometry(&mp, pointConverters...), nil
}

// returnConvertedGeometry converts the multipolygon with PointConverters (if supplied)
// Can be used to help convert to lat/long or any other co-ordinate system.
func returnConvertedGeometry(mp *geom.MultiPolygon, pointConverters ...PointConverter) *geom.Geometry {
        finalMultiPoly := convertCoords(mp, pointConverters...)
        g := finalMultiPoly.AsGeometry()
        return &g
}

// convertCoords converts the coordinates of a multipolygon using the supplied PointConverters.
func convertCoords(mp *geom.MultiPolygon, converters ...PointConverter) *geom.MultiPolygon {
        mp2 := mp.TransformXY(func(xy geom.XY) geom.XY {
                x := xy.X
                y := xy.Y
                // run through converters.
                for _, converter := range converters {
                        newX, newY := converter(x, y)
                        x = newX
                        y = newY
                }
                return geom.XY{X: x, Y: y}
        })
        return &mp2
}

// generateLineString generates a LineString from a slice of image.Points.
func generateLineString(points []image.Point) *geom.LineString {
        seq := pointsToSequence(points)

        if seq.Length() > 2 {
                ls := geom.NewLineString(seq)
                return &ls
        }

        return &geom.LineString{}
}

// convertContourToPolygons converts the contour to a set of polygons but does NOT convert to different co-ord systems.
// A polygon is discarded when its outer ring has fewer than minPoints vertices. minPoints == 0 disables the filter.
// Contours with an empty outer ring are always skipped, and empty hole rings are dropped so they don't poison the polygon.
func convertContourToPolygons(c *border.Contour, minPoints int, polygons *[]geom.Polygon) error {
        if c.BorderType == border.Outer && len(c.Points) > 0 && (minPoints == 0 || len(c.Points) >= minPoints) {

                lineStrings := []geom.LineString{}
                outerLS := generateLineString(c.Points)
                lineStrings = append(lineStrings, *outerLS)

                // holes — skip any with no points, which would otherwise be appended as an empty LineString.
                for _, child := range c.Children {
                        if !child.ParentCollision && child.Usable && len(child.Points) > 0 {
                                ls := generateLineString(child.Points)
                                lineStrings = append(lineStrings, *ls)
                        }
                }

                poly := geom.NewPolygon(lineStrings)
                *polygons = append(*polygons, poly)
        }

        for _, child := range c.Children {
                // only process child if no conflict with parent.
                if !child.ParentCollision && child.Usable {
                        err := convertContourToPolygons(child, minPoints, polygons)
                        if err != nil {
                                return err
                        }
                }
        }

        return nil
}

// pointsToSequence converts a slice of image.Points to a geom.Sequence.
func pointsToSequence(points []image.Point) geom.Sequence {
        s := len(points)*2 + 2
        seq := make([]float64, s, s)
        index := 0
        for _, origP := range points {
                x, y := float64(origP.X), float64(origP.Y)
                seq[index] = x
                seq[index+1] = y
                index += 2
        }

        seq[index] = seq[0]
        seq[index+1] = seq[1]
        return geom.NewSequence(seq, geom.DimXY)
}

// slippyCoordsToLongLat converts to lat/long... and requires the slippy offset of top left corner of area.
func slippyCoordsToLongLat(slippyXOffset float64, slippyYOffset float64, xTile float64, yTile float64, latLongN float64) (float64, float64) {
        x := xTile + slippyXOffset
        y := yTile + slippyYOffset

        longDeg := (x/latLongN)*360.0 - 180.0
        latRad := math.Atan(math.Sinh(math.Pi - (y/latLongN)*2*math.Pi))
        latDeg := latRad * (180.0 / math.Pi)

        return longDeg, latDeg
}

// generateSimplifyTolerance will mainly be used when we want to convert to geographical co-ordinates
// By default we will determine how many metres per pixel (for input scale/zoom) and double it.
func generateSimplifyTolerance(scale int) float64 {
        mtrPerPixel := metresPerPixel(scale)
        tolerance := mtrPerPixel * toleranceInMetres
        return tolerance
}

// tileSizeInMetres is the size of a tile in metres.
func tileSizeInMetres(scale int) float64 {
        return 2 * math.Pi * EarthRadius / float64(uint64(1)<<uint64(scale))
}

// metresPerPixel is number of metres for a given input pixel. This is based on the scale/zoom.
func metresPerPixel(scale int) float64 {
        return tileSizeInMetres(scale) / 256.0
}

// filterMultiPolygonFromGeometryCollection currently unused. Will be used in upcoming version.
func filterMultiPolygonFromGeometryCollection(col *geom.GeometryCollection) (*geom.MultiPolygon, error) {
        var mp geom.MultiPolygon
        var ok bool
        for i := 0; i < col.NumGeometries(); i++ {
                g := col.GeometryN(i)
                mp, ok = g.AsMultiPolygon()
                if ok {
                        return &mp, nil
                }
        }

        return nil, errors.New("no multipolygon found in geometry collection")
}

// NewPixelToLatLongConverter returns a function that converts pixel coordinates to lat/long.
// topLeftPixelLat/topLeftPixelLong are the geographic coordinates of the image's top-left pixel.
// Important to note that the converter function returned when executed will return
// (longitude,latitude) in that order.
// Process is:
//
// 1) get X,Y coordinates for the topleft pixel
// 2) For each x,y coords passed (which will be position within image), convert to global space (add globalX/globalY)
// 3) Convert each globally positioned pixel to lat/long using the precomputed scale-dependent constants.
func NewPixelToLatLongConverter(topLeftPixelLat float64, topLeftPixelLong float64, scale int) func(X float64, Y float64) (float64, float64) {

        // Scale-dependent constants — computed once here rather than per call.
        // PixelXYToLatLong recomputes these every invocation; the closure below is
        // called once per polygon vertex (often tens of thousands of times per
        // conversion), so hoisting the math.Exp2 and three divisions matters.
        pixelGlobeSize := 256.0 * math.Exp2(float64(scale))
        halfPixelGlobeSize := pixelGlobeSize / 2.0
        xPixelsToDegreesRatio := pixelGlobeSize / 360.0
        yPixelsToRadiansRatio := pixelGlobeSize / (2.0 * math.Pi)

        // global pixel position of top left corner.
        gX, gY := LatLongToPixelXY(topLeftPixelLat, topLeftPixelLong, scale)
        globalX := float64(gX)
        globalY := float64(gY)

        return func(x float64, y float64) (float64, float64) {
                pixelX := x + globalX
                pixelY := y + globalY
                longitude := (pixelX - halfPixelGlobeSize) / xPixelsToDegreesRatio
                latitude := (2*math.Atan(math.Exp((pixelY-halfPixelGlobeSize)/(-yPixelsToRadiansRatio))) - math.Pi/2.0) * radiansToDegreesRatio
                return longitude, latitude
        }
}

func PixelXYToLatLong(pixelX int64, pixelY int64, scale int) (float64, float64) {

        pixelTileSize := 256.0
        pixelGlobeSize := pixelTileSize * math.Pow(2, float64(scale))
        xPixelsToDegreesRatio := pixelGlobeSize / 360.0
        yPixelsToRadiansRatio := pixelGlobeSize / (2.0 * math.Pi)
        halfPixelGlobeSize := pixelGlobeSize / 2.0

        longitude := (float64(pixelX) - halfPixelGlobeSize) / xPixelsToDegreesRatio
        latitude := (2*math.Atan(math.Exp((float64(pixelY)-halfPixelGlobeSize)/(-yPixelsToRadiansRatio))) -
                math.Pi/2.0) * radiansToDegreesRatio

        return latitude, longitude
}

func LatLongToPixelXY(latitude float64, longitude float64, scale int) (int64, int64) {

        pixelTileSize := 256.0
        pixelGlobeSize := pixelTileSize * math.Pow(2, float64(scale))
        xPixelsToDegreesRatio := pixelGlobeSize / 360.0
        yPixelsToRadiansRatio := pixelGlobeSize / (2.0 * math.Pi)
        halfPixelGlobeSize := pixelGlobeSize / 2.0

        x := math.Round(halfPixelGlobeSize + (longitude * xPixelsToDegreesRatio))
        f := math.Min(math.Max(math.Sin(latitude*degreesToRadiansRatio), -0.9999), 0.9999)
        y := math.Round(halfPixelGlobeSize + 0.5*math.Log((1+f)/(1-f))*(-yPixelsToRadiansRatio))
        return int64(x), int64(y)

}

package image

import (
        "github.com/kpfaulkner/borders/common"
)

// Erode applies morphological erosion with a square (2*radius+1)×(2*radius+1)
// structuring element. The outermost 1-pixel border is treated as zero, so the
// result's outermost border is always zero. Input is not modified.
//
// Two separable 1-D passes — O(W·H·r) work instead of O(W·H·r²).
// https://en.wikipedia.org/wiki/Erosion_(morphology)
func Erode(img *common.SuzukiImage, radius int) (*common.SuzukiImage, error) {
        width, height := img.Width, img.Height
        tmp := common.NewSuzukiImage(width, height, img.HasPadding())
        out := common.NewSuzukiImage(width, height, img.HasPadding())

        // x is bounded so the window never touches column 0 or width-1, which are
        // treated as permanently zero — that leaves the outer band of the output
        // eroded against the zero edge.
        for y := 0; y < height; y++ {
                for x := radius + 1; x < width-1-radius; x++ {
                        v := 1
                        for k := x - radius; k <= x+radius; k++ {
                                if img.GetXY(k, y) != 1 {
                                        v = 0
                                        break
                                }
                        }
                        if v == 1 {
                                tmp.SetXY(x, y, 1)
                        }
                }
        }

        for y := radius + 1; y < height-1-radius; y++ {
                for x := 0; x < width; x++ {
                        v := 1
                        for k := y - radius; k <= y+radius; k++ {
                                if tmp.GetXY(x, k) != 1 {
                                        v = 0
                                        break
                                }
                        }
                        if v == 1 {
                                out.SetXY(x, y, 1)
                        }
                }
        }

        return out, nil
}

// Dilate applies morphological dilation with a square (2*radius+1)×(2*radius+1)
// structuring element. Out-of-bounds pixels are treated as zero. Input is not
// modified.
//
// Two separable 1-D passes — O(W·H·r) work instead of O(W·H·r²).
// https://en.wikipedia.org/wiki/Dilation_(morphology)
func Dilate(img *common.SuzukiImage, radius int) (*common.SuzukiImage, error) {
        width, height := img.Width, img.Height
        tmp := common.NewSuzukiImage(width, height, img.HasPadding())
        out := common.NewSuzukiImage(width, height, img.HasPadding())

        for y := 0; y < height; y++ {
                for x := 0; x < width; x++ {
                        lo, hi := x-radius, x+radius
                        if lo < 0 {
                                lo = 0
                        }
                        if hi >= width {
                                hi = width - 1
                        }
                        for k := lo; k <= hi; k++ {
                                if img.GetXY(k, y) == 1 {
                                        tmp.SetXY(x, y, 1)
                                        break
                                }
                        }
                }
        }

        for y := 0; y < height; y++ {
                lo, hi := y-radius, y+radius
                if lo < 0 {
                        lo = 0
                }
                if hi >= height {
                        hi = height - 1
                }
                for x := 0; x < width; x++ {
                        for k := lo; k <= hi; k++ {
                                if tmp.GetXY(x, k) == 1 {
                                        out.SetXY(x, y, 1)
                                        break
                                }
                        }
                }
        }

        return out, nil
}