// 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 octree

import (
	"bytes"
	"context"
	"database/sql"
	"encoding/binary"
	"fmt"
	"log"
	"math"
	"time"

	"github.com/tidwall/rtree"

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

type (
	ring []float64

	Polygon struct {
		rings   []ring
		indices []*rtree.RTree
	}

	IntersectionType byte

	lineSegment []float64
)

const (
	clippingPolygonSQL = `
WITH joined AS (
  SELECT
    sr.area      AS area,
    sr.date_info AS date_info
  FROM waterway.sounding_results sr
  WHERE sr.bottleneck_id = $1
)
SELECT ST_AsBinary(
  ST_intersection(
    (SELECT ST_Transform(area::geometry, $2::int) FROM joined WHERE date_info = $3::date),
    (SELECT ST_Transform(area::geometry, $2::int) FROM joined WHERE date_info = $4::date)
  )
  ) AS area
`
)

const (
	IntersectionInside IntersectionType = iota
	IntersectionOutSide
	IntersectionOverlaps
)

func LoadClippingPolygon(
	ctx context.Context,
	conn *sql.Conn,
	epsg uint32,
	bottleneck string,
	first, second time.Time,
) (*Polygon, error) {

	var clip []byte

	if err := conn.QueryRowContext(
		ctx, clippingPolygonSQL,
		bottleneck,
		epsg,
		first, second,
	).Scan(&clip); err != nil {
		return nil, err
	}

	var polygon Polygon
	if err := polygon.FromWKB(clip); err != nil {
		return nil, err
	}
	polygon.Indexify()
	return &polygon, nil
}

func (ls lineSegment) Rect(interface{}) ([]float64, []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[:]
}

func (p *Polygon) Indexify() {
	indices := make([]*rtree.RTree, len(p.rings))

	for i := range indices {
		index := rtree.New(nil)
		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]}
			}
			index.Insert(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
}

func (p *Polygon) IntersectionBox2D(box Box2D) IntersectionType {

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

	for _, index := range p.indices {
		var intersects bool
		index.Search(box, func(item rtree.Item) 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.
	point := []float64{box.X1, box.Y1}

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

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

func (p *Polygon) IntersectionWithTriangle(t *Triangle) IntersectionType {
	box := t.BBox()
	for _, index := range p.indices {
		var intersects bool
		index.Search(box, func(item rtree.Item) 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.
	point := []float64{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 hole.contains(point) {
				return IntersectionOutSide
			}
		}
	}

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

func (rng ring) isClosed() bool { return (len(rng) / 2) >= 3 }

func (rng ring) contains(point []float64) bool {
	if !rng.isClosed() {
		return false
	}

	end := len(rng)/2 - 1

	contains := intersectsWithRaycast(point, rng[:2], rng[end*2:end*2+2])

	for i := 2; i < len(rng); i += 2 {
		if intersectsWithRaycast(point, rng[i-2:i], rng[i:i+2]) {
			contains = !contains
		}
	}

	return contains
}

// 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(point, start, end []float64) bool {

	// Always ensure that the the first point
	// has a y coordinate that is less than the second point
	if start[1] > end[1] {
		// Switch the points if otherwise.
		start, end = end, start
	}

	// Move the point's y coordinate
	// outside of the bounds of the testing region
	// so we can start drawing a ray
	for point[1] == start[1] || point[1] == end[1] {
		y := math.Nextafter(point[1], math.Inf(1))
		point = []float64{point[0], y}
	}

	// If we are outside of the polygon, indicate so.
	if point[1] < start[1] || point[1] > end[1] {
		return false
	}

	if start[0] > end[0] {
		if point[0] > start[0] {
			return false
		}
		if point[0] < end[0] {
			return true
		}
	} else {
		if point[0] > end[0] {
			return false
		}
		if point[0] < start[0] {
			return true
		}
	}

	raySlope := (point[1] - start[1]) / (point[0] - start[0])
	diagSlope := (end[1] - start[1]) / (end[0] - start[0])

	return raySlope >= diagSlope
}

func (p *Polygon) NumVertices(ring int) int {
	if ring < 0 || ring >= len(p.rings) {
		return 0
	}
	return len(p.rings[ring]) / 2
}

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])
	}
}

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()
}

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.Printf("info: 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.Printf("info: 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
}

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}
}