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

import (
	"math"
	"runtime"
	"sync"
	"time"

	"github.com/fogleman/contourmap"

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

// Raster represents a 2D cell raster with a given extend.
// The cells are equally sized in X and Y direction.
type Raster struct {
	BBox     Box2D
	CellSize float64
	XCells   int
	YCells   int
	Cells    []float64
}

const noData = -math.MaxFloat64

// NewRaster creates a new Raster with the given extend
// and cell size.
func NewRaster(bbox Box2D, cellSize float64) *Raster {

	width, height := bbox.Size()

	log.Infof("raster extent: %.2f / %.2f", width, height)

	xCells := int(math.Ceil(width / cellSize))
	yCells := int(math.Ceil(height / cellSize))

	log.Infof("raster size: %d / %d\n", xCells, yCells)

	size := (xCells + 2) * (yCells + 2)
	cells := make([]float64, size)
	for i := range cells {
		cells[i] = noData
	}
	return &Raster{
		BBox:     bbox,
		CellSize: cellSize,
		XCells:   xCells,
		YCells:   yCells,
		Cells:    cells,
	}
}

// Rasterize fills the raster with the evaluation data
// received via the given eval function. The cell values
// 4x oversampled,
func (r *Raster) Rasterize(eval func(float64, float64) (float64, bool)) {
	var wg sync.WaitGroup

	rows := make(chan int)

	rasterRow := func() {
		defer wg.Done()
		quat := 0.25 * r.CellSize
		for i := range rows {
			pos := (i+1)*(r.XCells+2) + 1
			row := r.Cells[pos : pos+r.XCells]
			py := r.BBox.Y1 + float64(i)*r.CellSize + r.CellSize/2
			px := r.BBox.X1 + r.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 && !math.IsNaN(sum) {
					row[j] = sum / float64(n)
				} else {
					row[j] = noData
				}
				px += r.CellSize
			}
		}
	}

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

	for i := 0; i < r.YCells; i++ {
		rows <- i
	}
	close(rows)
	wg.Wait()
}

// Diff updates the cell values of the raster with difference
// of this raster and the values from the gien eval function.
// The data will 4x over-sampled.
func (r *Raster) Diff(eval func(float64, float64) (float64, bool)) {
	var wg sync.WaitGroup

	rows := make(chan int)

	rasterRow := func() {
		defer wg.Done()
		quat := 0.25 * r.CellSize
		for i := range rows {
			pos := (i+1)*(r.XCells+2) + 1
			row := r.Cells[pos : pos+r.XCells]
			py := r.BBox.Y1 + float64(i)*r.CellSize + r.CellSize/2
			px := r.BBox.X1 + r.CellSize/2
			y1 := py - quat
			y2 := py + quat
			for j, old := range row {
				// only diff where values
				if old == noData {
					px += r.CellSize
					continue
				}
				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 && !math.IsNaN(sum) {
					row[j] -= sum / float64(n)
				} else {
					row[j] = noData
				}

				px += r.CellSize
			}
		}
	}

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

	for i := 0; i < r.YCells; i++ {
		rows <- i
	}
	close(rows)
	wg.Wait()
}

// ZExtent returns the z-range of the raster.
// The last return value is false if the raster
// only consists of no data values.
func (r *Raster) ZExtent() (float64, float64, bool) {
	min, max := math.MaxFloat64, -math.MaxFloat64
	for _, v := range r.Cells {
		if v == noData {
			continue
		}
		if v < min {
			min = v
		}
		if v > max {
			max = v
		}
	}
	return min, max, min != math.MaxFloat64
}

// Trace generates contour lines of this raster given
// the given class breaks.
func (r *Raster) Trace(heights ClassBreaks) []wkb.MultiPolygonGeom {
	start := time.Now()

	tracer := contourmap.FromFloat64s(r.XCells+2, r.YCells+2, r.Cells)

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

	reprojX := common.Linear(0.5, r.BBox.X1, 1.5, r.BBox.X1+r.CellSize)
	reprojY := common.Linear(0.5, r.BBox.Y1, 1.5, r.BBox.Y1+r.CellSize)

	var wg sync.WaitGroup

	cnts := make(chan int)

	doContours := func() {
		defer wg.Done()
		for hIdx := range cnts {
			if c := tracer.Contours(heights[hIdx]); len(c) > 0 {
				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.Infof("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.Debugf("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
}