// ■簡便法による３次スプライン曲線の描画　Coded by Yutaka Shirai.
//   
//   詳しい式の説明については、拙著「ExcelとVBAによる実用数値解析入門」の
//   「5.4 スプライン補間 (6)もうひとつの簡便法」を参照してください。
//
//  【自然条件の場合の解くべき式】
//  ┌                                                                           ┐┌     　┐
//　│2σ0     σ0       0          0    …　 …　       …　               0　  ││ P0' 　│
//  │ σ1  2(σ0+σ1)  σ0         0    …　 …　       …　               0　  ││ P1' 　│
//  │  0      σ2     2(σ1+σ2)  σ1   …　 …　       …　               0　  ││ P2' 　│
//  │  :      ：　　　　：　　　　：    ：   ：         ： 　             ：    ││ ：    │
//  │  0       0        …　　　  …　  …　σ(n-2)  2(σ(n-3)+σ(n-2))  σ(n-3)││P(n-2)'│
//　│  0       0        …　　　  …　  …   0         σ(n-1)          2σ(n-1)││P(n-1)'│
//　└                                                                           ┘└       ┘
//    ┌                                                                   ┐
//　  │3(P1-P0)　                                                         │
//    │3{(P1-P0)(σ1/σ0)+(P2-P1)(σ0/σ1)}                               │
//  = │3{(P2-P1)(σ2/σ1)+(P3-P2)(σ1/σ2)}                               │
//    │  :      ：　　　　：　　　　：    ：   …         …　        0　 │
//    │3{(P(n-2)-P(n-3))(σ(n-2)/σ(n-3)+(P(n-1)-P(n-2))(σ(n-3)/σ(n-2)} │
//　  │3{(P(n-1)-P(n-2))                                                  │
//　  └                                                                   ┘

#include "myWin.h"
#include "stdlib.h"
static HDC hbuf; static HBITMAP hBM;                 //描画用ビットマップ
double currentX, currentY;int Frag;                  //直線描画座標値（プログラムで設定される） 
//--------設定を変えるときは以下の変更してください---------------------------------------------
#define NUMDATA 4                                    //通過点座標の数
double pnt[NUMDATA][2]={{0,0}, {0,2}, {2,0}, {2,2}}; //通過点座標（X,Y）の順に並べる
double Sgm[NUMDATA-1]={1,1,1};                       //σ(k+1)：点幅が極端に異なるときの調整用(通常1)
double ds=0.05;                                      //計算刻み幅
double P0d[2]={1,1};                                 //固定条件のときのP0'
double Pnd[2]={1,1};                                 //固定条件のときのPn'
int Cond=true;                                       //固定条件のときtrue, 自然条件のときFalse
//--------設定変更 ここまで---------------------------------------------------------------------
// ■P0，Pn 以外の 右辺の計算
double func(double P1,double P2, double P3,double sgm1, double sgm2){
  return 3 * ((P2 - P1) * sgm2 / sgm1 + (P3 - P2) * sgm1 / sgm2);
}
// ■三重対角方程式の解法
void TriDiag(double A[][3], double D[], double X[], int N){
   double CC[NUMDATA];
   double btemp = A[0][1]; X[0] = D[0] / btemp;
   for(int j=1;j<N;j++){
     CC[j] = A[j - 1][2] / btemp;
	 btemp = A[j][1] - A[j][0] * CC[j];
     X[j] = (D[j] - A[j][0] * X[j - 1]) / btemp;
   }
   for(int j=N-2;j>=0; j--) X[j] = X[j] - CC[j + 1] * X[j + 1];
}
// ■固定条件(P0'，Pn'を固定的に与える）
void fixCondition(double A[][3], double sgm[], double Pd[][NUMDATA], double Dd[][NUMDATA], int N){
   double Dtemp[NUMDATA], Ptemp[NUMDATA];//A[0][0], A[N-1][N-1]の値は意味がない。
   A[0    ][1] = 2 * (sgm[0] + sgm[1]);  A[0    ][2] = sgm[0];
   A[N - 3][0] = sgm[N - 2]           ;  A[N - 3][1] = 2 * (sgm[N - 3] + sgm[N - 2]);
   for(int k=2;k<N-2; k++){
     A[k - 1][0] = sgm[k];
     A[k - 1][1] = 2 * (sgm[k - 1] + sgm[k]);
     A[k - 1][2] = sgm[k - 1];
   }
   for(int j=0;j<2;j++){
     for(int k=1;k<N-1;k++){
       Dd[j][k-1] = func(pnt[k-1][j], pnt[k][j], pnt[k+1][j], sgm[k-1], sgm[k]);
	 }
     Dd[j][0  ] = Dd[j][0  ] - sgm[1  ] * P0d[j];
     Dd[j][N-3] = Dd[j][N-3] - sgm[N-3] * Pnd[j];
	 for(int k=0;k<N;k++) Dtemp[k] = Dd[j][k];
	 TriDiag(A, Dtemp,Ptemp,N - 2);                    
	 for(int k=0;k<N-1;k++) Pd[j][k + 1] = Ptemp[k];   
	 Pd[j][0] = P0d[j];  Pd[j][N-1] = Pnd[j]; 
   }
}
// ■自然条件(P0"=Pn"=0 : 始点・終点の二次微分値を0とする。このときP0'とP1'は意味をもたない）
void naturalCondition(double A[][3], double sgm[], double Pd[][NUMDATA], double Dd[][NUMDATA], int N){
   double Dtemp[NUMDATA], Ptemp[NUMDATA];//A[0][0], A[N-1][2]の値は意味がない。
    A[0  ][1] = 2 * sgm[1]    ; A[0  ][2] = sgm[1];
   A[N-1][0] = sgm[N - 2]; A[N-1][1] = 2 * sgm[N - 2]; A[N-1][2] = 0;
   for(int k=1;k<N-1;k++){
     A[k][0] = sgm[k];
     A[k][1] = 2 * (sgm[k - 1] + sgm[k]);
     A[k][2] = sgm[k - 1];
   }
   for(int j=0;j<2;j++){
     Dd[j][0  ] = 3.0 * (pnt[1  ][j] - pnt[0    ][j]);
     Dd[j][N-1] = 3.0 * (pnt[N-1][j] - pnt[N - 2][j]);
     for(int k=1;k<N-1;k++){
       Dd[j][k] = func(pnt[k - 1][j], pnt[k][j], pnt[k + 1][j], sgm[k - 1], sgm[k]);
     }
	 for(int k=0;k<N;k++) Dtemp[k] = Dd[j][k];
     TriDiag(A, Dtemp, Ptemp, N);
	 for(int k=0;k<N;k++) Pd[j][k] = Ptemp[k];
   }
}
// ■線の描画
void Draw(double X, double Y, int *Flag, int Color, int w){
  if(*Flag){
      int X1 = 50 + (int)(currentX * 100),Y1 = 250 - (int)(currentY * 100);
      int X2 = 50 + (int)(X * 100),  Y2 = 250 -(int)(Y * 100);
	  HPEN pen=CreatePen(PS_SOLID, w, Color);SelectObject(hbuf,pen);
	  MoveToEx(hbuf, X1,Y1,NULL);LineTo(hbuf,X2,Y2);
	  DeleteObject(pen);
  }
  currentX = X; currentY = Y;
  *Flag = true;
}
// ■スプライン曲線の計算と描画
void Spline(double pnt[][2], double sgm[], double P0d[], double Pnd[], 
	double ds, int N, int Cond){
  double Ak0[2], Ak1[2], Ak2[2], Ak3[2];
  double A[NUMDATA][3], Pd[2][NUMDATA],Dd[2][ NUMDATA];
  if(Cond) fixCondition    (A, sgm, Pd, Dd, N);    //固定条件
  else     naturalCondition(A, sgm, Pd, Dd, N);    //自然条件
  int Flag = false;
  Draw(pnt[0][0], pnt[0][1], &Flag,0xFF0000,2);    //初期位置に移動
  for(int k=0; k<N-1; k++){                        //各区間の計算と描画
  	for(int j=0;j<2;j++){// j=0 : X計算用,  j=1 : Y計算用 
		double sgm2 = sgm[k] * sgm[k],  sgm3 = sgm[k] * sgm2;
		Ak0[j] = 2.0 * (pnt[k][j] - pnt[k + 1][j]) / sgm3 +
			(Pd[j][k] + Pd[j][k + 1]) / sgm2;
		Ak1[j] = 3.0 * (pnt[k + 1][j] - pnt[k][j]) / sgm2 -
			(2.0 * Pd[j][k] + Pd[j][k + 1]) / sgm[k];
		Ak2[j] = Pd[j][k]; Ak3[j] = pnt[k][j];
	}
   	double sk = 0;// 1区間の描画(刻み幅 ds)
	while((sk + ds) <= sgm[k]){//
      sk += ds;
      double PX = ((Ak0[0] * sk + Ak1[0]) * sk + Ak2[0]) * sk + Ak3[0];
      double PY = ((Ak0[1] * sk + Ak1[1]) * sk + Ak2[1]) * sk + Ak3[1];
      Draw (PX, PY, &Flag,0xFF0000,2);
	}
  }
  Draw(pnt[N-1][0],pnt[N-1][1], &Flag,0xFF0000,2);
}
// ■ビットマップ作成とメモリデバイスコンテキスト作成
void initBitMap(HWND hw){
	HDC hdc=GetDC(hw); hBM=CreateCompatibleBitmap(hdc,330,330);
	hbuf=CreateCompatibleDC(hdc);SelectObject(hbuf,hBM);
	SelectObject(hbuf,GetStockObject(NULL_PEN));
	PatBlt(hbuf,0,0,330,330,WHITENESS);
	ReleaseDC(hw, hdc);
}
// ■通過点を直線で結んだ図を描く
void drawLinear(){
	int Flag =false;
	Draw(pnt[0][0], pnt[0][1], &Flag, 0xFF,1);
    for(int k=0; k<NUMDATA;k++){
      Draw(pnt[k][0], pnt[k][1], &Flag,0xFF,1);
	}
}
void procCreate(HWND hw, WPARAM wp, LPARAM lp){
	SetWindowText(hw,TEXT("簡便法による３次スプライン曲線の描画"));
	initBitMap(hw);	MoveWindow(hw,0,0,330,330,TRUE);
	drawLinear();                                  //通過点を直線で結んだ図を描く
	Spline(pnt, Sgm, P0d, Pnd, ds, NUMDATA, Cond); //スプライン曲線を描く
	InvalidateRect(hw,NULL,FALSE);
}
void procPaint(HWND hw, WPARAM wp, LPARAM lp){// ビットマップを画面に表示
	PAINTSTRUCT ps; HDC hdc=BeginPaint(hw, &ps);
	BitBlt(hdc,0,0,330,330,hbuf,0,0,SRCCOPY);
	EndPaint(hw,&ps);
}
LRESULT CALLBACK WndProc(HWND hw, UINT msg, WPARAM wp, LPARAM lp){
	switch(msg){
	case WM_DESTROY : PostQuitMessage(0)    ; return 0;
	case WM_CREATE  : procCreate(hw,wp,lp)  ; return 0;
	case WM_PAINT   : procPaint (hw, wp, lp); return 0;
	}
	return DefWindowProc(hw, msg, wp, lp);
}