// 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 (
	"fmt"
	"log"
	"math"
	"os"
)

// Tree is an Octree holding triangles.
type Tree struct {
	// EPSG is the projection.
	EPSG uint32

	vertices  []Vertex
	triangles [][]int32
	index     []int32

	// Min is the lower left corner of the bbox.
	Min Vertex
	// Max is the upper right corner of the bbox.
	Max Vertex
}

type boxFrame struct {
	pos int32
	Box2D
}

func (ot *Tree) Vertices() []Vertex {
	return ot.vertices
}

var scale = [4][4]float64{
	{0.0, 0.0, 0.5, 0.5},
	{0.5, 0.0, 1.0, 0.5},
	{0.0, 0.5, 0.5, 1.0},
	{0.5, 0.5, 1.0, 1.0},
}

func (ot *Tree) Dump() error {
	fake, err := os.Create("/tmp/dump.geojson")
	if err != nil {
		return err
	}
	defer fake.Close()

	fmt.Fprintln(fake,
		`{
"crs": {
    "type": "name",
    "properties": {
      "name": "urn:ogc:def:crs:EPSG::32633"
    }
  },

  "type": "FeatureCollection",
  "features": [`)

	first := true

	already := make(map[int32]struct{})

	stack := []int32{1}
	for len(stack) > 0 {
		pos := stack[len(stack)-1]
		stack = stack[:len(stack)-1]

		if pos > 0 { // node
			index := ot.index[pos:]
			if len(index) > 8 {
				index = index[:8]
			}
			for i := range index {
				if index[i] != 0 {
					stack = append(stack, index[i])
				}
			}
		} else {
			pos = -pos - 1
			n := ot.index[pos]
			indices := ot.index[pos+1 : pos+1+n]
			for _, index := range indices {
				if _, found := already[index]; found {
					continue
				}
				already[index] = struct{}{}
				tri := ot.triangles[index]
				t := Triangle{
					ot.vertices[tri[0]],
					ot.vertices[tri[1]],
					ot.vertices[tri[2]],
				}
				if first {
					first = false
				} else {
					fmt.Fprintln(fake, ",")
				}
				fmt.Fprintf(fake,
					`{ "type": "Feature",
  "properties": {
     "id": %d
  },
  "geometry": {
     "type": "Polygon",
     "coordinates": [[
       [ %.6f, %.6f ],
       [ %.6f, %.6f ],
       [ %.6f, %.6f ],
       [ %.6f, %.6f ]]
    ]
  }
}`, index,
					t[0].X, t[0].Y,
					t[1].X, t[1].Y,
					t[2].X, t[2].Y,
					t[0].X, t[0].Y)
			}
		}
	}

	fmt.Fprintln(fake,
		`]
}`)
	return nil
}

func (ot *Tree) eliminate(removed map[int32]IntersectionType) {

	var bad int

	stack := []int32{1}
	for len(stack) > 0 {
		pos := stack[len(stack)-1]
		stack = stack[:len(stack)-1]

		if pos > 0 { // node
			index := ot.index[pos:]
			if len(index) > 8 {
				index = index[:8]
			}
			for i := range index {
				if index[i] != 0 {
					stack = append(stack, index[i])
				}
			}
		} else {
			pos = -pos - 1
			n := ot.index[pos]
			indices := ot.index[pos+1 : pos+1+n]

			for i := len(indices) - 1; i >= 0; i-- {
				index := indices[i]
				if t, found := removed[index]; found && t != IntersectionInside {
					if i < len(indices)-1 {
						copy(indices[i:], indices[i+1:])
					}
					indices[len(indices)-1] = 0
					indices = indices[:len(indices)-1]
					bad++
				}
			}
			ot.index[pos] = int32(len(indices))
		}
	}

	log.Printf("bad triangles: %d\n", bad)
}

func (ot *Tree) Clip(p *Polygon) {

	/*
			fake, _ := os.Create("/tmp/removed.geojson")
			defer fake.Close()

			fmt.Fprintln(fake,
				`{
		"crs": {
		    "type": "name",
		    "properties": {
		      "name": "urn:ogc:def:crs:EPSG::32633"
		    }
		  },

		  "type": "FeatureCollection",
		  "features": [`)

			first := true
	*/

	log.Printf("info: num triangles: %d\n", len(ot.triangles))

	all := Box2D{ot.Min.X, ot.Min.Y, ot.Max.X, ot.Max.Y}

	stack := []boxFrame{{1, all}}

	triChecks := make(map[int32]IntersectionType)

	var nodeTests int
	var nodesClipped int
	var trianglesClipped int
	var nodesAllInside int

frames:
	for len(stack) > 0 {
		top := stack[len(stack)-1]
		stack = stack[:len(stack)-1]

		if top.pos > 0 { // node
			nodeTests++
			switch p.IntersectionBox2D(top.Box2D) {
			case IntersectionInside:
				// all inside so nothing to clip.
				nodesAllInside++
				continue frames
			case IntersectionOutSide:
				// all outside -> clip from tree.
				index := ot.index[top.pos:]
				if len(index) > 8 {
					index = index[:8]
				}
				for i := range index {
					if index[i] != 0 {
						index[i] = 0
						nodesClipped++
					}
				}
				continue frames
			default: // Overlaps
				if index := ot.index[top.pos:]; len(index) > 7 {
				children:
					for i := 0; i < 4; i++ {
						a := index[i]
						b := index[i+4]
						if a == 0 && b == 0 {
							continue
						}
						dx := top.X2 - top.X1
						dy := top.Y2 - top.Y1
						nbox := Box2D{
							dx*scale[i][0] + top.X1,
							dy*scale[i][1] + top.Y1,
							dx*scale[i][2] + top.X1,
							dy*scale[i][3] + top.Y1,
						}
						switch p.IntersectionBox2D(nbox) {
						case IntersectionInside:
							// all inside so nothing to clip.
							nodesAllInside++
							continue children
						case IntersectionOutSide:
							// all are ouside -> clip from tree.
							if a != 0 {
								index[i] = 0
								nodesClipped++
							}
							if b != 0 {
								index[i+4] = 0
								nodesClipped++
							}
							continue children
						default: // Overlaps
							if a != 0 {
								stack = append(stack, boxFrame{a, nbox})
							}
							if b != 0 {
								stack = append(stack, boxFrame{b, nbox})
							}
						}
					}
				}
			}
		} else { // leaf
			pos := -top.pos - 1
			n := ot.index[pos]
			indices := ot.index[pos+1 : pos+1+n]
		tris:
			for i := len(indices) - 1; i >= 0; i-- {
				triIndex := indices[i]
				what, found := triChecks[triIndex]
				if !found {
					tri := ot.triangles[triIndex]
					t := Triangle{
						ot.vertices[tri[0]],
						ot.vertices[tri[1]],
						ot.vertices[tri[2]],
					}
					what = p.IntersectionWithTriangle(&t)
					triChecks[triIndex] = what

					/*

											if what != IntersectionInside {
												if first {
													first = false
												} else {
													fmt.Fprintln(fake, ",")
												}
												fmt.Fprintf(fake,
													`{ "type": "Feature",
						  "properties": {
						    "id": %d
						  },
						  "geometry": {
						     "type": "Polygon",
						     "coordinates": [[
						       [ %.6f, %.6f ],
						       [ %.6f, %.6f ],
						       [ %.6f, %.6f ],
						       [ %.6f, %.6f ]]
						    ]
						  }
						}`, triIndex,
													t[0].X, t[0].Y,
													t[1].X, t[1].Y,
													t[2].X, t[2].Y,
													t[0].X, t[0].Y)
											}
					*/
				}
				switch what {
				case IntersectionInside:
					// triangle inside -> stay.
					continue tris
				default:
					trianglesClipped++
					// outside or not fully covered -> remove.
					if i < len(indices)-1 {
						copy(indices[i:], indices[i+1:])
					}
					indices[len(indices)-1] = 0
					indices = indices[:len(indices)-1]
				}
			}
			ot.index[pos] = int32(len(indices))
		}
	}
	/*
			fmt.Fprintln(fake,
				`]
		}`)
	*/
	log.Printf("info: node tests: %d\n", nodeTests)
	log.Printf("info: nodes clipped: %d\n", nodesClipped)
	log.Printf("info: nodes all inside: %d\n", nodesAllInside)
	log.Printf("info: triangle clipped: %d\n", trianglesClipped)
	log.Printf("info: tested triangles: %d\n", len(triChecks))

	ot.eliminate(triChecks)
}

func (ot *Tree) Value(x, y float64) (float64, bool) {

	// out of bounding box
	if x < ot.Min.X || ot.Max.X < x ||
		y < ot.Min.Y || ot.Max.Y < y {
		return 0, false
	}

	all := Box2D{ot.Min.X, ot.Min.Y, ot.Max.X, ot.Max.Y}

	stack := []boxFrame{{1, all}}

	for len(stack) > 0 {
		top := stack[len(stack)-1]
		stack = stack[:len(stack)-1]

		if top.pos > 0 { // node
			if index := ot.index[top.pos:]; len(index) > 7 {
				for i := 0; i < 4; i++ {
					a := index[i]
					b := index[i+4]
					if a == 0 && b == 0 {
						continue
					}
					dx := top.X2 - top.X1
					dy := top.Y2 - top.Y1
					nbox := Box2D{
						dx*scale[i][0] + top.X1,
						dy*scale[i][1] + top.Y1,
						dx*scale[i][2] + top.X1,
						dy*scale[i][3] + top.Y1,
					}
					if nbox.Contains(x, y) {
						if a != 0 {
							stack = append(stack, boxFrame{a, nbox})
						}
						if b != 0 {
							stack = append(stack, boxFrame{b, nbox})
						}
						break
					}
				}
			}
		} else { // leaf
			pos := -top.pos - 1
			n := ot.index[pos]
			indices := ot.index[pos+1 : pos+1+n]

			for _, idx := range indices {
				tri := ot.triangles[idx]
				t := Triangle{
					ot.vertices[tri[0]],
					ot.vertices[tri[1]],
					ot.vertices[tri[2]],
				}
				if t.Contains(x, y) {
					return t.Plane3D().Z(x, y), true
				}
			}
		}
	}

	return 0, false
}

// Vertical does a vertical cross cut from (x1, y1) to (x2, y2).
func (ot *Tree) Vertical(x1, y1, x2, y2 float64, fn func(*Triangle)) {

	box := Box2D{
		X1: math.Min(x1, x2),
		Y1: math.Min(y1, y2),
		X2: math.Max(x1, x2),
		Y2: math.Max(y1, y2),
	}

	// out of bounding box
	if box.X2 < ot.Min.X || ot.Max.X < box.X1 ||
		box.Y2 < ot.Min.Y || ot.Max.Y < box.Y1 {
		return
	}

	line := NewPlane2D(x1, y1, x2, y2)

	dupes := map[int32]struct{}{}

	all := Box2D{ot.Min.X, ot.Min.Y, ot.Max.X, ot.Max.Y}
	//log.Printf("area: %f\n", (box.x2-box.x1)*(box.y2-box.y1))
	//log.Printf("all: %f\n", (all.x2-all.x1)*(all.y2-all.y1))

	stack := []boxFrame{{1, all}}

	for len(stack) > 0 {
		top := stack[len(stack)-1]
		stack = stack[:len(stack)-1]

		if top.pos > 0 { // node
			if index := ot.index[top.pos:]; len(index) > 7 {
				for i := 0; i < 4; i++ {
					a := index[i]
					b := index[i+4]
					if a == 0 && b == 0 {
						continue
					}
					dx := top.X2 - top.X1
					dy := top.Y2 - top.Y1
					nbox := Box2D{
						dx*scale[i][0] + top.X1,
						dy*scale[i][1] + top.Y1,
						dx*scale[i][2] + top.X1,
						dy*scale[i][3] + top.Y1,
					}
					if nbox.Intersects(box) && nbox.IntersectsPlane(line) {
						if a != 0 {
							stack = append(stack, boxFrame{a, nbox})
						}
						if b != 0 {
							stack = append(stack, boxFrame{b, nbox})
						}
					}
				}
			}
		} else { // leaf
			pos := -top.pos - 1
			n := ot.index[pos]
			indices := ot.index[pos+1 : pos+1+n]

			for _, idx := range indices {
				if _, found := dupes[idx]; found {
					continue
				}
				tri := ot.triangles[idx]
				t := Triangle{
					ot.vertices[tri[0]],
					ot.vertices[tri[1]],
					ot.vertices[tri[2]],
				}

				v0 := line.Eval(t[0].X, t[0].Y)
				v1 := line.Eval(t[1].X, t[1].Y)
				v2 := line.Eval(t[2].X, t[2].Y)

				if onPlane(v0) || onPlane(v1) || onPlane(v2) ||
					sides(sides(sides(0, v0), v1), v2) == 3 {
					fn(&t)
				}
				dupes[idx] = struct{}{}
			}
		}
	}
}

// Horizontal does a horizontal cross cut.
func (ot *Tree) Horizontal(h float64, fn func(*Triangle)) {

	if h < ot.Min.Z || ot.Max.Z < h {
		return
	}

	type frame struct {
		pos int32
		min float64
		max float64
	}

	dupes := map[int32]struct{}{}

	stack := []frame{{1, ot.Min.Z, ot.Max.Z}}

	for len(stack) > 0 {
		top := stack[len(stack)-1]
		stack = stack[:len(stack)-1]

		pos := top.pos
		if pos == 0 {
			continue
		}
		min, max := top.min, top.max

		if pos > 0 { // node
			if mid := (max-min)*0.5 + min; h >= mid {
				pos += 4 // nodes with z-bit set
				min = mid
			} else {
				max = mid
			}
			if pos < int32(len(ot.index)) {
				if index := ot.index[pos:]; len(index) > 3 {
					stack = append(stack,
						frame{index[0], min, max},
						frame{index[1], min, max},
						frame{index[2], min, max},
						frame{index[3], min, max})
				}
			}
		} else { // leaf
			pos = -pos - 1
			n := ot.index[pos]
			//log.Printf("%d %d %d\n", pos, n, len(ot.index))
			indices := ot.index[pos+1 : pos+1+n]

			for _, idx := range indices {
				if _, found := dupes[idx]; found {
					continue
				}
				tri := ot.triangles[idx]
				t := Triangle{
					ot.vertices[tri[0]],
					ot.vertices[tri[1]],
					ot.vertices[tri[2]],
				}

				if !(math.Min(t[0].Z, math.Min(t[1].Z, t[2].Z)) > h ||
					math.Max(t[0].Z, math.Max(t[1].Z, t[2].Z)) < h) {
					dupes[idx] = struct{}{}
					fn(&t)
				}
			}
		}
	}
}