func (c *contour) extend(l line, ends endMap) (ok bool) {
	switch c.front() {
	case l.p1:
		c.backward = append(c.backward, l.p2)
		delete(ends, l.p1)
		ends[l.p2] = c
		return true
	case l.p2:
		c.backward = append(c.backward, l.p1)
		delete(ends, l.p2)
		ends[l.p1] = c
		return true
	}
	switch c.back() {
	case l.p1:
		c.forward = append(c.forward, l.p2)
		delete(ends, l.p1)
		ends[l.p2] = c
		return true
	case l.p2:
		c.forward = append(c.forward, l.p1)
		delete(ends, l.p2)
		ends[l.p1] = c
		return true
	}
	return false
}
extend adds the line l to the contour, updating the ends map. It returns
a boolean indicating whether the extension was successful.
package plotter
import (
	"image/color"
	"math"
	"sort"
	"gonum.org/v1/plot"
	"gonum.org/v1/plot/palette"
	"gonum.org/v1/plot/vg"
	"gonum.org/v1/plot/vg/draw"
)
type Contour struct {
	GridXYZ GridXYZ
	Levels []float64
	LineStyles []draw.LineStyle
	Palette palette.Palette
	Underflow color.Color
	Overflow  color.Color
	Min, Max float64
}
func NewContour(g GridXYZ, levels []float64, p palette.Palette) *Contour {
	var min, max float64
	type minMaxer interface {
		Min() float64
		Max() float64
	}
	switch g := g.(type) {
	case minMaxer:
		min, max = g.Min(), g.Max()
	default:
		min, max = math.Inf(1), math.Inf(-1)
		c, r := g.Dims()
		for i := 0; i < c; i++ {
			for j := 0; j < r; j++ {
				v := g.Z(i, j)
				if math.IsNaN(v) {
					continue
				}
				min = math.Min(min, v)
				max = math.Max(max, v)
			}
		}
	}
	if len(levels) == 0 {
		levels = quantilesR7(g, defaultQuantiles)
	}
	return &Contour{
		GridXYZ:    g,
		Levels:     levels,
		LineStyles: []draw.LineStyle{DefaultLineStyle},
		Palette:    p,
		Min:        min,
		Max:        max,
	}
}
var defaultQuantiles = []float64{0.01, 0.05, 0.25, 0.5, 0.75, 0.95, 0.99}
func quantilesR7(g GridXYZ, p []float64) []float64 {
	c, r := g.Dims()
	data := make([]float64, 0, c*r)
	for i := 0; i < c; i++ {
		for j := 0; j < r; j++ {
			if v := g.Z(i, j); !math.IsNaN(v) {
				data = append(data, v)
			}
		}
	}
	sort.Float64s(data)
	v := make([]float64, len(p))
	for j, q := range p {
		if q == 1 {
			v[j] = data[len(data)-1]
		}
		h := float64(len(data)-1) * q
		i := int(h)
		v[j] = data[i] + (h-math.Floor(h))*(data[i+1]-data[i])
	}
	return v
}
const naive = false
func (h *Contour) Plot(c draw.Canvas, plt *plot.Plot) {
	if h.Min > h.Max {
		panic("contour: invalid Z range: min greater than max")
	}
	if naive {
		h.naivePlot(c, plt)
		return
	}
	var pal []color.Color
	if h.Palette != nil {
		pal = h.Palette.Colors()
	}
	trX, trY := plt.Transforms(&c)
	cp := contourPaths(h.GridXYZ, h.Levels, trX, trY)
	ps := float64(len(pal)-1) / (h.Levels[len(h.Levels)-1] - h.Levels[0])
	if len(h.Levels) == 1 {
		ps = 0
	}
	for i, z := range h.Levels {
		if math.IsNaN(z) {
			continue
		}
		for _, pa := range cp[z] {
			if isLoop(pa) {
				pa.Close()
			}
			style := h.LineStyles[i%len(h.LineStyles)]
			var col color.Color
			switch {
			case z < h.Min:
				col = h.Underflow
			case z > h.Max:
				col = h.Overflow
			case len(pal) == 0:
				col = style.Color
			default:
				col = pal[int((z-h.Levels[0])*ps+0.5)]  
			}
			if col != nil && style.Width != 0 {
				c.SetLineStyle(style)
				c.SetColor(col)
				c.Stroke(pa)
			}
		}
	}
}
func (h *Contour) naivePlot(c draw.Canvas, plt *plot.Plot) {
	var pal []color.Color
	if h.Palette != nil {
		pal = h.Palette.Colors()
	}
	trX, trY := plt.Transforms(&c)
	sort.Float64s(h.Levels)
	ps := float64(len(pal)-1) / (h.Levels[len(h.Levels)-1] - h.Levels[0])
	if len(h.Levels) == 1 {
		ps = 0
	}
	levelMap := make(map[float64]int)
	for i, z := range h.Levels {
		levelMap[z] = i
	}
	var pa vg.Path
	conrec(h.GridXYZ, h.Levels, func(_, _ int, l line, z float64) {
		if math.IsNaN(z) {
			return
		}
		pa = pa[:0]
		x1, y1 := trX(l.p1.X), trY(l.p1.Y)
		x2, y2 := trX(l.p2.X), trY(l.p2.Y)
		pt1 := vg.Point{X: x1, Y: y1}
		pt2 := vg.Point{X: x2, Y: y2}
		if !c.Contains(pt1) || !c.Contains(pt2) {
			return
		}
		pa.Move(pt1)
		pa.Line(pt2)
		pa.Close()
		style := h.LineStyles[levelMap[z]%len(h.LineStyles)]
		var col color.Color
		switch {
		case z < h.Min:
			col = h.Underflow
		case z > h.Max:
			col = h.Overflow
		case len(pal) == 0:
			col = style.Color
		default:
			col = pal[int((z-h.Levels[0])*ps+0.5)]  
		}
		if col != nil && style.Width != 0 {
			c.SetLineStyle(style)
			c.SetColor(col)
			c.Stroke(pa)
		}
	})
}
func (h *Contour) DataRange() (xmin, xmax, ymin, ymax float64) {
	c, r := h.GridXYZ.Dims()
	return h.GridXYZ.X(0), h.GridXYZ.X(c - 1), h.GridXYZ.Y(0), h.GridXYZ.Y(r - 1)
}
func (h *Contour) GlyphBoxes(plt *plot.Plot) []plot.GlyphBox {
	c, r := h.GridXYZ.Dims()
	b := make([]plot.GlyphBox, 0, r*c)
	for i := 0; i < c; i++ {
		for j := 0; j < r; j++ {
			b = append(b, plot.GlyphBox{
				X: plt.X.Norm(h.GridXYZ.X(i)),
				Y: plt.Y.Norm(h.GridXYZ.Y(j)),
				Rectangle: vg.Rectangle{
					Min: vg.Point{X: -2.5, Y: -2.5},
					Max: vg.Point{X: +2.5, Y: +2.5},
				},
			})
		}
	}
	return b
}
func isLoop(p vg.Path) bool {
	s := p[0]
	e := p[len(p)-1]
	return s.Pos == e.Pos
}
func contourPaths(m GridXYZ, levels []float64, trX, trY func(float64) vg.Length) map[float64][]vg.Path {
	sort.Float64s(levels)
	ends := make(map[float64]endMap)
	conts := make(contourSet)
	conrec(m, levels, func(_, _ int, l line, z float64) {
		paths(l, z, ends, conts)
	})
	ends = nil
	for c := range conts {
		c.exciseLoops(conts, true)
	}
	paths := make(map[float64][]vg.Path)
	for c := range conts {
		paths[c.z] = append(paths[c.z], c.path(trX, trY))
	}
	return paths
}
type contourSet map[*contour]struct{}
type endMap map[point]*contour
func paths(l line, z float64, ends map[float64]endMap, conts contourSet) {
	zEnds, ok := ends[z]
	if !ok {
		zEnds = make(endMap)
		ends[z] = zEnds
		c := newContour(l, z)
		zEnds[l.p1] = c
		zEnds[l.p2] = c
		conts[c] = struct{}{}
		return
	}
	c1, ok1 := zEnds[l.p1]
	c2, ok2 := zEnds[l.p2]
	if !ok1 && !ok2 {
		c := newContour(l, z)
		zEnds[l.p1] = c
		zEnds[l.p2] = c
		conts[c] = struct{}{}
		return
	}
	if ok1 {
		if !c1.extend(l, zEnds) {
			panic("internal link")
		}
	} else if ok2 {
		if !c2.extend(l, zEnds) {
			panic("internal link")
		}
	}
	if c1 == c2 {
		return
	}
	if ok1 && ok2 {
		if !c1.connect(c2, zEnds) {
			panic("internal link")
		}
		delete(conts, c2)
	}
}
type path []point
type contour struct {
	z float64
	backward path
	forward  path
}
func newContour(l line, z float64) *contour {
	return &contour{z: z, forward: path{l.p1}, backward: path{l.p2}}
}
func (c *contour) path(trX, trY func(float64) vg.Length) vg.Path {
	var pa vg.Path
	p := c.front()
	pa.Move(vg.Point{X: trX(p.X), Y: trY(p.Y)})
	for i := len(c.backward) - 2; i >= 0; i-- {
		p = c.backward[i]
		pa.Line(vg.Point{X: trX(p.X), Y: trY(p.Y)})
	}
	for _, p := range c.forward {
		pa.Line(vg.Point{X: trX(p.X), Y: trY(p.Y)})
	}
	return pa
}
func (c *contour) front() point { return c.backward[len(c.backward)-1] }
func (c *contour) back() point { return c.forward[len(c.forward)-1] }