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