// ■陰線消去
//   �@左クリックで図形が変わります。
//   �A右クリックで陰線消去の手法が変わります。
//
#include "myWin.h"
#include "math.h"
#define numX 51
#define numY 51
static HDC hBuff; static HBITMAP hBM;
double PI=3.14159265358979;     // π
double dlx, dx, dy;	// 表示刻み幅、XY軸方向の単位メッシュ幅
double alpha,beta;	// X(α)Y(β)軸との水平軸との角度
double Z [numX+1][numY+1];	// 高さデータ
double dZ[numX+1][numY+1];	// 傾斜分布
double *Ymax, *Ymin;// 上下浮動水平線
double Xlen, Ylen;  //　表示上のXY方向長さ
double curX, curY;	//　現在ペン位置
double dxCosA;		// dx*cos(α)
double dyCosB;		// dy*cos(β)
double dxSinA;		// dx*sin(α)
double dySinB;		// dy*sin(β)
double dlxTanA;		// dlx*tan(α)
double dlxTanB;		// dlx*tan(β)
int NR;				// 浮動水平線用配列の大きさ
int cnt=100;		// 形状変更用
int method=1;		// 手法変更用
int softFocus=true; // softFocusを行うかどうかのフラグ(行わないときfalseを指定する)
static double (*fun[10])(double A, double B);
void softFocusing(){
	if(!softFocus) return;
	int Px[9],R,G,B;
	for(int X=20;X<400;X++)for(int Y=20;Y<400;Y++){
		Px[0]=GetPixel(hBuff,X-1,Y-1);Px[1]=GetPixel(hBuff,X,Y-1);Px[2]=GetPixel(hBuff,X+1,Y-1);
		Px[3]=GetPixel(hBuff,X-1,Y  );Px[4]=GetPixel(hBuff,X,Y  );Px[5]=GetPixel(hBuff,X+1,Y  );
		Px[6]=GetPixel(hBuff,X-1,Y+1);Px[7]=GetPixel(hBuff,X,Y+1);Px[8]=GetPixel(hBuff,X+1,Y+1);
		R=Px[4]& 0xFF; G=(Px[4]>>8)& 0xFF; B=(Px[4]>>16) & 0xFF;;
		for(int i=0;i<9;i++){
			R+= Px[i]& 0xFF; G+= (Px[i]>>8)& 0xFF;  B+= (Px[i]>>16) & 0xFF;
		}
		R /= 10; G /=10; B /=10; SetPixel(hBuff,X,Y,RGB(R,G,B));
	}
}
void setDZ(){
	for(int i=0;i<=numX; i++) dZ[i][0]=0;
	for(int i=0;i<=numY; i++) dZ[0][i]=0;
	for(int i=1;i<=numX; i++) for(int j=1;j<=numX; j++){
		dZ[i][j]=(abs(Z[i][j]-Z[i][j-1])+abs(Z[i][j]-Z[i-1][j]))/40;
	}
}
void drawLine(double X1, double Y1, double X2, double Y2, int width, int color){
	HPEN hpen=CreatePen(PS_SOLID,width, color); SelectObject(hBuff, hpen);
	MoveToEx(hBuff,20+(int) (X1*2), 400 - (int)(Y1*2), NULL);
	LineTo  (hBuff,20+(int) (X2*2), 400 - (int)(Y2*2));
	DeleteObject(hpen);
}
int hiddenDraw(double px, double py,int p, int Visible, int Update){
	int R=false;
	if((py>=Ymax[p]) ||(py<=Ymin[p])){
		if(Update && (py>=Ymax[p])) Ymax[p]=py;
		if(Update && (py<=Ymin[p])) Ymin[p]=py;
		if(Visible) drawLine(curX, curY, px, py, 1, 0xFF0000);
		R = true;
	}
	curX=px; curY=py;return R;
}
double GroundX(int j, int k){// X0=Y0=0のときのX座標
	return ((double)j * dxCosA-(double)k * dyCosB);
}
double GroundY(int j, int k){// X0=Y0=0のときのX座標
	return ((double)j * dxSinA+(double)k * dySinB);
}
int DrawPos(int j, int k){// 表示位置
	return (int)(0.5+(Ylen+GroundX(j,k)/dlx));
}
double setPx(int j, int k,double X0, double PH){
	return PH* dlx + GroundX(j,k)+X0;
}
double setPy(int j, int k,double Y0, double PH, double fp){
	return PH* dlxTanA + GroundY(j,k) + fp + Y0;
}
double setPy2(int j, int k,double Y0, double PH, double fp){
	return -PH* dlxTanA + GroundY(j,k) + fp + Y0;
}
double xyExpM(double X, double Y){//高さの値を決める関数 (x-y)*Exp(-(x*x+y*y))
	X /= 4; Y /= -4;
	return (X-Y)*exp(-(X*X+Y*Y))*1.5;
}
double hSphere(double X, double Y){//高さの値を決める関数 球
	double R=sqrt(X*X+Y*Y);
	if(R<=7) return sqrt(49-R*R)*0.2;
	else return 0;
}
double mntFuji(double X, double Y){//高さの値を決める関数 富士山モデル
	double R=sqrt(X*X+Y*Y);
	if(R<=1) return 1;
	R--;return exp(-R/3);
}
double sinByR(double X, double Y) { //高さの値を決める関数 sin(R)/R
	double R = sqrt(X * X + Y * Y);
	if(R==0.0) return 1; else return sin(R)/R;
}
double cone(double X,double Y){ //　高さの値を決める関数 円錐モデル
	double R = sqrt(X * X + Y * Y);
	if(R<5) return 1-R*2/5; else return -1.0;
}
void setData(double (*fun)(double A, double B)){//データを設定
	double DNX2=(double) numX / 2.0, DNY2=(double) numY / 2.0, X;
	for(int j=0;j<numX; j++){
		X=0.3 * ((double)j -DNX2);
		for(int k=0;k<numY;k++)
			Z[j][k]=40.0*(*fun)(X,0.3 * ((double)k - DNY2));
	}
}
void initialize(){
	fun[0]= sinByR; fun[1] = xyExpM; fun[2] =hSphere;
	fun[3] =mntFuji; fun[4] =cone;
	dlx = 0.1; alpha = PI / 12; beta = PI / 8; dx = 2; dy = 1.4;
	dxCosA  = dx * cos(alpha)  ; dyCosB  = dy * cos(beta);
	dxSinA  = dx * sin(alpha)  ; dySinB  = dy * sin(beta);
	Xlen    = (numX-1) * dxCosA; Ylen    = (numY-1) * dyCosB; 
    dlxTanA = dlx * tan(alpha) ; dlxTanB = dlx * tan(beta);
	NR=(int)((Xlen+Ylen)/dlx)+2;
	Ymax=(double *)malloc(sizeof(double)*NR);
	Ymin=(double *)malloc(sizeof(double)*NR);
	curX=0; curY=0;
}
double getZ(int j, int k){
	if(Z[j][k]>=-999) return Z[j][k];
	double R=1000, N=0;
	for(int jj = j - 1; jj <= j + 1; jj += 2) 
		for(int kk = k - 1; kk <= k + 1; kk += 2){
		if(Z[jj][kk]>-999 && Z[jj][kk]<R) R=Z[jj][kk];
	}
	return R;
}

int screenX(int j, int k){
    return (int)((0.5 + (Ylen + GroundX(j, k)) / dlx));
}

int setColor(double Z[], int N, double CC){// もっと良い色設定があるかもしれない。
		double W=127; double SZ=0; for(int i=0; i<N; i++)SZ += Z[i];
	int R=(int)((SZ/160 + 1.001)*256);
	if(R>=512) return RGB(255        , (int)(CC*W), (int)(CC*W));
	if(R>=384) return RGB(255        , 511- R * 2 , (int)(CC*W));
	if(R>=256) return RGB((R-256)*2  , 255        , (int)(CC*W));
	if(R>=128) return RGB((int)(CC*W), R          , (int)(CC*W));
	if(R<0)    return RGB( 0         , 0          , 127        );
	           return RGB((int)(CC*R), R          , 255-R*2    );	

}
int setColor2(double Z[], int N, double CC){
	double W=127; double SZ=0; for(int i=0; i<N; i++)SZ += Z[i];
	int R=(int)((SZ/160 + 1.001)*256);
	if(R>=512) return RGB(127          , 0                , 0);
	if(R>=384) return RGB((int)(CC*127), (int)(CC*(511-R)), 0);
	if(R>=256) return RGB((int)(CC*(R-256)), (int)(CC*128), 0);
	if(R>=128) return RGB(0                , (int)(CC*(R-128))  ,(int)(CC*(256-R)));
	if(R<0)    return RGB((int)(CC*127)      , 0            , (int)(CC*128));
	           return RGB((int)(CC*(128-R))      , 0 , (int)(CC*(256-R)));	

}
int setColor3(double Z[], int N, double CC){
	double W=127; double SZ=0; for(int i=0; i<N; i++)SZ += Z[i];
	int R=(int)((SZ/160 + 1.001)*256);
	if(R>=512) return RGB(255          , 0                , 0);
	if(R>=384) return RGB((int)(CC*255), (int)(CC*(511-R)*2), 0);
	if(R>=256) return RGB((int)(CC*(R-256)*2), (int)(CC*255), 0);
	if(R>=128) return RGB(0                , (int)(CC*(R-128)*2)  ,(int)(CC*(256-R)*2));
	if(R<0)    return RGB((int)(CC*255)      , 0            , (int)(CC*255));
	           return RGB((int)(CC*(128-R)*2)      , 0 , (int)(CC*(R*2)));	

}
int setColor(double Z[], int N){
	double SZ=0; for(int i=0; i<N; i++)SZ += Z[i];
	int R=(int)((SZ/(N*40) + 1.001)*256);
	if(R>=512) return RGB(255      , 0         , 0  );
	if(R>=384) return RGB(255      , 511- R * 2, 0  );
	if(R>=256) return RGB((R-256)*2, 255       , 0  );
	if(R>=128) return RGB( 0       , R         , 0  );
	if(R<0)    return RGB( 0       , 0         , 64);
	return RGB( 0       ,    R   , 177-R);	
}

double cmpArea(double AX[],double AY[], int N){ //面積計算
	double A=0, X=(double)AX[0], Y=(double)AY[0];
	for(int i=1;i<N;i++){ A += X*AY[i] - Y*AX[i]; X=AX[i]; Y=AY[i];}
	return abs(A);	
}
void drawPoly(double X[], double Y[],int N, int color){
	POINT p[5];
	for(int i=0;i<N;i++) {p[i].x= 20+(int)(X[i]*2+0.5);p[i].y=400-(int)(Y[i]*2+0.5);}
	HBRUSH brush=CreateSolidBrush(color); SelectObject(hBuff,brush);
	Polygon(hBuff, p,N);DeleteObject(brush);
}
void setSquareXYZ(double X0, double Y0, int j, int k, double AX[], double AY[], double Z[]){
	int jj=j-1, kk=k-1;
	Z[0]=getZ(j,k), Z[1]=getZ(j,kk),Z[2]=getZ(jj,kk),Z[3]=getZ(jj,k);
	AX[0]=setPx(j , k , X0, 0); AY[0]=setPy(j , k , Y0, 0, Z[0]);
	AX[1]=setPx(j ,kk , X0, 0); AY[1]=setPy(j , kk, Y0, 0, Z[1]);
	AX[2]=setPx(jj,kk , X0, 0); AY[2]=setPy(jj, kk, Y0, 0, Z[2]);
	AX[3]=setPx(jj,k  , X0, 0); AY[3]=setPy(jj, k , Y0, 0, Z[3]);
}
void hiddenLine2(){//傾斜の度合いで明るさを変える（線を描く）（プロッタでは使えない）
	HPEN pen=CreatePen(PS_SOLID, 1,0x000000); SelectObject(hBuff,pen);
	double X0 = 80, Y0=100, maxS=0, AX[4], AY[4], AZ[4];
	// 奥方向から四角形を描く
	for(int k = numY-1; k >0; k--) for(int j = numX-1; j >0; j--)
		if(Z[j][k]>-999 || Z[j-1][k]>-999 || Z[j][k-1]>-999 || Z[j-1][k-1]>-999){
			setSquareXYZ(X0, Y0, j, k, AX, AY, AZ);
			//--------------------------------------------------------
			// 四角形で表示する方法
			//drawPoly(AX,AY,4, setColor(AZ, 4));
			//--------------------------------------------------------
			// 三角形に分ける方法（その１）
			//double TX[3], TY[3], TZ[3];
			//TX[0]=AX[0]; TY[0]=AY[0]; TZ[0]=AZ[0];
			//TX[1]=AX[2]; TY[1]=AY[2]; TZ[1]=AZ[2];
			//TX[2]=AX[1]; TY[2]=AY[1]; TZ[2]=AZ[1];
			//drawPoly(TX,TY,3, setColor(TZ, 3));
			//TX[2]=AX[3]; TY[2]=AY[3]; TZ[2]=AZ[3];
			//drawPoly(TX,TY,3, setColor(TZ, 3, CC));
			//--------------------------------------------------------
			// 三角形に分ける方法（その２）
			double TX[3], TY[3], TZ[3];
			double CX = 0, CY = 0, CZ=0;	
			for(int i=0;i<4; i++) { CX+=AX[i]; CY+=AY[i]; CZ+=AZ[i];}
			CX /= 4; CY /= 4; CZ /=4;
			TX[0]=CX; TY[0]=CY; TZ[0]=CZ;
			int IB=3;
			SelectObject(hBuff,GetStockObject(NULL_PEN));
			double CC=1-dZ[k][j];if(CC>1) CC=1;
	
			for(int i=0; i<4;i++) {
				TX[1]=AX[IB]; TX[2]= AX[i];
				TY[1]=AY[IB]; TY[2]= AY[i];
				TZ[1]=AZ[IB]; TZ[2]= AZ[i];
				drawPoly(TX,TY,3, setColor3(TZ, 3,CC));
				IB=i;
			}
			SelectObject(hBuff,pen);SetROP2(hBuff,R2_MASKPEN);// 画面とAndで線を描く
			drawPoly(AX,AY,4, 0xFFFFFF);SetROP2(hBuff,R2_COPYPEN);
		}
	DeleteObject(pen);
	softFocusing();
}
void hiddenLine3(){//傾斜の度合いで明るさを変える（線を描かない）
	SelectObject(hBuff,GetStockObject(NULL_PEN));
	double X0 = 80, Y0=100, AX[4], AY[4], AZ[4], CC;
	// 表示面積で色調を変える場合、最大値を求める
	//double maxS=0;
	//for(int k = numY-1; k >0; k--)	for(int j = numX-1; j >0; j--)
	//	if(Z[j][k]>-999 || Z[j-1][k]>-999 || Z[j][k-1]>-999 || Z[j-1][k-1]>-999){
	//		setSquareXYZ(X0, Y0, j, k, AX, AY, AZ);
	//		CC= cmpArea(AX,AY,4);if(maxS<CC)maxS=CC;
	//	}
	// 奥方向から四角形を描く
	for(int k = numY-1; k >0; k--) for(int j = numX-1; j >0; j--)
		if(Z[j][k]>-999 || Z[j-1][k]>-999 || Z[j][k-1]>-999 || Z[j-1][k-1]>-999){
			setSquareXYZ(X0, Y0, j, k, AX, AY, AZ);
			//CC= cmpArea(AX,AY,4)/maxS;if(CC>1)CC=1;
			//--------------------------------------------------------
			// 四角形で表示する方法
			//drawPoly(AX,AY,4, setColor(AZ,4, CC));
			//--------------------------------------------------------
			// 三角形に分ける方法（その１）
			//double TX[3], TY[3], TZ[3];
			//drawPoly(AX,AY,4, setColor(AZ,4, CC));
			//TX[0]=AX[0]; TY[0]=AY[0]; TZ[0]=AZ[0];
			//TX[1]=AX[2]; TY[1]=AY[2]; TZ[1]=AZ[2];
			//TX[2]=AX[1]; TY[2]=AY[1]; TZ[2]=AZ[1];
			//drawPoly(TX,TY,3, setColor(TZ, 3));
			//TX[2]=AX[3]; TY[2]=AY[3]; TZ[2]=AZ[3];
			//drawPoly(TX,TY,3, setColor(TZ,3,CC));
			//--------------------------------------------------------
			// 三角形に分ける方法（その２）
			double TX[3], TY[3], TZ[3];
			double CX = 0, CY = 0, CZ=0;	
			for(int i=0;i<4; i++) { CX+=AX[i]; CY+=AY[i]; CZ+=AZ[i];}
			CX /= 4; CY /= 4; CZ /=4;
			TX[0]=CX; TY[0]=CY; TZ[0]=CZ;
			CC=1-dZ[k][j]*1.5;if(CC>1) CC=1;if(CC<0.02) CC=0.02;
			int IB=3;
			for(int i=0; i<4;i++) {
				TX[1]=AX[IB]; TX[2]= AX[i];
				TY[1]=AY[IB]; TY[2]= AY[i];
				TZ[1]=AZ[IB]; TZ[2]= AZ[i];
				drawPoly(TX,TY,3, setColor3(TZ,3,CC));
				IB=i;
			}
		}
		softFocusing();
}

void hiddenLine(){
	for(int j=0; j<NR; j++)			//浮動水平線初期化
		{Ymax[j] = -1E20; Ymin[j] = 1E20;}	
	double X0 = 80, Y0=100;			//表示始点位置
	int Vflag, p1, p2;double H, DH, PH, px, py, fp;
	for(int k = 0; k < numY; k++){
		Vflag=false;
		for(int j = 0; j < numX - 1; j++) // X軸方向描画
			if(Z[j][k]>-999 || Z[j+1][k]>-999){
				double Z1=getZ(j,k), Z2=getZ(j+1,k);
				p1 = screenX(j,k); p2=screenX(j+1,k);
				H = Z1; DH = (Z2 - Z1) * dlx/ dxCosA;
				for(int p = p1; p <= p2; p++){
					PH=(double)(p - p1); fp = H + DH * PH;
					px = setPx(j, k, X0, PH); py = setPy(j, k, Y0, PH, fp);
					if((j < numX - 2 && p < p2) || (j == numX-2))
						Vflag=hiddenDraw(px,py,p,Vflag, true);
				}
			}
		if(k < numY - 1)    
			for (int j = 0; j < numX ; j++)    // Y軸方向描画
				if(Z[j][k]>-999 || Z[j][k+1]>-999){
					double Z1=getZ(j,k), Z2=getZ(j,k+1);
					Vflag = false;p1=screenX(j,k);p2=screenX(j,k+1);
					for (int p = p1; p >= p2; p--){               // 補間
						PH = (double)(p - p1);
						fp = Z1 - (Z2 - Z1) * PH * dlx / dyCosB;
						px = setPx(j, k, X0, PH); py = setPy2(j, k, Y0, PH, fp);
						Vflag= hiddenDraw(px, py, p, Vflag, p != p2);
					}
				}
	}
	softFocusing();
}
void procPaint(HWND hw, WPARAM wp,LPARAM lp){
	PAINTSTRUCT ps;	HDC hdc=BeginPaint(hw, &ps);
	BitBlt(hdc, 0,0,500,500,hBuff,0,0,SRCCOPY);
	EndPaint(hw,&ps);
}
void crearBitmap(int width, int height){
	SelectObject(hBuff,GetStockObject(NULL_PEN));
	PatBlt(hBuff,0, 0,width, height, WHITENESS);
}
void procLButtonDown(HWND hw, WPARAM wp,LPARAM lp){
	crearBitmap(500,500);
	cnt++;if(cnt>4)cnt=0;
	setData(fun[cnt]);setDZ();
	if     (method==0){SetWindowText(hw,TEXT("陰線消去 (その１）"));hiddenLine(); }
	else if(method==1){ SetWindowText(hw,TEXT("陰線消去 (その２）")); hiddenLine2();}
	else { SetWindowText(hw,TEXT("陰線消去 (その３）線を描かなければ…")); hiddenLine3();}
	InvalidateRect(hw,NULL,FALSE);
}
void procRButtonDown(HWND hw, WPARAM wp,LPARAM lp){
	method++;if(method>2) method=0;
	cnt--; procLButtonDown(hw, wp,lp);
}
void genBitmap(HWND hw, int width,int height){
	HDC hdc=GetDC(hw); hBM=CreateCompatibleBitmap(hdc,500,400);
	hBuff=CreateCompatibleDC(hdc); SelectObject(hBuff,hBM);
	crearBitmap(500,500); ReleaseDC(hw,hdc);
}
void setWindowSize(HWND hw, int width,int height){
	RECT r;	GetWindowRect(hw, &r);
	MoveWindow(hw,r.left,r.top, width,height,TRUE);
}
void procCreate(HWND hw, WPARAM wp,LPARAM lp){
	setWindowSize(hw,452,400);	genBitmap(hw, 500,500);
	initialize(); 
	procLButtonDown(hw, wp,lp);
}
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;
	case WM_LBUTTONDOWN   : procLButtonDown(hw,wp,lp)  ; break;
	case WM_RBUTTONDOWN   : procRButtonDown(hw,wp,lp)  ; break;
	}
	return DefWindowProc(hw,msg,wp,lp);
}