// 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"
	"os"
	"strconv"
	"strings"

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

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 storeXYZ(points octree.MultiPointZ, w io.Writer) error {
	out := bufio.NewWriter(w)
	for i := range points {
		fmt.Fprintf(out, "%.5f,%.5f,%.5f\n",
			points[i].X, points[i].Y, points[i].Z)
	}
	return out.Flush()
}

func simplify(points octree.MultiPointZ, tolerance float64) octree.MultiPointZ {

	if len(points) < 2 {
		return points
	}

	if tolerance < 0 {
		tolerance = -tolerance
	}

	min := octree.Vertex{X: math.MaxFloat64, Y: math.MaxFloat64, Z: math.MaxFloat64}
	max := octree.Vertex{X: -math.MaxFloat64, Y: -math.MaxFloat64, Z: -math.MaxFloat64}

	var maxIdx int

	for i, v := range points {
		min.Minimize(v)

		if v.X < min.X {
			min.X = v.X
		}
		if v.X > max.X {
			max.X = v.X
		}
		if v.Y < min.Y {
			min.Y = v.Y
		}
		if v.Y > max.Y {
			max.Y = v.Y
		}
		if v.Z < min.Z {
			min.Z = v.Z
		}
		if v.Z > max.Z {
			max.Z = v.Z
			maxIdx = i
		}

		max.Maximize(v)
	}

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

	below := min.Z - 3*tolerance
	xMin := min.X - tolerance
	xMax := max.X + tolerance
	yMin := min.Y - tolerance
	yMax := max.Y + tolerance

	corners := []octree.Vertex{
		{xMin, yMin, below},
		{xMax, yMin, below},
		{xMax, yMax, below},
		{xMin, yMax, below},
	}

	top := points[maxIdx]

	tris := make([]octree.Triangle, len(corners))
	planes := make([]octree.Plane3D, len(corners))

	for i, v1 := range corners {
		v2 := corners[(i+1)%len(corners)]
		tris[i] = octree.Triangle{v1, v2, top}
		planes[i] = tris[i].Plane3D()
	}

	parts := make([][]octree.Vertex, len(tris))

	maxDists := make([]float64, len(planes))
	maxIdxs := make([]int, len(planes))

nextPoint:
	for i, v := range points {
		if i == maxIdx {
			continue
		}

		for j := range tris {
			if tris[j].Contains(v.X, v.Y) {
				if dist := math.Abs(planes[j].Eval(v)); dist > maxDists[j] {
					maxDists[j] = dist
					maxIdxs[j] = len(parts[j])
				}
				parts[j] = append(parts[j], v)
				continue nextPoint
			}
		}
	}

	result := make([]octree.Vertex, 0, len(points))

	var handleTriangle func(*octree.Triangle, float64, int, []octree.Vertex) bool

	handleTriangle = func(
		t *octree.Triangle,
		maxDist float64, maxIdx int,
		points []octree.Vertex,
	) bool {
		if maxDist <= tolerance {
			return false
		}

		if len(points) == 1 {
			result = append(result, points[0])
		}

		var (
			tris     [3]octree.Triangle
			planes   [3]octree.Plane3D
			maxDists [3]float64
			maxIdxs  [3]int
			parts    [3][]octree.Vertex
		)

		top := points[maxIdx]
		for i := 0; i < 3; i++ {
			tris[i] = octree.Triangle{t[i], t[(i+1)%3], top}
			planes[i] = tris[i].Plane3D()
		}

	nextPoint:
		for i, v := range points {
			if i == maxIdx {
				continue
			}

			for j := range tris {
				if tris[j].Contains(v.X, v.Y) {
					if dist := math.Abs(planes[j].Eval(v)); dist > maxDists[j] {
						maxDists[j] = dist
						maxIdxs[j] = len(parts[j])
					}
					parts[j] = append(parts[j], v)
					continue nextPoint
				}
			}
		}

		var found bool
		for i, part := range parts {
			if len(part) > 0 && handleTriangle(&tris[i], maxDists[i], maxIdxs[i], part) {
				found = true
			}
		}

		if found {
			result = append(result, top)
		}

		return found
	}

	var found bool
	for i, part := range parts {
		if len(part) > 0 && handleTriangle(&tris[i], maxDists[i], maxIdxs[i], part) {
			found = true
		}
	}

	if found {
		result = append(result, top)
	}

	return result
}

func main() {

	var tolerance float64

	flag.Float64Var(&tolerance, "t", 0.1, "accepted tolerance (shorthand)")
	flag.Float64Var(&tolerance, "tolerance", 0.1, "accepted tolerance")

	flag.Parse()

	points, err := loadXYZ(os.Stdin)
	if err != nil {
		log.Fatalf("err: %v\n", err)
	}

	points = simplify(points, tolerance)

	if err := storeXYZ(points, os.Stdout); err != nil {
		log.Fatalf("err: %v\n", err)
	}
}