//アキマ(秋間 浩)の補間法２ for WIN32 API

#include"Akima_NI.h"

double akima::zet(double d)		//絶対値(zettaichi)
{
	return (d<0)?d*(-1):d;
}

double akima::pow(double base,int n)
{
	double d=1;
	for(int i=0;i<n;i++)
		d*=base;
	return d;
}

double akima::m(int j)
{
	return (f[j+1]-f[j]) / (x[j+1]-x[j]);
}

double akima::t(int j)
{
	double d=zet(m(j+1)-m(j))+zet(m(j-1)-m(j-2));
	if(d)		//d!=0
		return (zet(m(j+1)-m(j))*m(j-1)+zet(m(j-1)-m(j-2))*m(j)) / d;
	else		//d==0
		return (m(j-1)+m(j)) / 2;
}

double akima::h(int j)
{
	return x[j+1]-x[j];
}

double akima::a(int j,int k)
{
	switch(k)
	{
		case 0:		return f[j];
		case 1:		return t(j);
		case 2:		return (3*m(j)-2*t(j)-t(j+1)) / h(j);
		case 3:		return (t(j)+t(j+1)-2*m(j)) / pow(h(j),2);
	}
	return -1;		//error
}

double akima::P(int j,double xx)
{
	double d=0;
	for(int i=0;i<=3;i++)
		d+=a(j,i)*pow(xx-x[j],i);		//↓展開式
	return d;
	//a(j,0) + a(j,1)*(xx-x[j]) + a(j,2)*(xx-x[j])*(xx-x[j]) + a(j,3)*(xx-x[j])*(xx-x[j])*(xx-x[j]);
}

double akima::A(int k,double xx,int n)
{
	double d=1;
	for(int i=n-2;i<=n;i++)
		if(i!=k)
			d*=(xx-x[i]) / (x[k]-x[i]);
	return d;
}

double akima::q(double xx,int n)
{
	double d=0;
	for(int k=n-2;k<=n;k++)
		d+=A(k,xx,n)*f[k];
	return d;
}

void akima::add_point(int pt[],int py[],int i)		//使える分点は利用、使えないのなら算出
{	
	if(x[2]<x[4]){
		if(i>=1 && pt[i-1]<x[2]){
			x[1]=pt[i-1];	f[1]=py[i-1];
			if(i>=2 && pt[i-2]<pt[i-1]){
				x[0]=pt[i-2];	f[0]=py[i-2];
			}else{
				x[0]=x[1]+x[2]-x[3];	f[0]=q(x[0],3);
			}
		}else{
			x[0]=x[2]+x[2]-x[4];	f[0]=q(x[0],4);
			x[1]=x[2]+x[3]-x[4];	f[1]=q(x[1],4);
		}
		if(i<=N-4 && x[4]<pt[i+3]){
			x[5]=pt[i+3];	f[5]=py[i+3];
			if(i<=N-5 && pt[i+3]<pt[i+4]){
				x[6]=pt[i+4];	f[6]=py[i+4];
			}else{
				x[6]=x[5]+x[4]-x[3];	f[6]=q(x[6],5);
			}
		}else{
			x[5]=x[4]+x[3]-x[2];	f[5]=q(x[5],4);
			x[6]=x[4]+x[4]-x[2];	f[6]=q(x[6],4);
		}
	}else{
		if(i>=1 && pt[i-1]>x[2]){
			x[1]=pt[i-1];	f[1]=py[i-1];
			if(i>=2 && pt[i-2]>pt[i-1]){
				x[0]=pt[i-2];	f[0]=py[i-2];
			}else{
				x[0]=x[1]+x[2]-x[3];	f[0]=q(x[0],3);
			}
		}else{
			x[0]=x[2]+x[2]-x[4];	f[0]=q(x[0],4);
			x[1]=x[2]+x[3]-x[4];	f[1]=q(x[1],4);
		}
		if(i<=N-4 && x[4]>pt[i+3]){
			x[5]=pt[i+3];	f[5]=py[i+3];
			if(i<=N-5 && pt[i+3]>pt[i+4]){
				x[6]=pt[i+4];	f[6]=py[i+4];
			}else{
				x[6]=x[5]+x[4]-x[3];	f[6]=q(x[6],5);
			}
		}else{
			x[5]=x[4]+x[3]-x[2];	f[5]=q(x[5],4);
			x[6]=x[4]+x[4]-x[2];	f[6]=q(x[6],4);
		}
	}
}

int akima::count_point(int pt[],int i)		//使える分点の数を数える
{
	int n=0;
	if(pt[i]<pt[i+2]){
		if(i>=1 && pt[i-1]<pt[i]){		//側近の分点(前)
			n+=3;		//重み大
			if(i>=2 && pt[i-2]<pt[i-1])		//外側の分点(前)
				n+=2;		//重み中
		}
		if(i<=N-4 && pt[i+2]<pt[i+3]){		//側近の分点(後)
			n+=2;		//重み中
			if(i<=N-5 && pt[i+3]<pt[i+4])		//外側の分点(後)
				n++;		//重み小
		}
	}else{
		if(i>=1 && pt[i-1]>pt[i]){
			n+=3;
			if(i>=2 && pt[i-2]>pt[i-1])
				n+=2;
		}
		if(i<=N-4 && pt[i+2]>pt[i+3]){
			n+=2;
			if(i<=N-5 && pt[i+3]>pt[i+4])
				n++;
		}
	}
	return n;
}

int akima::select_axis(int i)		//x方向を基準とするか、y方向を基準とするか
{
	int xs=PX[i],xc=PX[i+1],xe=PX[i+2],ys=PF[i],yc=PF[i+1],ye=PF[i+2];
	if(xs>xe){xs=xe-xs;	xe-=xs;	xs+=xe;}		//2値交換
	if(ys>ye){ys=ye-ys;	ye-=ys;	ys+=ye;}		//xs<xe , ys<ye とする

	if(xs<xc && xc<xe && ys<yc && yc<ye)		//境界は含まない
		if(count_point(PX,i)<count_point(PF,i))		//どちらの方が多く分点を使えるか
			return 1;		//y方向を基準とする事を示す値は0以外
		else if(count_point(PX,i)>count_point(PF,i))
			return 0;		//x方向を基準とする事を示す値は0
		else		//範囲条件より分母が0になる事は有り得ない
			//return zet((PF[i+1]-PF[i])/(PX[i+1]-PX[i]))<1 ? 0:1;		//部分的には滑らかだが、全体としては不自然
			return (ye-ys)/(xe-xs)<1 ? 0:1;		//部分的にはカクカクしてしまうが、全体としては自然
			//return (xe-xs)?(ye-ys)/(xe-xs):1;		//int型という前提があれば
	else if(xs<xc && xc<xe)		//境界は含まない
		return 0;
	else if(ys<yc && yc<ye)		//境界は含まない
		return 1;

	return -1;		//error
}

int akima::set_point(int i)		//計算用の配列に分点の座標を設定
{
	if(select_axis(i)){		//基準方向判定
		f[2]=PX[i];		x[2]=PF[i];		//y方向を基準とする
		f[3]=PX[i+1];	x[3]=PF[i+1];
		f[4]=PX[i+2];	x[4]=PF[i+2];
		add_point(PF,PX,i);
		return 1;
	}else{
		x[2]=PX[i];		f[2]=PF[i];		//x方向を基準とする
		x[3]=PX[i+1];	f[3]=PF[i+1];
		x[4]=PX[i+2];	f[4]=PF[i+2];
		add_point(PX,PF,i);
		return 0;
	}
}

void akima::set(int px[],int pf[],int n)
{
	for(int i=0;i<n;i++){
		PX[i]=px[i];
		PF[i]=pf[i];
	}
	N=n;
}

void akima::draw(HDC hdc,BOOL maru)		//第2引数は○を描画するか否か
{
	int i,j,k;
	HPEN hPen,hOldPen;
	HBRUSH hBrush,hOldBrush;

	for(i=0;i<N-2;i++){
		k=set_point(i);		//計算用の配列に分点の座標を設定　基準方向判定
		if(x[2]<x[3])
			for(j=(int)x[2];j<=(int)x[3];j++)
				if(j==(int)x[2])
					if(k)	MoveToEx(hdc,(int)P(2,j),j,NULL);		//y方向を基準とする
					else	MoveToEx(hdc,j,(int)P(2,j),NULL);		//x方向を基準とする
				else
					if(k)	LineTo(hdc,(int)P(2,j),j);
					else	LineTo(hdc,j,(int)P(2,j));
		else
			for(j=(int)x[2];j>=(int)x[3];j--)		//大小符号に注目
				if(j==(int)x[2])
					if(k)	MoveToEx(hdc,(int)P(2,j),j,NULL);
					else	MoveToEx(hdc,j,(int)P(2,j),NULL);
				else
					if(k)	LineTo(hdc,(int)P(2,j),j);
					else	LineTo(hdc,j,(int)P(2,j));
	}
	if(i==N-2 && N>=3)		//最後だけ
		if(x[3]<x[4])
            for(j=(int)x[3];j<=(int)x[4];j++)
				if(k)	LineTo(hdc,(int)P(3,j),j);
				else	LineTo(hdc,j,(int)P(3,j));
		else
            for(j=(int)x[3];j>=(int)x[4];j--)
				if(k)	LineTo(hdc,(int)P(3,j),j);
				else	LineTo(hdc,j,(int)P(3,j));

	if(maru && N>=3){		//計算に使った分点を表示
		hPen=CreatePen(PS_SOLID,1,RGB(0,255,0));
		hOldPen=(HPEN)SelectObject(hdc,hPen);
 		hBrush=(HBRUSH)GetStockObject(NULL_BRUSH);
		hOldBrush=(HBRUSH)SelectObject(hdc,hBrush);
        for(i=0;i<7;i++)
			if(k)	Ellipse(hdc,(int)f[i]-10,(int)x[i]-10,(int)f[i]+10,(int)x[i]+10);
			else	Ellipse(hdc,(int)x[i]-10,(int)f[i]-10,(int)x[i]+10,(int)f[i]+10);
		DeleteObject(hPen);
		SelectObject(hdc,hOldPen);
		DeleteObject(hBrush);
		SelectObject(hdc,hOldBrush);
	}
}

void akima::force(int px[],int pf[],int *n)		//不正な分点は無効とする
{
	int xs,xc,xe,ys,yc,ye;

	if(*n>=3){
		xs=px[*n-3];	xc=px[*n-2];	xe=px[*n-1];
		ys=pf[*n-3];	yc=pf[*n-2];	ye=pf[*n-1];
 		if(xs>xe){xs=xe-xs;	xe-=xs;	xs+=xe;}		//2値交換
		if(ys>ye){ys=ye-ys;	ye-=ys;	ys+=ye;}		//xs<xe , ys<ye とする
		if((xs<xc && xc<xe) || (ys<yc && yc<ye));		//正常描画可能範囲
		else	--*n;
	}else if(*n>=2)
		if(px[*n-2]==px[*n-1] && pf[*n-2]==pf[*n-1])
			--*n;
}