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

			// We need to bring it back to the
			// none raster coordinate system.
			a := make(wkb.MultiPolygonGeom, len(c))

			for i, pl := range c {
				shell := make(wkb.LinearRingGeom, len(pl))
				for j, pt := range pl {
					shell[j] = wkb.PointGeom{
						X: reprojX(pt.X),
						Y: reprojY(pt.Y),
					}
				}
				a[i] = wkb.PolygonGeom{shell}
			}

			areas[hIdx] = a
		}
	}

	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
}