func (ss *CubicSpline) Polygonize() SDF2 {
	p := NewPolygon()
	p.Add(ss.xmin, 0)
	p.Add(ss.xmax, 0)
	dx := (ss.xmax - ss.xmin) / float64(N_SAMPLES-1)
	x := ss.xmax
	for i := 0; i < N_SAMPLES; i++ {
		p.Add(x, ss.Function(x))
		x -= dx
	}
	p.Render("spline.dxf")
	return Polygon2D(p.Vertices())
}
Return a 2D polygon approximating the cubic spline.
package sdf
import (
	"fmt"
	"math"
)
func TriDiagonal(m []V3, d []float64) []float64 {
	n := len(m)
	if len(d) != n {
		panic("bad sizes rows(m) != rows(d)")
	}
	if m[0].X != 0 || m[n-1].Z != 0 {
		panic("bad values for tridiagonal matrix")
	}
	if m[0].Y == 0 {
		panic("m[0].Y == 0")
	}
	cp := make([]float64, n)  
	x := make([]float64, n)   
	cp[0] = m[0].Z / m[0].Y
	x[0] = d[0] / m[0].Y
	for i := 1; i < n; i++ {
		denom := m[i].Y - m[i].X*cp[i-1]
		if denom == 0 {
			panic("denom == 0")
		}
		cp[i] = m[i].Z / denom
		x[i] = (d[i] - m[i].X*x[i-1]) / denom
	}
	for i := n - 2; i >= 0; i-- {
		x[i] -= cp[i] * x[i+1]
	}
	return x
}
type CS struct {
	p0, p1     V2       
	idx        int      
	k          float64  
	a, b, c, d float64  
}
type CubicSpline struct {
	xmin, xmax float64
	spline     []CS
}
func NewCubicSpline(data []V2) CubicSpline {
	n := len(data)
	m := make([]V3, n)
	d := make([]float64, n)
	for i := 1; i < n-1; i++ {
		m[i] = V3{1, 4, 1}
		d[i] = 3 * (data[i+1].Y - data[i-1].Y)
	}
	m[0] = V3{0, 2, 1}
	d[0] = 3 * (data[1].Y - data[0].Y)
	m[n-1] = V3{1, 2, 0}
	d[n-1] = 3 * (data[n-1].Y - data[n-2].Y)
	x := TriDiagonal(m, d)
	spline := make([]CS, n-1)
	for i := 0; i < n-1; i++ {
		x0 := data[i].X
		x1 := data[i+1].X
		y0 := data[i].Y
		y1 := data[i+1].Y
		D0 := x[i]
		D1 := x[i+1]
		spline[i].idx = i
		spline[i].p0 = data[i]
		spline[i].p1 = data[i+1]
		spline[i].k = 1.0 / (x1 - x0)
		spline[i].a = y0
		spline[i].b = D0
		spline[i].c = 3*(y1-y0) - 2*D0 - D1
		spline[i].d = 2*(y0-y1) + D0 + D1
	}
	return CubicSpline{data[0].X, data[n-1].X, spline}
}
func (s *CS) XtoT(x float64) float64 {
	return s.k * (x - s.p0.X)
}
func (s *CS) Function(t float64) float64 {
	return s.a + t*(s.b+t*(s.c+s.d*t))
}
func (s *CS) FirstDerivative(t float64) float64 {
	return s.b + t*(2*s.c+3*s.d*t)
}
func (s *CS) SecondDerivative(t float64) float64 {
	return 2*s.c + 6*s.d*t
}
func (s *CS) Max() float64 {
	ymax := Max(s.p0.Y, s.p1.Y)
	t, n := quadratic(3*s.d, 2*s.c, s.b)
	if n == ONE_SOLN {
		ymax = Max(ymax, s.Function(t[0]))
	} else if n == TWO_SOLN {
		ymax = Max(ymax, s.Function(t[0]))
		ymax = Max(ymax, s.Function(t[1]))
	}
	return ymax
}
func (ss CubicSpline) Find(x float64) *CS {
	n := len(ss.spline)
	if n == 0 {
		panic("no splines")
	}
	lo := 0
	hi := n
	for hi-lo > 1 {
		mid := (lo + hi) >> 1
		if ss.spline[mid].p0.X < x {
			lo = mid
		} else {
			hi = mid
		}
	}
	return &ss.spline[lo]
}
func (ss CubicSpline) Function(x float64) float64 {
	s := ss.Find(x)
	return s.Function(s.XtoT(x))
}
func (ss *CubicSpline) Dist2(x float64, p V2) float64 {
	dx := x - p.X
	dy := ss.Function(x) - p.Y
	return dx*dx + dy*dy
}
const N_SAMPLES = 30