// 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"
	"fmt"
	"io"
	"log"
	"math"
	"math/bits"
	"math/rand"
	"os"
	"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)
	}
}

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

	type indexedCut struct {
		cut   octree.LineStringZ
		index int32
	}

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

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

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 {
				cuts[j] = append(cuts[j], indexedCut{cut, i / 3})
			}
		}
	}

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

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

		var pb polygonBuilder

		var isolated, inside int
	allInClass:
		for _, idx := range c {
			ns := neighbors(tri, idx)
			w := where(ns, c)
			switch bits.OnesCount8(w) {
			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
			}
		}

		log.Printf("\t%d: inside: %d / isolated: %d\n",
			i, inside, isolated)

		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
}

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 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 = 1000
	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))
}