// 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"
	"flag"
	"fmt"
	"io"
	"log"
	"math"
	"math/rand"
	"os"
	"runtime"
	"strconv"
	"strings"
	"sync"
	"time"

	svg "github.com/ajstarks/svgo"
	"github.com/fogleman/contourmap"

	"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 check(err error) {
	if err != nil {
		log.Fatalf("error: %v\n", err)
	}
}

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 dumpContoursSVG(
	w io.Writer,
	contours [][]contourmap.Contour,
) (err error) {

	var (
		minX = math.MaxFloat64
		minY = math.MaxFloat64
		maxX = -math.MaxFloat64
		maxY = -math.MaxFloat64
	)

	for _, c := range contours {
		for _, r := range c {
			for _, p := range r {
				if p.X < minX {
					minX = p.X
				}
				if p.Y < minY {
					minY = p.Y
				}
				if p.X > maxX {
					maxX = p.X
				}
				if p.Y > maxY {
					maxY = p.Y
				}
			}
		}
	}
	ratio := (maxX - minX) / (maxY - minY)

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

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

	px := linear(minX, 0, maxX, width)
	py := linear(minY, float64(height), maxY, 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 _, c := range contours {

		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 _, r := range c {
			x := make([]int, len(r))
			y := make([]int, len(r))
			for i, p := range r {
				x[i] = int(math.Round(px(p.X)))
				y[i] = int(math.Round(py(p.Y)))
			}

			canvas.Polygon(x, y, style)
		}
	}

	canvas.End()
	return
}

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 dumpArcsSVG(
	w io.Writer,
	min, max octree.Vertex,
	cuts [][]indexedArc,
	arcs []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 _, cut := range cuts {

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

		r := byte(rnd.Intn(256))
		g := byte(rnd.Intn(256))
		b := byte(rnd.Intn(256))

		style := fmt.Sprintf("fill:none;stroke:#%02x%02x%02x", r, g, b)

		for i := range cut {
			idx := cut[i].arc
			if _, already := usedArcs[idx]; already {
				continue
			}
			usedArcs[idx] = struct{}{}
			arc := arcs[idx]
			x := make([]int, len(arc))
			y := make([]int, len(arc))
			for i, p := range arc {
				x[i] = int(math.Round(px(p.X)))
				y[i] = int(math.Round(py(p.Y)))
			}

			canvas.Polyline(x, y, style)
		}
	}

	canvas.End()

	return nil
}

func main() {

	cellSize := float64(0.5)

	flag.Float64Var(&cellSize, "cellsize", 0.5, "width/height of raster cell [m]")

	flag.Parse()

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

	start = time.Now()
	tri, err := octree.Triangulate(xyz)
	check(err)
	log.Printf("triangulation took %v\n", time.Since(start))
	tooLongEdge := tri.EstimateTooLong()
	log.Printf("Eliminate triangles with edges longer than %.2f meters.", tooLongEdge)

	start = time.Now()
	polygon, removed := tri.ConcaveHull(tooLongEdge)

	polygonArea := polygon.Area()
	if polygonArea < 0.0 { // counter clockwise
		polygonArea = -polygonArea
		polygon.Reverse()
	}

	var clippingPolygon octree.Polygon
	clippingPolygon.FromLineStringZ(polygon)
	clippingPolygon.Indexify()

	tin := tri.Tin()
	// tin.EPSG = epsg

	var str octree.STRTree
	str.Build(tin)

	removed = str.Clip(&clippingPolygon)

	builder := octree.NewBuilder(tin)
	builder.Build(removed)

	tree := builder.Tree()

	min, max := tin.Min, tin.Max

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

	width := max.X - min.X
	height := max.Y - min.Y

	log.Printf("width/height: %.2f / %.2f\n", width, height)

	xcells := int(math.Ceil(width / cellSize))
	ycells := int(math.Ceil(height / cellSize))

	log.Printf("raster size: (%d, %d)\n", xcells, ycells)

	// Add border for closing
	raster := make([]float64, (xcells+2)*(ycells+2))

	const closed = -math.MaxFloat64
	for i := range raster {
		raster[i] = closed
	}

	start = time.Now()

	var wg sync.WaitGroup

	rows := make(chan int)

	rasterRow := func() {
		defer wg.Done()
		for i := range rows {
			pos := (i+1)*(xcells+2) + 1
			row := raster[pos : pos+xcells]
			//log.Printf("len: %d\n", len(row))
			py := min.Y + float64(i)*cellSize + cellSize/2
			px := min.X + cellSize/2
			for j := range row {
				v, ok := tree.Value(px, py)
				if ok {
					row[j] = v
				}
				px += cellSize
			}
		}
	}

	start = time.Now()

	for n := runtime.NumCPU(); n >= 1; n-- {
		wg.Add(1)
		go rasterRow()
	}

	for i := 0; i < ycells; i++ {
		rows <- i
	}
	close(rows)

	wg.Wait()

	log.Printf("Rasterizing took %v.\n", time.Since(start))

	start = time.Now()
	cm := contourmap.FromFloat64s(xcells+2, ycells+2, raster)

	start = time.Now()

	contours := make([][]contourmap.Contour, len(heights))

	cnts := make(chan int)

	doContours := func() {
		defer wg.Done()
		for i := range cnts {
			contours[i] = cm.Contours(heights[i])
		}
	}

	for n := runtime.NumCPU(); n >= 1; n-- {
		wg.Add(1)
		go doContours()
	}

	for i := range heights {
		cnts <- i
	}
	close(cnts)

	wg.Wait()

	log.Printf("Calculating contours took %v.\n", time.Since(start))
	check(dumpContoursSVG(os.Stdout, contours))

	/*

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

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