local bezier = {}
 
local function length(n, func, ...)
	local sum, ranges, sums = 0, {}, {}

	for i = 0, n-1 do
		local p1, p2 = func(i/n, ...), func((i+1)/n, ...)
		local dist = vec2.dist(p2, p1)
		ranges[sum] = {dist, p1, p2}
		table.insert(sums, sum)
		sum = sum + dist
	end
	return sum, ranges, sums
end
 
local function quadBezier(t, p0, p1, p2)
	local v1 = p0 * (1 - t)^2 
	local v2 = p1 * 2 * (1 - t) * t
	local v3 = p2 * t^2
	return v1 + v2 + v3
end

local function cubicBezier(t, p0, p1, p2, p3)
	local l1 = lerp(p0, p1, t)
	local l2 = lerp(p1, p2, t)
	local l3 = lerp(p2, p3, t)
	local a = lerp(l1, l2, t)
	local b = lerp(l2, l3, t)
	local cubic = lerp(a, b, t)
	return cubic
end

function bezier.new(func, n, ...)
	local pts = {...}
	local func = func or (#pts == 8 and cubicBezier or quadBezier)
	local self = setmetatable({}, {__index = bezier})
	local sum, ranges, sums = length(n, func, ...)
	self.func = func 
	self.n = n
	self.points = pts
	self.length = sum
	self.ranges = ranges
	self.sums = sums
	return self
end
 
function bezier:setPoints(...)
	-- only update the length when the control points are changed
	local sum, ranges, sums = length(self.n, self.func, ...)
	self.points = {...}
	self.length = sum
	self.ranges = ranges
	self.sums = sums
end

function bezier:setControlPoint(i,x,y)
	-- print(inspect(self.points))
	-- print_info(3)
	self.points[i].x = x
	self.points[i].y = y
	local sum, ranges, sums = length(self.n, self.func, unpack(self.points))
	self.length = sum
	self.ranges = ranges
	self.sums = sums
end
 
function bezier:calc(t)
	-- if you don't need t to be a percentage of distance
	return self.func(t, unpack(self.points))
end

function bezier:evaluate(t)
	return self:calc(t):unpack()
end
 
function bezier:calcFixed(t)
	local T, near = t * self.length, 0
	for _, n in next, self.sums do
		if (T - n) < 0 then break end
		near = n
	end
	local set = self.ranges[near]
	local percent = (T - near)/set[1]
	return set[2], set[3], percent
end

-- https://www.particleincell.com/2013/cubic-line-intersection/
-- function cubicRoots(P)
--     var a=P[0];
--     var b=P[1];
--     var c=P[2];
--     var d=P[3];
 
--     var A=b/a;
--     var B=c/a;
--     var C=d/a;
 
--     var Q, R, D, S, T, Im;
 
--     var Q = (3*B - Math.pow(A, 2))/9;
--     var R = (9*A*B - 27*C - 2*Math.pow(A, 3))/54;
--     var D = Math.pow(Q, 3) + Math.pow(R, 2);    // polynomial discriminant
 
--     var t=Array();
 
--     if (D >= 0)                                 // complex or duplicate roots
--     {
--         var S = sgn(R + Math.sqrt(D))*Math.pow(Math.abs(R + Math.sqrt(D)),(1/3));
--         var T = sgn(R - Math.sqrt(D))*Math.pow(Math.abs(R - Math.sqrt(D)),(1/3));
 
--         t[0] = -A/3 + (S + T);                    // real root
--         t[1] = -A/3 - (S + T)/2;                  // real part of complex root
--         t[2] = -A/3 - (S + T)/2;                  // real part of complex root
--         Im = Math.abs(Math.sqrt(3)*(S - T)/2);    // complex part of root pair   
 
--         /*discard complex roots*/
--         if (Im!=0)
--         {
--             t[1]=-1;
--             t[2]=-1;
--         }
 
--     }
--     else                                          // distinct real roots
--     {
--         var th = Math.acos(R/Math.sqrt(-Math.pow(Q, 3)));
 
--         t[0] = 2*Math.sqrt(-Q)*Math.cos(th/3) - A/3;
--         t[1] = 2*Math.sqrt(-Q)*Math.cos((th + 2*Math.PI)/3) - A/3;
--         t[2] = 2*Math.sqrt(-Q)*Math.cos((th + 4*Math.PI)/3) - A/3;
--         Im = 0.0;
--     }
 
--     /*discard out of spec roots*/
--     for (var i=0;i<3;i++) 
--         if (t[i]<0 || t[i]>1.0) t[i]=-1;
 
--     /*sort but place -1 at the end*/
--     t=sortSpecial(t);
 
--     console.log(t[0]+" "+t[1]+" "+t[2]);
--     return t;
-- }

-- --computes intersection between a cubic spline and a line segment
-- function computeIntersections(px,py,lx,ly)
-- {
--     var X=Array();
 
--     var A=ly[1]-ly[0];      //A=y2-y1
--     var B=lx[0]-lx[1];      //B=x1-x2
--     var C=lx[0]*(ly[0]-ly[1]) + 
--           ly[0]*(lx[1]-lx[0]);  //C=x1*(y1-y2)+y1*(x2-x1)
 
--     var bx = bezierCoeffs(px[0],px[1],px[2],px[3]);
--     var by = bezierCoeffs(py[0],py[1],py[2],py[3]);
 
--     var P = Array();
--     P[0] = A*bx[0]+B*by[0];   /*t^3*/
--     P[1] = A*bx[1]+B*by[1];   /*t^2*/
--     P[2] = A*bx[2]+B*by[2];   /*t*/
--     P[3] = A*bx[3]+B*by[3] + C; /*1*/
 
--     var r=cubicRoots(P);
 
--     /*verify the roots are in bounds of the linear segment*/  
--     for (var i=0;i<3;i++)
--     {
--         t=r[i];
 
--         X[0]=bx[0]*t*t*t+bx[1]*t*t+bx[2]*t+bx[3];
--         X[1]=by[0]*t*t*t+by[1]*t*t+by[2]*t+by[3];            
 
--         /*above is intersection point assuming infinitely long line segment,
--           make sure we are also in bounds of the line*/
--         var s;
--         if ((lx[1]-lx[0])!=0)           /*if not vertical line*/
--             s=(X[0]-lx[0])/(lx[1]-lx[0]);
--         else
--             s=(X[1]-ly[0])/(ly[1]-ly[0]);
 
--         /*in bounds?*/
--         if (t<0 || t>1.0 || s<0 || s>1.0)
--         {
--             X[0]=-100;  /*move off screen*/
--             X[1]=-100;
--         }
 
--         /*move intersection point*/
--         I[i].setAttributeNS(null,"cx",X[0]);
--         I[i].setAttributeNS(null,"cy",X[1]);
--     }    
-- }

 -- float[] BezierCoeffs(float P0, float P1, float P2, float P3) {
 --        float[] Z = new float[4];
 --        Z[0] = -P0 + 3 * P1 + -3 * P2 + P3;
 --        Z[1] = 3 * P0 - 6 * P1 + 3 * P2;
 --        Z[2] = -3 * P0 + 3 * P1;
 --        Z[3] = P0;
 --        return Z;
 --    }
return bezier