// 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) 2018 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 (
	"log"
	"math"
	"runtime"
	"sync"
	"time"

	"github.com/fogleman/contourmap"

	"gemma.intevation.de/gemma/pkg/common"
	"gemma.intevation.de/gemma/pkg/wkb"
)

func GenerateRandomVertices(
	n int,
	min, max Vertex,
	eval func(float64, float64) (float64, bool),
	callback func([]Vertex),
) {
	var wg sync.WaitGroup

	jobs := make(chan int)
	out := make(chan []Vertex)
	done := make(chan struct{})

	cpus := runtime.NumCPU()

	free := make(chan []Vertex, cpus)

	for i := 0; i < cpus; i++ {
		wg.Add(1)
		go func() {
			defer wg.Done()
			xRange := common.Random(min.X, max.X)
			yRange := common.Random(min.Y, max.Y)

			for size := range jobs {
				var vertices []Vertex
				select {
				case vertices = <-free:
				default:
					vertices = make([]Vertex, 0, 1000)
				}
				for len(vertices) < size {
					x, y := xRange(), yRange()
					if z, ok := eval(x, y); ok {
						vertices = append(vertices, Vertex{X: x, Y: y, Z: z})
					}
				}
				out <- vertices
			}
		}()
	}

	go func() {
		defer close(done)
		for vertices := range out {
			callback(vertices)
			select {
			case free <- vertices[:0]:
			default:
			}
		}
	}()

	for remain := n; remain > 0; {
		if remain > 1000 {
			jobs <- 1000
			remain -= 1000
		} else {
			jobs <- remain
			remain = 0
		}
	}
	close(jobs)
	wg.Wait()
	close(out)
	<-done
}

func TraceAreas(
	heights []float64,
	cellSize float64,
	min, max Vertex,
	eval func(float64, float64) (float64, bool),
) []wkb.MultiPolygonGeom {

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

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

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

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

	start := time.Now()

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

	// prefill for no data
	const nodata = -math.MaxFloat64
	for i := range raster {
		raster[i] = nodata
	}

	// rasterize the height model

	var wg sync.WaitGroup

	rows := make(chan int)

	rasterRow := func() {
		defer wg.Done()
		quat := 0.25 * cellSize
		for i := range rows {
			pos := (i+1)*(xcells+2) + 1
			row := raster[pos : pos+xcells]
			py := min.Y + float64(i)*cellSize + cellSize/2
			px := min.X + cellSize/2
			y1 := py - quat
			y2 := py + quat
			for j := range row {
				var n int
				var sum float64

				if v, ok := eval(px-quat, y1); ok {
					sum = v
					n = 1
				}
				if v, ok := eval(px-quat, y2); ok {
					sum += v
					n++
				}
				if v, ok := eval(px+quat, y1); ok {
					sum += v
					n++
				}
				if v, ok := eval(px+quat, y2); ok {
					sum += v
					n++
				}

				if n > 0 {
					row[j] = sum / float64(n)
				}
				px += cellSize
			}
		}
	}

	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("info: Rastering took %v\n", time.Since(start))

	start = time.Now()

	tracer := contourmap.FromFloat64s(xcells+2, ycells+2, raster)

	areas := make([]wkb.MultiPolygonGeom, len(heights))

	// TODO: Check if this correct!
	reprojX := common.Linear(0.5, min.X, 1.5, min.X+cellSize)
	reprojY := common.Linear(0.5, min.Y, 1.5, min.Y+cellSize)

	cnts := make(chan int)

	doContours := func() {
		defer wg.Done()
		for hIdx := range cnts {
			c := tracer.Contours(heights[hIdx])

			if len(c) == 0 {
				continue
			}

			areas[hIdx] = buildMultipolygon(c, reprojX, reprojY)
		}
	}

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

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

	wg.Wait()
	log.Printf("info: Tracing areas took %v\n", time.Since(start))

	return areas
}

type contour []contourmap.Point

type bboxNode struct {
	box      Box2D
	cnt      contour
	children []*bboxNode
}

func (cnt contour) contains(o contour) bool {
	return contains(cnt, o[0].X, o[0].Y) ||
		contains(cnt, o[len(o)/2].X, o[len(o)/2].Y)
}

func (cnt contour) length() int {
	return len(cnt)
}

func (cnt contour) point(i int) (float64, float64) {
	return cnt[i].X, cnt[i].Y
}

func (cnt contour) bbox() Box2D {
	minX, minY := math.MaxFloat64, math.MaxFloat64
	maxX, maxY := -math.MaxFloat64, -math.MaxFloat64

	for _, p := range cnt {
		if p.X < minX {
			minX = p.X
		}
		if p.X > maxX {
			maxX = p.X
		}
		if p.Y < minY {
			minY = p.Y
		}
		if p.Y > maxY {
			maxY = p.Y
		}
	}
	return Box2D{X1: minX, X2: maxX, Y1: minY, Y2: maxY}
}

func (bn *bboxNode) insert(cnt contour, box Box2D) {
	// check if children are inside new
	var nr *bboxNode

	for i, r := range bn.children {
		if r.box.Inside(box) && cnt.contains(r.cnt) {
			if nr == nil {
				nr = &bboxNode{box: box, cnt: cnt}
			}
			nr.children = append(nr.children, r)
			bn.children[i] = nil
		}
	}

	// we have a new child
	if nr != nil {
		// compact the list
		for i := len(bn.children) - 1; i >= 0; i-- {
			if bn.children[i] == nil {
				if i < len(bn.children)-1 {
					copy(bn.children[i:], bn.children[i+1:])
				}
				bn.children[len(bn.children)-1] = nil
				bn.children = bn.children[:len(bn.children)-1]
			}
		}
		bn.children = append(bn.children, nr)
		return
	}

	// check if new is inside an old
	for _, r := range bn.children {
		if box.Inside(r.box) && r.cnt.contains(cnt) {
			r.insert(cnt, box)
			return
		}
	}

	// its a new child node.
	nr = &bboxNode{box: box, cnt: cnt}
	bn.children = append(bn.children, nr)
}

func (bn *bboxNode) insertRoot(cnt contour) {
	bn.insert(cnt, cnt.bbox())
}

type bboxOutFunc func(contour, []contour)

func (bn *bboxNode) generate(out bboxOutFunc) {

	var grands []*bboxNode

	holes := make([]contour, len(bn.children))

	for i, ch := range bn.children {
		holes[i] = ch.cnt
		grands = append(grands, ch.children...)
	}
	out(bn.cnt, holes)

	// the grand children are new polygons.
	for _, grand := range grands {
		grand.generate(out)
	}
}

func (bn *bboxNode) generateRoot(out bboxOutFunc) {
	for _, r := range bn.children {
		r.generate(out)
	}
}

func buildMultipolygon(
	cnts []contourmap.Contour,
	reprojX, reprojY func(float64) float64,
) wkb.MultiPolygonGeom {

	var forest bboxNode

	for _, cnt := range cnts {
		forest.insertRoot(contour(cnt))
	}

	//log.Printf("cnts: %d roots: %d\n", len(cnts), len(bf.roots))

	var mp wkb.MultiPolygonGeom

	out := func(sh contour, hls []contour) {

		polygon := make(wkb.PolygonGeom, 1+len(hls))

		// Handle shell
		shell := make(wkb.LinearRingGeom, len(sh))
		for j, pt := range sh {
			shell[j] = wkb.PointGeom{
				X: reprojX(pt.X),
				Y: reprojY(pt.Y),
			}
		}
		if shell.CCW() {
			shell.Reverse()
		}
		polygon[0] = shell

		// handle holes
		for i, hl := range hls {
			hole := make(wkb.LinearRingGeom, len(hl))
			for j, pt := range hl {
				hole[j] = wkb.PointGeom{
					X: reprojX(pt.X),
					Y: reprojY(pt.Y),
				}
			}
			if !hole.CCW() {
				hole.Reverse()
			}
			polygon[1+i] = hole
		}

		mp = append(mp, polygon)
	}

	forest.generateRoot(out)

	return mp
}