// Copyright (C) 2018 Michael Fogleman
//
// Permission is hereby granted, free of charge, to any person obtaining
// a copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

package mesh

import (
	"math"

	"gemma.intevation.de/gemma/pkg/log"
	"gonum.org/v1/gonum/stat"
)

// Triangulation represents a triangulation
// consisting of its vertices, the convex hull,
// the triangles and a half edge representation.
type Triangulation struct {
	Points     []Vertex
	ConvexHull []Vertex
	Triangles  []int32
	Halfedges  []int32
}

// Triangulate returns a Delaunay triangulation of the provided points.
func Triangulate(points []Vertex) (*Triangulation, error) {
	t := newTriangulator(points)
	err := t.triangulate()
	return &Triangulation{points, t.convexHull(), t.triangles, t.halfedges}, err
}

// EstimateTooLong estimates an edge length which is too long
// to be a real triangle in the triangulation.
// It is assumed that the triangles in this mesh are
// more or less equally sized. If an edge length is more
// than the 3.5 of the standard dev of the medium length
// it is considered as too long.
func (t *Triangulation) EstimateTooLong() float64 {

	num := len(t.Triangles) / 3

	lengths := make([]float64, 0, num)

	points := t.Points

tris:
	for i := 0; i < num; i++ {
		idx := i * 3
		var max float64
		vs := t.Triangles[idx : idx+3]
		for j, vj := range vs {
			if t.Halfedges[idx+j] < 0 {
				continue tris
			}
			k := (j + 1) % 3
			if l := points[vj].Distance2D(points[vs[k]]); l > max {
				max = l
			}
		}
		lengths = append(lengths, max)
	}

	std := stat.StdDev(lengths, nil)
	return 3.5 * std
}

// ConcaveHull constructs a concave hull for this mesh,
// Triangles that are considered as too large based on
// the EstimateTooLong estimation are removed from the border.
func (t *Triangulation) ConcaveHull(tooLong float64) (LineStringZ, map[int32]struct{}) {

	if tooLong <= 0 {
		tooLong = t.EstimateTooLong()
	}

	tooLong *= tooLong

	var candidates []int32

	points := t.Points

	for i, num := 0, len(t.Triangles)/3; i < num; i++ {
		idx := i * 3
		var max float64
		vs := t.Triangles[idx : idx+3]
		for j, vj := range vs {
			k := (j + 1) % 3
			if l := points[vj].SquaredDistance2D(points[vs[k]]); l > max {
				max = l
			}
		}
		if max > tooLong {
			candidates = append(candidates, int32(i))
		}
	}

	removed := map[int32]struct{}{}

	isBorder := func(n int32) bool {
		n *= 3
		for i := int32(0); i < 3; i++ {
			e := n + i
			o := t.Halfedges[e]
			if o < 0 {
				return true
			}
			if _, found := removed[o/3]; found {
				return true
			}
		}
		return false
	}

	var newCandidates []int32

	log.Infof("candidates: %d\n", len(candidates))
	for len(candidates) > 0 {

		oldRemoved := len(removed)

		for _, i := range candidates {

			if isBorder(i) {
				removed[i] = struct{}{}
			} else {
				newCandidates = append(newCandidates, i)
			}
		}

		if oldRemoved == len(removed) {
			break
		}

		candidates = newCandidates
		newCandidates = newCandidates[:0]
	}

	log.Infof("candidates left: %d\n", len(candidates))
	log.Infof("triangles: %d\n", len(t.Triangles)/3)
	log.Infof("info: triangles to remove: %d\n", len(removed))

	type edge struct {
		a, b       int32
		prev, next *edge
	}

	isClosed := func(e *edge) bool {
		for curr := e.next; curr != nil; curr = curr.next {
			if curr == e {
				return true
			}
		}
		return false
	}

	open := map[int32]*edge{}
	var rings []*edge

	for i, num := int32(0), int32(len(t.Triangles)/3); i < num; i++ {
		if _, found := removed[i]; found {
			continue
		}
		n := i * 3
		for j := int32(0); j < 3; j++ {
			e := n + j
			f := t.Halfedges[e]
			if f >= 0 {
				if _, found := removed[f/3]; !found {
					continue
				}
			}
			a := t.Triangles[e]
			b := t.Triangles[n+(j+1)%3]

			curr := &edge{a: a, b: b}

			if old := open[a]; old != nil {
				delete(open, a)
				if old.a == a {
					old.prev = curr
					curr.next = old
				} else {
					old.next = curr
					curr.prev = old
				}

				if isClosed(old) {
					rings = append(rings, old)
				}
			} else {
				open[a] = curr
			}

			if old := open[b]; old != nil {
				delete(open, b)
				if old.b == b {
					old.next = curr
					curr.prev = old
				} else {
					old.prev = curr
					curr.next = old
				}

				if isClosed(old) {
					rings = append(rings, old)
				}
			} else {
				open[b] = curr
			}
		}
	}

	if len(open) > 0 {
		log.Warnf("open vertices left: %d\n", len(open))
	}

	if len(rings) == 0 {
		log.Warnln("no ring found")
		return nil, removed
	}

	curr := rings[0]

	polygon := LineStringZ{
		points[curr.a],
		points[curr.b],
	}

	for curr = curr.next; curr != rings[0]; curr = curr.next {
		polygon = append(polygon, points[curr.b])
	}

	polygon = append(polygon, t.Points[rings[0].a])

	log.Infof("length of boundary: %d\n", len(polygon))

	return polygon, removed
}

// TriangleSlices returns the indices of the triangles
// of the mesh as slices og length three.
func (t *Triangulation) TriangleSlices() [][]int32 {
	sl := make([][]int32, len(t.Triangles)/3)
	var j int
	for i := range sl {
		sl[i] = t.Triangles[j : j+3]
		j += 3
	}
	return sl
}

// Tin creates a TIN out of this triangulation.
func (t *Triangulation) Tin() *Tin {

	min := Vertex{math.MaxFloat64, math.MaxFloat64, math.MaxFloat64}
	max := Vertex{-math.MaxFloat64, -math.MaxFloat64, -math.MaxFloat64}

	vertices := t.Points

	for _, v := range vertices {
		min.Minimize(v)
		max.Maximize(v)
	}

	return &Tin{
		Vertices:  vertices,
		Triangles: t.TriangleSlices(),
		Min:       min,
		Max:       max,
	}
}