// This is Free Software under GNU Affero General Public License v >= 3.0
// without warranty, see README.md and license for details.
//
// SPDX-License-Identifier: AGPL-3.0-or-later
// License-Filename: LICENSES/AGPL-3.0.txt
//
// Copyright (C) 2018 by via donau
//   – Österreichische Wasserstraßen-Gesellschaft mbH
// Software engineering by Intevation GmbH
//
// Author(s):
//  * Sascha L. Teichmann <sascha.teichmann@intevation.de>

package mesh

import (
	"bytes"
	"encoding/binary"
	"fmt"
	"math"

	"github.com/tidwall/rtree"

	"gemma.intevation.de/gemma/pkg/log"
	"gemma.intevation.de/gemma/pkg/wkb"
)

type (
	ring []float64

	// Polygon has a border and holes.
	// The line segments are spatially indexed.
	Polygon struct {
		rings   []ring
		indices []*rtree.RTree
	}

	// IntersectionType represents an enum
	// of the type of intersection.
	IntersectionType byte

	lineSegment []float64
)

const (
	// IntersectionInside is inside the polygon.
	IntersectionInside IntersectionType = iota
	// IntersectionOutSide is outside the polygon.
	IntersectionOutSide
	// IntersectionOverlaps overlaps the polygon.
	IntersectionOverlaps
)

func (ls lineSegment) Rect() ([2]float64, [2]float64) {

	var min, max [2]float64

	if ls[0] < ls[2] {
		min[0] = ls[0]
		max[0] = ls[2]
	} else {
		min[0] = ls[2]
		max[0] = ls[0]
	}

	if ls[1] < ls[3] {
		min[1] = ls[1]
		max[1] = ls[3]
	} else {
		min[1] = ls[3]
		max[1] = ls[1]
	}

	return min, max
}

// Indexify creates a spatial index for thw polygon.
func (p *Polygon) Indexify() {
	indices := make([]*rtree.RTree, len(p.rings))

	for i := range indices {
		index := new(rtree.RTree)
		indices[i] = index

		rng := p.rings[i]

		for i := 0; i < len(rng); i += 2 {
			var ls lineSegment
			if i+4 <= len(rng) {
				ls = lineSegment(rng[i : i+4])
			} else {
				ls = []float64{rng[i], rng[i+1], rng[0], rng[1]}
			}
			min, max := ls.Rect()
			index.Insert(min, max, ls)
		}
	}

	p.indices = indices
}

func (ls lineSegment) intersects(a Box2D) bool {
	p1x := ls[0]
	p1y := ls[1]
	p2x := ls[2]
	p2y := ls[3]

	left := a.X1
	right := a.X2
	top := a.Y1
	bottom := a.Y2

	// The direction of the ray
	dx := p2x - p1x
	dy := p2y - p1y

	min, max := 0.0, 1.0

	var t0, t1 float64

	// Left and right sides.
	// - If the line is parallel to the y axis.
	if dx == 0 {
		if p1x < left || p1x > right {
			return false
		}
	} else {
		// - Make sure t0 holds the smaller value by checking the direction of the line.
		if dx > 0 {
			t0 = (left - p1x) / dx
			t1 = (right - p1x) / dx
		} else {
			t1 = (left - p1x) / dx
			t0 = (right - p1x) / dx
		}

		if t0 > min {
			min = t0
		}
		if t1 < max {
			max = t1
		}
		if min > max || max < 0 {
			return false
		}
	}

	// The top and bottom side.
	// - If the line is parallel to the x axis.
	if dy == 0 {
		if p1y < top || p1y > bottom {
			return false
		}
	} else {
		// - Make sure t0 holds the smaller value by checking the direction of the line.
		if dy > 0 {
			t0 = (top - p1y) / dy
			t1 = (bottom - p1y) / dy
		} else {
			t1 = (top - p1y) / dy
			t0 = (bottom - p1y) / dy
		}

		if t0 > min {
			min = t0
		}
		if t1 < max {
			max = t1
		}
		if min > max || max < 0 {
			return false
		}
	}

	// The point of intersection
	// ix = p1x + dx*min
	// iy = p1y + dy*min
	return true
}

func (ls lineSegment) intersectsLineSegment(o lineSegment) bool {

	p0 := ls[:2]
	p1 := ls[2:4]
	p2 := o[:2]
	p3 := o[2:4]

	s10x := p1[0] - p0[0]
	s10y := p1[1] - p0[1]
	s32x := p3[0] - p2[0]
	s32y := p3[1] - p2[1]

	den := s10x*s32y - s32x*s10y

	if den == 0 {
		return false
	}

	denPos := den > 0

	s02x := p0[0] - p2[0]
	s02y := p0[1] - p2[1]

	sNum := s10x*s02y - s10y*s02x
	if sNum < 0 == denPos {
		return false
	}

	tNum := s32x*s02y - s32y*s02x
	if tNum < 0 == denPos {
		return false
	}

	if sNum > den == denPos || tNum > den == denPos {
		return false
	}

	// t := tNum / den
	// intersection at( p0[0] + (t * s10x), p0[1] + (t * s10y) )
	return true
}

// IntersectionBox2D checks the type of intersection of the
// given box.
func (p *Polygon) IntersectionBox2D(box Box2D) IntersectionType {

	if len(p.rings) == 0 {
		return IntersectionOutSide
	}

	min, max := box.Rect()

	for _, index := range p.indices {
		var intersects bool
		index.Search(min, max, func(_, _ [2]float64, item any) bool {
			if item.(lineSegment).intersects(box) {
				intersects = true
				return false
			}
			return true
		})
		if intersects {
			return IntersectionOverlaps
		}
	}

	// No intersection -> check inside or outside
	// if an abritrary point  is inside or not.

	// Check holes first: inside a hole means outside.
	if len(p.rings) > 1 {
		for _, hole := range p.rings[1:] {
			if contains(hole, box.X1, box.Y1) {
				return IntersectionOutSide
			}
		}
	}

	// Check shell
	if contains(p.rings[0], box.X1, box.Y1) {
		return IntersectionInside
	}
	return IntersectionOutSide
}

// IntersectionWithTriangle checks the intersection type
// for the given triangle.
func (p *Polygon) IntersectionWithTriangle(t *Triangle) IntersectionType {
	box := t.BBox()
	min, max := box.Rect()
	for _, index := range p.indices {
		var intersects bool
		index.Search(min, max, func(_, _ [2]float64, item any) bool {
			ls := item.(lineSegment)
			other := make(lineSegment, 4)
			for i := range t {
				other[0] = t[i].X
				other[1] = t[i].Y
				other[2] = t[(i+1)%len(t)].X
				other[3] = t[(i+1)%len(t)].Y
				if ls.intersectsLineSegment(other) {
					intersects = true
					return false
				}
			}
			return true
		})
		if intersects {
			return IntersectionOverlaps
		}
	}
	// No intersection -> check inside or outside
	// if an abritrary point  is inside or not.
	pX, pY := t[0].X, t[0].Y

	// Check holes first: inside a hole means outside.
	if len(p.rings) > 1 {
		for _, hole := range p.rings[1:] {
			if contains(hole, pX, pY) {
				return IntersectionOutSide
			}
		}
	}

	// Check shell
	if contains(p.rings[0], pX, pY) {
		return IntersectionInside
	}
	return IntersectionOutSide
}

func (rng ring) length() int {
	return len(rng) / 2
}

func (rng ring) point(i int) (float64, float64) {
	i *= 2
	return rng[i], rng[i+1]
}

type segments interface {
	length() int
	point(int) (float64, float64)
}

func contains(s segments, pX, pY float64) bool {

	n := s.length()
	if n < 3 {
		return false
	}

	sX, sY := s.point(0)
	eX, eY := s.point(n - 1)

	const eps = 0.0000001

	if math.Abs(sX-eX) > eps || math.Abs(sY-eY) > eps {
		// It's not closed!
		return false
	}

	var inside bool

	for i := 1; i < n; i++ {
		eX, eY := s.point(i)
		if intersectsWithRaycast(pX, pY, sX, sY, eX, eY) {
			inside = !inside
		}
		sX, sY = eX, eY
	}

	return inside
}

// Using the raycast algorithm, this returns whether or not the passed in point
// Intersects with the edge drawn by the passed in start and end points.
// Original implementation: http://rosettacode.org/wiki/Ray-casting_algorithm#Go
func intersectsWithRaycast(pX, pY, sX, sY, eX, eY float64) bool {

	// Always ensure that the the first point
	// has a y coordinate that is less than the second point
	if sY > eY {
		// Switch the points if otherwise.
		sX, sY, eX, eY = eX, eY, sX, sY
	}

	// Move the point's y coordinate
	// outside of the bounds of the testing region
	// so we can start drawing a ray
	for pY == sY || pY == eY {
		pY = math.Nextafter(pY, math.Inf(1))
	}

	// If we are outside of the polygon, indicate so.
	if pY < sY || pY > eY {
		return false
	}

	if sX > eX {
		if pX > sX {
			return false
		}
		if pX < eX {
			return true
		}
	} else {
		if pX > eX {
			return false
		}
		if pX < sX {
			return true
		}
	}

	raySlope := (pY - sY) / (pX - sX)
	diagSlope := (eY - sY) / (eX - sX)

	return raySlope >= diagSlope
}

// NumVertices returns the number of vertices of a given ring.
func (p *Polygon) NumVertices(ring int) int {
	if ring < 0 || ring >= len(p.rings) {
		return 0
	}
	return len(p.rings[ring]) / 2
}

// Vertices passes the vertices of a given ring
// to the given fn function.
func (p *Polygon) Vertices(ring int, fn func(float64, float64)) {
	if ring < 0 || ring >= len(p.rings) {
		return
	}
	rng := p.rings[ring]
	for i := 0; i < len(rng); i += 2 {
		fn(rng[i+0], rng[i+1])
	}
}

// AsWKB serializes the polygon as WKB.
func (p *Polygon) AsWKB() []byte {

	size := 1 + 4 + 4
	for _, r := range p.rings {
		size += 4 + len(r)*8
	}

	buf := bytes.NewBuffer(make([]byte, 0, size))

	binary.Write(buf, binary.LittleEndian, wkb.NDR)
	binary.Write(buf, binary.LittleEndian, wkb.Polygon)
	binary.Write(buf, binary.LittleEndian, uint32(len(p.rings)))

	for _, r := range p.rings {
		binary.Write(buf, binary.LittleEndian, uint32(len(r)/2))
		for i := 0; i < len(r); i += 2 {
			binary.Write(buf, binary.LittleEndian, math.Float64bits(r[i+0]))
			binary.Write(buf, binary.LittleEndian, math.Float64bits(r[i+1]))
		}
	}

	return buf.Bytes()
}

// FromWKB deserializes a polygon from WKB.
func (p *Polygon) FromWKB(data []byte) error {

	r := bytes.NewReader(data)

	endian, err := r.ReadByte()

	var order binary.ByteOrder

	switch {
	case err != nil:
		return err
	case endian == wkb.NDR:
		order = binary.LittleEndian
	case endian == wkb.XDR:
		order = binary.BigEndian
	default:
		return fmt.Errorf("unknown byte order %x", endian)
	}

	var geomType uint32
	err = binary.Read(r, order, &geomType)

	switch {
	case err != nil:
		return err
	case geomType != wkb.Polygon:
		return fmt.Errorf("unknown geometry type %x", geomType)
	}

	var numRings uint32
	if err = binary.Read(r, order, &numRings); err != nil {
		return err
	}

	rngs := make([]ring, numRings)

	log.Infof("number of rings: %d\n", len(rngs))

	for rng := uint32(0); rng < numRings; rng++ {
		var numVertices uint32
		if err = binary.Read(r, order, &numVertices); err != nil {
			return err
		}

		log.Infof("number of vertices in ring %d: %d\n", rng, numVertices)

		numVertices *= 2
		vertices := make([]float64, numVertices)

		for v := uint32(0); v < numVertices; v += 2 {
			var lat, lon uint64
			if err = binary.Read(r, order, &lat); err != nil {
				return err
			}
			if err = binary.Read(r, order, &lon); err != nil {
				return err
			}
			vertices[v] = math.Float64frombits(lat)
			vertices[v+1] = math.Float64frombits(lon)
		}

		rngs[rng] = vertices
	}

	p.rings = rngs

	return nil
}

// FromLineStringZ creates a polygon from a given linestring z.
func (p *Polygon) FromLineStringZ(ls LineStringZ) {
	r := make([]float64, 2*len(ls))
	var pos int
	for i := range ls {
		r[pos+0] = ls[i].X
		r[pos+1] = ls[i].Y
		pos += 2
	}
	p.rings = []ring{r}
}