// 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 octree

import (
	"math"
)

func handleTriangle(
	t *Triangle,
	maxDist, tolerance float64,
	maxIdx int,
	points MultiPointZ,
	result *MultiPointZ,
) bool {
	if maxDist <= tolerance {
		return false
	}

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

	var (
		tris     [3]Triangle
		planes   [3]Plane3D
		maxDists [3]float64
		maxIdxs  [3]int
		parts    [3]MultiPointZ
	)

	top := points[maxIdx]
	for i := 0; i < 3; i++ {
		tris[i] = 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], tolerance,
			maxIdxs[i],
			part,
			result,
		) {
			found = true
		}
	}

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

	return found
}

func (points MultiPointZ) Simplify(tolerance float64) MultiPointZ {

	if len(points) < 2 {
		return points
	}

	if tolerance < 0 {
		tolerance = -tolerance
	}

	min := Vertex{X: math.MaxFloat64, Y: math.MaxFloat64, Z: math.MaxFloat64}
	max := 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 := []Vertex{
		{xMin, yMin, below},
		{xMax, yMin, below},
		{xMax, yMax, below},
		{xMin, yMax, below},
	}

	top := points[maxIdx]

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

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

	parts := make([][]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(MultiPointZ, 0, len(points))

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

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

	return result
}