// 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) 2019 by via donau
//   – Österreichische Wasserstraßen-Gesellschaft mbH
// Software engineering by Intevation GmbH
//
// Author(s):
//  * Sascha L. Teichmann <sascha.teichmann@intevation.de>

package main

import (
	"bufio"
	"container/list"
	"fmt"
	"io"
	"log"
	"math"
	"math/bits"
	"math/rand"
	"os"
	"sort"
	"strconv"
	"strings"
	"time"

	svg "github.com/ajstarks/svgo"

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

const classBreaks = `1:#ff00dd,1.5,1.7,1.9,2.1,2.3,2.5:#f25f20,2.7,2.9,3.1:#f7e40e,3.3,3.5,4:#8ad51a,4.5,5,5.5,6,6.5,7:#1414ff`

func loadXYZ(r io.Reader) (octree.MultiPointZ, error) {

	scanner := bufio.NewScanner(r)

	points := make(octree.MultiPointZ, 0, 2000000)

	var x, y, z float64
	var err error

	for scanner.Scan() {
		line := strings.TrimSpace(scanner.Text())
		if len(line) == 0 || strings.HasPrefix(line, "#") {
			continue
		}
		parts := strings.SplitN(line, " ", 3)
		if len(parts) != 3 {
			continue
		}

		if x, err = strconv.ParseFloat(parts[0], 64); err != nil {
			return nil, err
		}
		if y, err = strconv.ParseFloat(parts[1], 64); err != nil {
			return nil, err
		}
		if z, err = strconv.ParseFloat(parts[2], 64); err != nil {
			return nil, err
		}
		points = append(points, octree.Vertex{X: x, Y: y, Z: z})
	}

	return points, nil
}

func minMax(points []octree.Vertex) (octree.Vertex, octree.Vertex) {
	if len(points) == 0 {
		return octree.Vertex{}, octree.Vertex{}
	}

	min, max := points[0], points[0]

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

	return min, max
}

func check(err error) {
	if err != nil {
		log.Fatalf("error: %v\n", err)
	}
}

type indexedArc struct {
	arc   int32
	index int32
}

func glue(a, b octree.LineStringZ) octree.LineStringZ {

	a1, a2 := a[0], a[len(a)-1]
	b1, b2 := b[0], b[len(b)-1]

	if a1.EpsEquals(b2) {
		c := make(octree.LineStringZ, len(a)-1+len(b))
		copy(c, b)
		copy(c[len(b):], a[1:])
		return c
	}

	if a2.EpsEquals(b1) {
		c := make(octree.LineStringZ, len(a)-1+len(b))
		copy(c, a)
		copy(c[len(a):], b[1:])
		return c
	}

	if a1.EpsEquals(b1) {
		c := make(octree.LineStringZ, len(a)-1+len(b))
		copy(c, b)
		c[:len(b)].Reverse()
		copy(c[len(b):], a[1:])
		return c
	}

	if a2.EpsEquals(b2) {
		c := make(octree.LineStringZ, len(a)-1+len(b))
		copy(c, a)
		copy(c[len(a):], b[:len(b)-1])
		c[len(a):].Reverse()
		return c
	}

	return nil
}

func connectArcs(tri *octree.Triangulation, cuts []indexedArc, arcs *[]octree.LineStringZ) {

	unique := map[int32]struct{}{}
	for i := range cuts {
		unique[cuts[i].arc] = struct{}{}
	}
	before := len(unique)

	origLen := int32(len(*arcs))

	for i := range cuts {
		if cuts[i].arc >= origLen {
			// already has a connected arc assigned.
			continue
		}

		line := (*arcs)[cuts[i].arc]

		var modified []int
		visited := map[int32]struct{}{}

		var recursive func(int32)
		recursive = func(idx int32) {
			visited[idx] = struct{}{}

			ns := neighbors(tri, idx)
			for _, n := range ns {
				n /= 3
				if _, already := visited[n]; already {
					continue
				}

				arcIdx := findArc(n, cuts)
				if arcIdx == -1 {
					continue
				}
				arc := cuts[arcIdx].arc
				if arc >= origLen {
					// already a new arc.
					continue
				}
				oline := (*arcs)[arc]

				nline := glue(line, oline)
				if nline == nil {
					// not joinable
					continue
				}
				line = nline
				modified = append(modified, arcIdx)
				recursive(cuts[arcIdx].index)
			}
		}
		recursive(cuts[i].index)
		if len(modified) > 0 {
			// alloc a new line an fix the references.
			nidx := int32(len(*arcs))
			*arcs = append(*arcs, line)
			cuts[i].arc = nidx
			for _, j := range modified {
				cuts[j].arc = nidx
			}
		}
	}

	unique = map[int32]struct{}{}
	for i := range cuts {
		unique[cuts[i].arc] = struct{}{}
	}
	log.Printf("unique arcs: before: %d after %d\n",
		before, len(unique))
}

func intersectTriangles(tri *octree.Triangulation, heights []float64) [][]octree.LineStringZ {

	cuts := make([][]indexedArc, len(heights))

	classes := make([][]int32, len(heights)+1)

	var arcs []octree.LineStringZ

all:
	for i := int32(0); i < int32(len(tri.Triangles)); i += 3 {
		ti := tri.Triangles[i : i+3]
		v0 := tri.Points[ti[0]]
		v1 := tri.Points[ti[1]]
		v2 := tri.Points[ti[2]]
		min := math.Min(v0.Z, math.Min(v1.Z, v2.Z))

		if min > heights[len(heights)-1] {
			classes[len(classes)-1] = append(classes[len(classes)-1], i/3)
			continue
		}
		max := math.Max(v0.Z, math.Max(v1.Z, v2.Z))

		for j, h := range heights {

			var l float64
			if j > 0 {
				l = heights[j-1]
			} else {
				l = -math.MaxFloat64
			}

			if l < min && max < h {
				classes[j] = append(classes[j], i/3)
				continue all
			}
			if min > h || max < h {
				continue
			}
			t := octree.Triangle{v0, v1, v2}
			cut := t.IntersectHorizontal(h)
			if len(cut) >= 2 {
				arc := int32(len(arcs))
				arcs = append(arcs, cut)
				cuts[j] = append(cuts[j], indexedArc{arc: arc, index: i / 3})
			}
		}
	}

	for l, c := range classes {
		sorted := sort.SliceIsSorted(c, func(i, j int) bool {
			return c[i] < c[j]
		})
		log.Printf("class[%d] sorted: %t\n", l, sorted)
	}

	for l, c := range cuts {
		sorted := sort.SliceIsSorted(c, func(i, j int) bool {
			return c[i].index < c[j].index
		})
		log.Printf("cut[%d] sorted: %t\n", l, sorted)
	}

	// connect the arcs in a cut list to longer arcs.

	for _, c := range cuts {
		connectArcs(tri, c, &arcs)
	}

	result := make([][]octree.LineStringZ, len(classes))

	log.Println("inside classes:")
	for i, c := range classes {

		pb := polygonBuilder{open: list.New()}

		usedArcs := map[int32]struct{}{}
		var dupes int

		var isolated, inside, found int
	allInClass:
		for _, idx := range c {
			ns := neighbors(tri, idx)
			mask := where(ns, c)
			switch bits.OnesCount8(mask) {
			case 0:
				// Totally insides do not contribute to the geometries.
				inside++
				continue allInClass
			case 3:
				// Isolated are areas w/o connections to the rest.
				isolated++
				ti := tri.Triangles[idx*3 : idx*3+3]
				pb.addTriangle(
					tri.Points[ti[0]],
					tri.Points[ti[1]],
					tri.Points[ti[2]])
				continue allInClass
			default:
				ti := tri.Triangles[idx*3 : idx*3+3]
				for j := 0; j < 3; j++ {
					if (mask & (1 << j)) == 0 {

						var curr octree.LineStringZ

						for l := i - 1; l <= i; l++ {
							if l < 0 || l >= len(cuts) {
								continue
							}
							arcIdx := findArc(ns[j]/3, cuts[l])
							if arcIdx == -1 {
								continue
							}
							found++
							aIdx := cuts[l][arcIdx].arc
							if _, already := usedArcs[aIdx]; already {
								dupes++
								continue
							}
							usedArcs[aIdx] = struct{}{}

							curr = arcs[aIdx]
							break
						}

						if curr == nil {
							k := (j + 1) % 3
							front := tri.Points[ti[j]]
							back := tri.Points[ti[k]]
							curr = octree.LineStringZ{front, back}
						}

					segment:
						for e := pb.open.Front(); e != nil; e = e.Next() {
							line := e.Value.(octree.LineStringZ)
							if nline := glue(curr, line); nline != nil {
								curr = nline
								pb.open.Remove(e)
								goto segment
							}
						} // all open

						// check if closed
						if curr[0].EpsEquals(curr[len(curr)-1]) {
							if !curr.CCW() {
								curr.Reverse()
							}
							pb.polygons = append(pb.polygons, curr)
						} else {
							pb.open.PushFront(curr)
						}
					} // for all border parts
				}
			}
		}

		log.Printf("\t%d: inside: %d / isolated: %d open: %d closed: %d dupes: %d found: %d\n",
			i, inside, isolated, pb.open.Len(), len(pb.polygons), dupes, found)

		/*
			for e := pb.open.Front(); e != nil; e = e.Next() {
				line := e.Value.(octree.LineStringZ)
				pb.polygons = append(pb.polygons, line)
			}
		*/

		result[i] = pb.polygons
	}

	log.Println("cuts:")
	for i, c := range cuts {
		log.Printf("\t%.3f: %d\n", heights[i], len(c))
	}

	return result

	// TODO: sew segments together.

}

type polygonBuilder struct {
	polygons []octree.LineStringZ

	open *list.List
}

func (pb *polygonBuilder) addTriangle(v0, v1, v2 octree.Vertex) {
	polygon := octree.LineStringZ{v0, v1, v2, v0}
	pb.polygons = append(pb.polygons, polygon)
}

func neighbors(t *octree.Triangulation, idx int32) []int32 {
	idx *= 3
	return t.Halfedges[idx : idx+3]
}

func findArc(needle int32, haystack []indexedArc) int {
	lo, hi := 0, len(haystack)-1
	for lo <= hi {
		mid := (hi-lo)/2 + lo
		switch v := haystack[mid].index; {
		case v < needle:
			lo = mid + 1
		case v > needle:
			hi = mid - 1
		default:
			return mid
		}
	}
	return -1
}

func contains(needle int32, haystack []int32) bool {
	lo, hi := 0, len(haystack)-1
	for lo <= hi {
		mid := (hi-lo)/2 + lo
		switch v := haystack[mid]; {
		case v < needle:
			lo = mid + 1
		case v > needle:
			hi = mid - 1
		default:
			return true
		}
	}
	return false
}

func where(neighbors, indices []int32) byte {
	var mask byte
	for i, n := range neighbors {
		if n < 0 || !contains(n/3, indices) {
			mask |= 1 << i
		}
	}
	return mask
}

func linear(x1, y1, x2, y2 float64) func(float64) float64 {
	// f(x1) = y1
	// f(x2) = y2
	// y1 = x1*a + b <=> b = y1 - x1*a
	// y2 = x2*a + b

	// y1 - y2 = a*(x1 - x2)
	// a = (y1-y2)/(x1 - x2) for x1 != x2

	if x1 == x2 {
		return func(float64) float64 {
			return 0.5 * (y1 + y2)
		}
	}
	a := (y1 - y2) / (x1 - x2)
	b := y1 - x1*a
	return func(x float64) float64 {
		return x*a + b
	}
}

func dumpPolygonsSVG(
	w io.Writer,
	min, max octree.Vertex,
	classPolygons [][]octree.LineStringZ,
) (err error) {

	ratio := (max.X - min.X) / (max.Y - min.Y)

	log.Printf("ratio: %.2f\n", ratio)

	const width = 50000
	height := int(math.Ceil(width * ratio))

	px := linear(min.X, 0, max.X, width)
	py := linear(min.Y, float64(height), max.Y, 0)

	out := bufio.NewWriter(w)
	defer func() { err = out.Flush() }()

	canvas := svg.New(out)
	canvas.Start(width, height)

	rnd := rand.New(rand.NewSource(42))

	for _, polygons := range classPolygons {

		r := byte(rnd.Intn(256))
		g := byte(rnd.Intn(256))
		b := byte(rnd.Intn(256))
		style := fmt.Sprintf("fill:#%02x%02x%02x", r, g, b)

		for _, polygon := range polygons {

			x := make([]int, len(polygon))
			y := make([]int, len(polygon))
			for i, p := range polygon {
				x[i] = int(math.Round(px(p.X)))
				y[i] = int(math.Round(py(p.Y)))
			}

			canvas.Polygon(x, y, style)
		}

	}

	canvas.End()

	return nil
}

func main() {

	heights, err := octree.ParseClassBreaks(classBreaks)
	check(err)
	log.Printf("num classes: %d\n", len(heights))

	start := time.Now()
	xyz, err := loadXYZ(os.Stdin)
	check(err)
	log.Printf("loading took %v\n", time.Since(start))

	log.Printf("num vertices: %d\n", len(xyz))

	min, max := minMax(xyz)

	log.Printf("(%.2f, %.2f, %.2f) - (%.2f, %.2f, %.2f)\n",
		min.X, min.Y, min.Z,
		max.X, max.Y, max.Z)

	start = time.Now()
	tri, err := octree.Triangulate(xyz)
	check(err)
	log.Printf("triangulation took %v\n", time.Since(start))

	log.Printf("number of triangles: %d\n", len(tri.Triangles)/3)

	start = time.Now()
	polygons := intersectTriangles(tri, heights)
	log.Printf("intersecting triangles took %v.\n", time.Since(start))

	check(dumpPolygonsSVG(os.Stdout, min, max, polygons))
}