#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>
//#include "util/box.h"
#include "util/vec3f.h"
#include "bspline.h"
//B-Splines for use in field modeling shape functions
//For a nice reference, see 
//Finite Element Methods with B-Splines, Klaus Hollig


/*
Define a univariate n-degree b spline.
Here the grid size is one, we scale below.s
n=1 gives linear hat function,
n=2 quadratic,
*/
float bspline_base(const float xx, const int n)
{
	const float x = abs(xx); //symmetric about origin
	if (n==1) {
		if (x<1) return 1-x;
		else return 0;
	} else if (n==2) {
		if (x<.5) return 1+x-(x+.5)*(x+.5);
		else if (x<1.5) return 1-x+.5*(x-.5)*(x-.5);
		else return 0.;
	} else {
		assert(0); //bspline degree not implemented
		return 0;
	}
}

/*
Define corresponding derivative
*/
float d_bspline_base(const float x, const int n)
{
	if (n==1){
		if (0<=x<=1) return -1;
		else if (-1<=x<0) return 1;
		else return 0;
	} else if (n==2){
		if (-.5<=x<=.5) return -2*x;
		else if (.5<x<=1.5) return -1.5+x;
		else if (-1.5<=x<-.5) return 1.5+x;
		else return 0.;
	} else{
		assert(0);
		return 0;
	}
}

float bspline(const float x, const int n, const float x0, const float h){
	return bspline_base( (x-x0)/h, n);
}
float d_bspline(const float x, const int n, const float x0, const float h){
	return (1./h)*d_bspline_base( (x-x0)/h, n);
}
float bspline_3(const Vec3f X, const int n, const Vec3f X0, const float h){
	return bspline(X.x,n,X0.x,h)*bspline(X.y,n,X0.y,h)*bspline(X.z,n,X0.z,h);
}
float dx_bspline_3(const Vec3f X, const int n, const Vec3f X0, const float h){
	return d_bspline(X.x,n,X0.x,h)*bspline(X.y,n,X0.y,h)*bspline(X.z,n,X0.z,h);
}
float dy_bspline_3(const Vec3f X, const int n, const Vec3f X0, const float h){
	return bspline(X.x,n,X0.x,h)*d_bspline(X.y,n,X0.y,h)*bspline(X.z,n,X0.z,h);
}
float dz_bspline_3(const Vec3f X, const int n, const Vec3f X0, const float h){
	return bspline(X.x,n,X0.x,h)*bspline(X.y,n,X0.y,h)*d_bspline(X.z,n,X0.z,h);
}
//convenience
Vec3f grad_bspline_3(const Vec3f X, const int n, const Vec3f X0, const float h){
	return (Vec3f){
		d_bspline(X.x,n,X0.x,h)*bspline(X.y,n,X0.y,h)*bspline(X.z,n,X0.z,h),
		bspline(X.x,n,X0.x,h)*d_bspline(X.y,n,X0.y,h)*bspline(X.z,n,X0.z,h),
		bspline(X.x,n,X0.x,h)*bspline(X.y,n,X0.y,h)*d_bspline(X.z,n,X0.z,h)
	};
}

Box bspline_support(const int n, const Vec3f X0, const float h)
{
	//const Vec3f rad = (Vec3f){.5*(n-1)*h,.5*(n-1)*h,.5*(n-1)*h}; //radius of support
	const Vec3f rad = (Vec3f){.5*(n+1)*h,.5*(n+1)*h,.5*(n+1)*h}; //radius of support
	return (Box){vdiff(X0,rad), vsum(X0,rad)};
}




/*
//Here's a hook for a bspline centered and stretched on an Interval
float b_spline_k(const float x, const int n, const Interval I)
{
	const float middle = .5*(I.upper+I.lower);
	const float radius = .5*(I.upper-I.lower);
	return b_spline( (x-middle)/radius, n );
}
float d_b_spline_k(const float x, const int n, const Interval I)
{
	const float middle = .5*(I.upper+I.lower);
	const float radius = .5*(I.upper-I.lower);
	return (1./radius)*d_b_spline( (x-middle)/radius, n );
}


//Now for multivariate, we just do products of univariates
float b_spline_3k(const Vec3f X, const int n, const Interval Ix, const Interval Iy, const Interval Iz)
{
	return b_spline_k(X.x,n,Ix)*b_spline_k(X.y,n,Iy)*b_spline_k(X.z,n,Iz);
}

float dx_b_spline_3k(const Vec3f X, const int n, const Interval Ix, const Interval Iy, const Interval Iz)
{return d_b_spline_k(X.x,n,Ix)*b_spline_k(X.y,n,Iy)*b_spline_k(X.z,n,Iz);}

float dy_b_spline_3k(const Vec3f X, const int n, const Interval Ix, const Interval Iy, const Interval Iz)
{return b_spline_k(X.x,n,Ix)*d_b_spline_k(X.y,n,Iy)*b_spline_k(X.z,n,Iz);}

float dz_b_spline_3k(const Vec3f X, const int n, const Interval Ix, const Interval Iy, const Interval Iz)
{return b_spline_k(X.x,n,Ix)*b_spline_k(X.y,n,Iy)*d_b_spline_k(X.z,n,Iz);}
*/