// Copyright ©2015 The gonum Authors. All rights reserved.
// Use of this source code is governed by a BSD-style
// license that can be found in the LICENSE file.

package conrec

import (
	"math"
	"sort"

	"github.com/gonum/matrix/mat64"
)

type point struct {
	X, Y float64
}

type line struct {
	p1, p2 point
}

func sect(h, p [5]float64, v1, v2 int) float64 {
	return (h[v2]*p[v1] - h[v1]*p[v2]) / (h[v2] - h[v1])
}

// ConrecLine performs an operation with a line at a given height derived
// from data over the 2D interval (i, j) to (i+1, j+1).
type ConrecLine func(i, j int, l line, height float64)

// Conrec is a Go translation of C version of CONREC by Paul Bourke:
//  http://paulbourke.net/papers/conrec/conrec.c
//
// Conrec takes d, an m×n matrix, an m-long row coordinates slice, an
// n-long column coordinates slice, a slice of contour heights and
// a ConrecLine function.
// Conrec panics if len(row) != m or len(col) != n.
//
// For full details of the algorithm, see the paper at
// http://paulbourke.net/papers/conrec/
func Conrec(d mat64.Matrix, row, col, heights []float64, fn ConrecLine) {
	r, c := d.Dims()
	if len(row) != r {
		panic("conrec: row coordinates mismatch matrix row dimension")
	}
	if len(col) != c {
		panic("conrec: col coordinates mismatch matrix col dimension")
	}

	sort.Float64s(heights)

	var (
		p1, p2 point

		h      [5]float64
		sh     [5]int
		xh, yh [5]float64

		im = [4]int{0, 0, 1, 1}
		jm = [4]int{0, 1, 1, 0}

		cases = [3][3][3]int{
			{{0, 0, 8}, {0, 2, 5}, {7, 6, 9}},
			{{0, 3, 4}, {1, 3, 1}, {4, 3, 0}},
			{{9, 6, 7}, {5, 2, 0}, {8, 0, 0}},
		}
	)

	for i, rv := range row[:len(row)-1] {
		for j, cv := range col[:len(col)-1] {
			dmin := math.Min(
				math.Min(d.At(i, j), d.At(i, j+1)),
				math.Min(d.At(i+1, j), d.At(i+1, j+1)),
			)

			dmax := math.Max(
				math.Max(d.At(i, j), d.At(i, j+1)),
				math.Max(d.At(i+1, j), d.At(i+1, j+1)),
			)

			if dmax < heights[0] || heights[len(heights)-1] < dmin {
				continue
			}

			for k := 0; k < len(heights); k++ {
				if heights[k] < dmin || dmax < heights[k] {
					continue
				}
				for m := 4; m >= 0; m-- {
					if m > 0 {
						h[m] = d.At(i+im[m-1], j+jm[m-1]) - heights[k]
						xh[m] = col[j+jm[m-1]]
						yh[m] = row[i+im[m-1]]
					} else {
						h[0] = 0.25 * (h[1] + h[2] + h[3] + h[4])
						xh[0] = 0.50 * (cv + col[j+1])
						yh[0] = 0.50 * (rv + row[i+1])
					}
					if h[m] > 0 {
						sh[m] = 1
					} else if h[m] < 0 {
						sh[m] = -1
					} else {
						sh[m] = 0
					}
				}

				/*
				   Note: at this stage the relative heights of the corners and the
				   centre are in the h array, and the corresponding coordinates are
				   in the xh and yh arrays. The centre of the box is indexed by 0
				   and the 4 corners by 1 to 4 as shown below.
				   Each triangle is then indexed by the parameter m, and the 3
				   vertices of each triangle are indexed by parameters m1,m2,and m3.
				   It is assumed that the centre of the box is always vertex 2
				   though this isimportant only when all 3 vertices lie exactly on
				   the same contour level, in which case only the side of the box
				   is drawn.

				      vertex 4 +-------------------+ vertex 3
				               | \               / |
				               |   \    m-3    /   |
				               |     \       /     |
				               |       \   /       |
				               |  m=2    X   m=2   |       the centre is vertex 0
				               |       /   \       |
				               |     /       \     |
				               |   /    m=1    \   |
				               | /               \ |
				      vertex 1 +-------------------+ vertex 2
				*/

				/* Scan each triangle in the box */
				for m := 1; m <= 4; m++ {
					m1 := m
					const m2 = 0
					var m3 int
					if m != 4 {
						m3 = m + 1
					} else {
						m3 = 1
					}
					switch cases[sh[m1]+1][sh[m2]+1][sh[m3]+1] {
					case 0:
						continue
					case 1: // Line between vertices 1 and 2
						p1 = point{X: xh[m1], Y: yh[m1]}
						p2 = point{X: xh[m2], Y: yh[m2]}

					case 2: // Line between vertices 2 and 3
						p1 = point{X: xh[m2], Y: yh[m2]}
						p2 = point{X: xh[m3], Y: yh[m3]}

					case 3: // Line between vertices 3 and 1
						p1 = point{X: xh[m3], Y: yh[m3]}
						p2 = point{X: xh[m1], Y: yh[m1]}

					case 4: // Line between vertex 1 and side 2-3
						p1 = point{X: xh[m1], Y: yh[m1]}
						p2 = point{X: sect(h, xh, m2, m3), Y: sect(h, yh, m2, m3)}

					case 5: // Line between vertex 2 and side 3-1
						p1 = point{X: xh[m2], Y: yh[m2]}
						p2 = point{X: sect(h, xh, m3, m1), Y: sect(h, yh, m3, m1)}

					case 6: // Line between vertex 3 and side 1-2
						p1 = point{X: xh[m3], Y: yh[m3]}
						p2 = point{X: sect(h, xh, m1, m2), Y: sect(h, yh, m1, m2)}

					case 7: // Line between sides 1-2 and 2-3
						p1 = point{X: sect(h, xh, m1, m2), Y: sect(h, yh, m1, m2)}
						p2 = point{X: sect(h, xh, m2, m3), Y: sect(h, yh, m2, m3)}

					case 8: // Line between sides 2-3 and 3-1
						p1 = point{X: sect(h, xh, m2, m3), Y: sect(h, yh, m2, m3)}
						p2 = point{X: sect(h, xh, m3, m1), Y: sect(h, yh, m3, m1)}

					case 9: // Line between sides 3-1 and 1-2
						p1 = point{X: sect(h, xh, m3, m1), Y: sect(h, yh, m3, m1)}
						p2 = point{X: sect(h, xh, m1, m2), Y: sect(h, yh, m1, m2)}

					default:
						panic("cannot reach")
					}

					fn(i, j, line{p1: p1, p2: p2}, heights[k])
				}
			}
		}
	}
}