Ryoichi Mizuno - Scientific Computer Graphics | Upper Directory (English) | icon

//////////////////////////////////////////////////////////
/* Bezier Surface, Copyright 2002 Ryoichi Mizuno        */
/* ryoichi[at]mizuno.org                                */
/* Dept. of Complexity Science and Engineering          */
/* at The University of Tokyo                           */
//////////////////////////////////////////////////////////
spacer
/* 解説 */
// nxm次ベジエ曲線を描画する
// ベジエ曲面の次数および分割数を任意に設定することができる
// 視点および視野を任意に変更することができる
// 制御点をユーザの直感に基づき任意に移動することができる
// ソフトウェアにより透視投影および陰面(線)消去を行っている
// 陰面(線)消去には奥行ソート法を用いているので正常に描画されない形状もある

/* 詳細な解説など */
// <http://www.mizuno.org/gl/bs/>を参照

/* 既知のバグまたは不備 */
// 1. アニメーション中(animationFlag!=0)glui1またはglui2よりmain_windowにクリックアンドドラッグするとフォーカスがglui1またはglui2に移動する

/* 既知のバグまたは不備の対処法 */
// 1. アニメーション中はglui1またはglui2よりmain_windowにクリックアンドドラッグしない、フォーカスが移動してしまった場合ウィンドウを一旦最小化する

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <malloc.h>
#include <time.h>
#include <GL/glut.h>
#include <GL/glui.h>

struct vector2D // 2次元ベクトルの構造体
{
	float x;
	float y;
};

struct vector3D // 3次元ベクトルの構造体
{
	float x;
	float y;
	float z;
};

struct vector5D // 5次元ベクトルの構造体
{
	float x;
	float y;
	float z;
	float p; // 透視変換後のx座標
	float q; // 透視変換後のy座標
};

struct index2D // 2次元の配列のインデックスの構造体
{
	int k;
	int l;
};

struct normalVectorAndZvalue // 法線ベクトル、表裏を識別するフラッグ、奥行値の構造体
{
	struct index2D basicpoint; // 基準となる点のインデックス
	struct vector3D n; // 法線ベクトル
	unsigned int visible; // 表(1)裏(0)を識別するフラッグ
	float Z; // 奥行値
	float height; // 重心の高さ(z座標)
};

struct color
{
	float r;
	float g;
	float b;
};

/* 数学関係の定数 */
const float PI=3.141592853f;

/* ベジエ曲面関係の変数 */
float size=100.0f; // 基準となる長さ(固定)
int n=3; // uに関しての次数(n>0)
int m=3; // vに関しての次数(m>0)
int usubdivi=25; // u方向の分割数
int vsubdivi=25; // u方向の分割数
struct vector5D **ctrlPoint; // 制御点の座標のポインタ
struct vector5D **bezierPoint; // ベジエ曲面上の座標のポインタ
struct normalVectorAndZvalue *bezierSurface; // ベジエ曲面上の平面の法線ベクトルと奥行値のポインタ

/* 透視変換関係の変数 */
float theta=0.0f*PI/180.0f; // 方位角
float phi=0.0f*PI/180.0f; // 仰角
float d0=200.0f; // 視点から注視点までの距離の初期値
float d=d0; // 視点から注視点までの距離
struct vector3D eye; // 視点の座標
float mat[9]; // 変換行列
int omegadeg=25; // 画角(度数法)
float omegarad=(float)omegadeg*PI/180.0f; // 画角(弧度法)
float radius=50.0f; // 視野半径
float f=radius/(float)tan(omegarad/2.0f); // 視距離

/* 制御点の移動関係の変数 */
int ctrlPointSelectedFlag=0; // 選択された制御点の有無
struct index2D selectedCtrlPointIndex; // 選択された制御点のインデックス

/* フラッグ関係の変数 */
int showCtrlPointsFlag=1; // 制御点
int showBezierWiresFlag=1; // ベジエ曲線
int showBezierSurfacesFlag=1; // ベジエ曲面
int hiddenLineRemovalFlag=1; // 陰線消去
int sourceShapeFlag=0; // 元の形状 (0)凸曲面(1)正弦曲面(2)平面
int animationFlag=0; // アニメーションの種類(0)手動(1)θ回転(2)φ回転(3)θ+φ回転

/* ベンチマークのための変数 */
int counter=0;
float timer=0.0f;
float fps=0.0f;

/* OpenGL関係の変数 */
int main_window;
GLUI *glui1,*glui2;

// 動的な配列の作成の失敗時の処理
void processMemoryOver(void)
{
	printf("memory over\n");
	exit(0);
}

// 制御点の座標を格納する配列の作成
void makeCtrlPointsArray(void)
{
	int i;

	if((ctrlPoint=(struct vector5D **)malloc((n+1)*sizeof(struct vector5D *)))==NULL) processMemoryOver();
	for(i=0;i<n+1;i++)
	{
		if((ctrlPoint[i]=(struct vector5D *)malloc((m+1)*sizeof(struct vector5D)))==NULL) processMemoryOver();
	}
}

// ベジエ曲面上の点の座標を格納する配列の作成
void makeBezierPointsArray(void)
{
	int k;

	if((bezierPoint=(struct vector5D **)malloc(usubdivi*sizeof(struct vector5D *)))==NULL) processMemoryOver();
	for(k=0;k<usubdivi;k++)
	{
		if((bezierPoint[k]=(struct vector5D *)malloc(vsubdivi*sizeof(struct vector5D)))==NULL) processMemoryOver();
	}
}

// ベジエ曲面上の平面の法線ベクトルと奥行値を格納する配列の作成
void makeBezierSurfacesArray(void)
{
	if((bezierSurface=(struct normalVectorAndZvalue *)malloc((usubdivi-1)*(vsubdivi-1)*sizeof(struct normalVectorAndZvalue)))==NULL) processMemoryOver();
}

// 制御点の初期化
void processInitCtrlPoints(void)
{
	int i,j;
	float k,l; // [-0.5,+0.5]の変数

	// 制御点の初期座標を設定(等間隔)
	for(i=0;i<n+1;i++)
	{
		for(j=0;j<m+1;j++)
		{
			k=(float)i/(float)n-0.5f;
			l=(float)j/(float)m-0.5f;
			ctrlPoint[i][j].x=size*k; // [-0.5*size,+0.5*size]
			ctrlPoint[i][j].y=size*l; // [-0.5*size,+0.5*size]
			if(sourceShapeFlag==0)ctrlPoint[i][j].z=size*(k*k+l*l-0.25f)*2.0f; // 凸曲面 [-0.5*size,+0.5*size]
			else if(sourceShapeFlag==1)ctrlPoint[i][j].z=size*((float)sin(2.0f*PI*k)+(float)sin(2.0f*PI*l))/4.0f; // 正弦曲面 [-0.5*size,+0.5*size]
			else if(sourceShapeFlag==2)ctrlPoint[i][j].z=0.0f; // 平面
			//printf("ctrlPoint (world): %f %f %f\n",ctrlPoint[i][j].x,ctrlPoint[i][j].y,ctrlPoint[i][j].z);
		}
	}
}

// 透視変換(注視点は原点(0,0,0)に固定)
struct vector2D processPerspectiveTransformation(struct vector5D input)
{
	struct vector3D diff; // 視点とx, y, zとの差
	struct vector3D inte; // 計算の便宜上の中間的な座標 
	struct vector2D output; // 変換後の座標

	diff.x=input.x-eye.x; diff.y=input.y-eye.y; diff.z=input.z-eye.z;

	inte.x=mat[0]*diff.x+mat[1]*diff.y+mat[2]*diff.z;
	inte.y=mat[3]*diff.x+mat[4]*diff.y+mat[5]*diff.z;
	inte.z=mat[6]*diff.x+mat[7]*diff.y+mat[8]*diff.z;
	//printf("inte: %f %f %f\n",inte.x,inte.y,inte.z);

	output.x=f*inte.y/inte.x;
	output.y=f*inte.z/inte.x;

	return output;
}

// 制御点の移動に伴う逆透視変換(奥行値を固定)
void processInversePerspectiveTransformation(void)
{
	struct vector5D selectedCtrlPoint;
	struct vector3D diff; // 視点と制御点との差
	struct vector3D inte; // 計算の便宜上の中間的な座標 
	float coefficient; // 計算の便宜上の係数
	float Z; // 奥行値

	selectedCtrlPoint=ctrlPoint[selectedCtrlPointIndex.k][selectedCtrlPointIndex.l];
	// 奥行値の取得
	diff.x=eye.x-selectedCtrlPoint.x; diff.y=eye.y-selectedCtrlPoint.y; diff.z=eye.z-selectedCtrlPoint.z;
	Z=(float)sqrt(diff.x*diff.x+diff.y*diff.y+diff.z*diff.z);
	// 計算の便宜上の中間的な座標の取得
	inte.x=mat[0]+mat[3]*selectedCtrlPoint.p/f+mat[6]*selectedCtrlPoint.q/f;
	inte.y=mat[1]+mat[4]*selectedCtrlPoint.p/f+mat[7]*selectedCtrlPoint.q/f;
	inte.z=mat[2]+mat[5]*selectedCtrlPoint.p/f+mat[8]*selectedCtrlPoint.q/f;
	// 計算の便宜上の係数の取得
	coefficient=-Z/(float)sqrt(inte.x*inte.x+inte.y*inte.y+inte.z*inte.z);
	// x, y, z座標の計算
	inte.x*=coefficient; inte.y*=coefficient; inte.z*=coefficient;
	selectedCtrlPoint.x=inte.x+eye.x; selectedCtrlPoint.y=inte.y+eye.y; selectedCtrlPoint.z=inte.z+eye.z;
	ctrlPoint[selectedCtrlPointIndex.k][selectedCtrlPointIndex.l]=selectedCtrlPoint;
}

// 制御点の描画
void drawCtrlPoints(void)
{
	if(showCtrlPointsFlag)
	{
		int i,j;

		glColor3f(0.0f,0.0f,0.0f);
		// 点の描画
		glPointSize(5.0f);
		glBegin(GL_POINTS);
		for(i=0;i<n+1;i++)
			{
				for(j=0;j<m+1;j++)
				{
					glVertex2f(ctrlPoint[i][j].p,ctrlPoint[i][j].q);
				}
		}
		glEnd();
		//選択された制御点の描画
		if(ctrlPointSelectedFlag)
		{
			i=selectedCtrlPointIndex.k; j=selectedCtrlPointIndex.l;
			glColor3f(0.0f,1.0f,0.2f);
			glBegin(GL_POINTS);
			glVertex2f(ctrlPoint[i][j].p,ctrlPoint[i][j].q);
			glEnd();
			glColor3f(0.0f,0.0f,0.0f);
		}
		// 線の描画
		for(i=0;i<n+1;i++)
		{
			glBegin(GL_LINE_STRIP);
			for(j=0;j<m+1;j++)
			{
				glVertex2f(ctrlPoint[i][j].p,ctrlPoint[i][j].q);
			}
			glEnd();
		}
		for(j=0;j<(m+1);j++)
		{
			glBegin(GL_LINE_STRIP);
			for(i=0;i<n+1;i++)
			{
				glVertex2f(ctrlPoint[i][j].p,ctrlPoint[i][j].q);
			}
			glEnd();
		}
	}
}

// 規格化された値([0.0,1.0])からRGBを返す関数
struct color processStandardizedValue2RGB(float x)
{
	struct color y;
	float h,f,p,q,t;
	int n;

	// 輝度hをB?R(HSB色空間)で取得
	h=(1.0f-x)*240.0f;
	// hの範囲を[0.0,4.0)に規格化
	h/=60.0f;
	// hの最大整数を取得
	n=(int)floor(h);
	// hの小数部分を取得
	f=h-(float)n;
	p=0.0f; q=1.0f-(1.0f*f); t=1.0f-(1.0f-f);
	switch(n)
	{
		case 0:
			y.r=1.0f; y.g=t; y.b=p;
			break;
		case 1:
			y.r=q; y.g=1.0f; y.b=p;
			break;
		case 2:
			y.r=p; y.g=1.0f; y.b=t;
			break;
		case 3:
			y.r=p; y.g=q; y.b=1.0f;
			break;
		case 4:
			y.r=t; y.g=p; y.b=1.0f;
			break;
		case 5:
			y.r=1.0f; y.g=p; y.b=q;
			break;
	}

	return y;
}

// ベジエ曲面の描画
void drawBezierSurfaces(void)
{
	int k,l,f;
	float standardizedheight; // [0.0,1.0]の変数
	struct color surfacecolor;
	
	// ベジエ曲面上の平面の描画
	for(f=0;f<(usubdivi-1)*(vsubdivi-1);f++)
	{
		k=bezierSurface[f].basicpoint.k; l=bezierSurface[f].basicpoint.l;
		if(showBezierSurfacesFlag || hiddenLineRemovalFlag)
		{
			if(showBezierSurfacesFlag)
			{
				// 重心の高さ(z座標)を1に規格化
				standardizedheight=(bezierSurface[f].height+size*0.5f)/size;
				if(standardizedheight<0.0f) standardizedheight=0.0f;
				else if(standardizedheight>1.0f) standardizedheight=1.0f;
				if(bezierSurface[f].visible)
				{
					standardizedheight*=-1.0f; standardizedheight+=1.0f;
				}
				surfacecolor=processStandardizedValue2RGB(standardizedheight);
				glColor3f(surfacecolor.r,surfacecolor.g,surfacecolor.b);
			}
			else glColor3f(1.0f,1.0f,1.0f);
			// 面の塗りつぶし
			glBegin(GL_QUADS);
			glVertex2f(bezierPoint[k][l].p,bezierPoint[k][l].q);
			glVertex2f(bezierPoint[k+1][l].p,bezierPoint[k+1][l].q);
			glVertex2f(bezierPoint[k+1][l+1].p,bezierPoint[k+1][l+1].q);
			glVertex2f(bezierPoint[k][l+1].p,bezierPoint[k][l+1].q);
			glEnd();
		}
		if(showBezierWiresFlag)
		{
			// 辺の描画
			if(bezierSurface[f].visible)glColor3f(0.0f,0.0f,1.0f);
			else glColor3f(1.0f,0.0f,0.0f);
			glBegin(GL_LINE_LOOP);
			glVertex2f(bezierPoint[k][l].p,bezierPoint[k][l].q);
			glVertex2f(bezierPoint[k+1][l].p,bezierPoint[k+1][l].q);
			glVertex2f(bezierPoint[k+1][l+1].p,bezierPoint[k+1][l+1].q);
			glVertex2f(bezierPoint[k][l+1].p,bezierPoint[k][l+1].q);
			glEnd();
		}
	}
}

// 視点と変換行列の取得
void getEyeCoordinateAndTransformationMatrix(void)
{
	float st,ct,sp,cp; // 左から順にsin(θ), cos(θ), sin(φ), cos(φ)

	st=(float)sin(theta); ct=(float)cos(theta);
	sp=(float)sin(phi); cp=(float)cos(phi);

	eye.x=cp*ct*d; eye.y=cp*st*d; eye.z=-sp*d;
	//printf("eye: %f %f %f\n",eye.x,eye.y,eye.z);

	mat[0]=ct*cp; mat[1]=st*cp; mat[2]=-sp;
	mat[3]=-st; mat[4]=ct; mat[5]=0.0f;
	mat[6]=ct*sp; mat[7]=st*sp; mat[8]=cp;
}

// 制御点の透視変換後のx (p), y (q)座標の取得
void getCtrlPointsTransformatedCoordinate(void)
{
	int i,j;
	struct vector2D buffer; // 透視変換後のx (p), y (q)座標を一時的に格納する構造体

	for(i=0;i<n+1;i++)
		{
			for(j=0;j<m+1;j++)
			{
				// 透視変換
				buffer=processPerspectiveTransformation(ctrlPoint[i][j]);
				ctrlPoint[i][j].p=buffer.x; ctrlPoint[i][j].q=buffer.y;
			}
	}
}

// ベジエ曲面上の点のx, y, z座標の取得
void getBezierPointsCoordinate(void)
{
	int i,j,k,l;
	float u,v; // [0.0,1.0]の変数
	float *bu,*bv; // uとvに関するバーンスタイン関数のポインタ

	// バーンスタイン関数を格納する配列の作成
	if((bu=(float *)malloc((n+2)*sizeof(float)))==NULL) processMemoryOver();
	if((bv=(float *)malloc((m+2)*sizeof(float)))==NULL) processMemoryOver();

	for(k=0;k<usubdivi;k++)
	{
		u=(float)k/(float)(usubdivi-1);
		// uに関するバーンスタイン関数を0で初期化
		for(i=0;i<(n+2);i++) bu[i]=0.0f;
		// deCasteljauのアルゴリズムによりuに関するバーンスタイン関数を計算
		bu[1]=1.0f;
		for(i=1;i<n+1;i++)
		{
			for(j=i+1;j>0;j--)
			{
				bu[j]=(1-u)*bu[j]+u*bu[j-1];
			}
		}
		for(l=0;l<vsubdivi;l++)
		{
			v=(float)l/(float)(vsubdivi-1);
			// uに関するバーンスタイン関数を0.0で初期化
			for(j=0;j<(m+2);j++) bv[j]=0.0f;
			// deCasteljauのアルゴリズムによりvに関するバーンスタイン関数を計算
			bv[1]=1.0f;
			for(j=1;j<m+1;j++)
			{
				for(i=j+1;i>0;i--)
				{
					bv[i]=(1-v)*bv[i]+v*bv[i-1];
				}
			}
			// x, y, z座標を0.0で初期化
			bezierPoint[k][l].x=0.0f; bezierPoint[k][l].y=0.0f; bezierPoint[k][l].z=0.0f;
			for(i=0;i<n+1;i++)
			{
				for(j=0;j<m+1;j++)
				{
					// x, y, z座標の取得
					bezierPoint[k][l].x+=ctrlPoint[i][j].x*bu[i+1]*bv[j+1];
					bezierPoint[k][l].y+=ctrlPoint[i][j].y*bu[i+1]*bv[j+1];
					bezierPoint[k][l].z+=ctrlPoint[i][j].z*bu[i+1]*bv[j+1];
				}
			}
		}
	}
	// バーンスタイン関数を格納する配列の解放
	free(bu); free(bv);
}

// ベジエ曲面上の点の透視変換後のx (p), y (q)座標の取得
void getBezierPointsTransformatedCoordinate(void)
{
	int k,l;
	struct vector2D buffer; // 透視変換後のx (p), y (q)座標を一時的に格納する構造体

	for(k=0;k<usubdivi;k++)
		{
			for(l=0;l<vsubdivi;l++)
			{
				// 透視変換
				buffer=processPerspectiveTransformation(bezierPoint[k][l]);
				bezierPoint[k][l].p=buffer.x; bezierPoint[k][l].q=buffer.y;
			}
	}
}

// ベジエ曲面上の平面の法線ベクトルの取得
void getBezierSurfacesNormalVector(void)
{
	int k,l,f;
	struct vector3D p1,p2,p3;

	f=0;
	for(k=0;k<usubdivi-1;k++)
	{
		for(l=0;l<vsubdivi-1;l++)
		{
			// 配列のインデックスの格納
			bezierSurface[f].basicpoint.k=k; bezierSurface[f].basicpoint.l=l;

			// 法線ベクトルの格納
			p1.x=bezierPoint[k+1][l].x; p1.y=bezierPoint[k+1][l].y; p1.z=bezierPoint[k+1][l].z;
			p2.x=bezierPoint[k][l].x; p2.y=bezierPoint[k][l].y; p2.z=bezierPoint[k][l].z;
			p3.x=bezierPoint[k][l+1].x; p3.y=bezierPoint[k][l+1].y; p3.z=bezierPoint[k][l+1].z;
			bezierSurface[f].n.x=(p2.y-p1.y)*(p3.z-p2.z)-(p2.z-p1.z)*(p3.y-p2.y);
			bezierSurface[f].n.y=(p2.z-p1.z)*(p3.x-p2.x)-(p2.x-p1.x)*(p3.z-p2.z);
			bezierSurface[f].n.z=(p2.x-p1.x)*(p3.y-p2.y)-(p2.y-p1.y)*(p3.x-p2.x);
			f++;
		}
	}
}

// ベジエ曲面上の平面の表裏の識別及び奥行値の取得
void getBezierSurfacesZvalue(void)
{
	int k,l,f;
	struct vector3D gravityCenter;
	struct vector3D diff; // 視点と重心との差

	for(f=0;f<(usubdivi-1)*(vsubdivi-1);f++)
	{
		k=bezierSurface[f].basicpoint.k; l=bezierSurface[f].basicpoint.l;
		// 表(1)裏(0)の識別
		if((eye.x-bezierPoint[k][l].x)*bezierSurface[f].n.x+(eye.y-bezierPoint[k][l].y)*bezierSurface[f].n.y+(eye.z-bezierPoint[k][l].z)*bezierSurface[f].n.z>=0)
		{
			bezierSurface[f].visible=1;
		}
		else
		{
			bezierSurface[f].visible=0;
		}
		// 重心の座標の取得
		gravityCenter.x=(bezierPoint[k][l].x+bezierPoint[k+1][l+1].x)/2.0f;
		gravityCenter.y=(bezierPoint[k][l].y+bezierPoint[k+1][l+1].y)/2.0f;
		gravityCenter.z=(bezierPoint[k][l].z+bezierPoint[k+1][l+1].z)/2.0f;
		bezierSurface[f].height=gravityCenter.z;
		// 視点から重心までの距離(の自乗)の取得
		diff.x=eye.x-gravityCenter.x; diff.y=eye.y-gravityCenter.y; diff.z=eye.z-gravityCenter.z;
		bezierSurface[f].Z=diff.x*diff.x+diff.y*diff.y+diff.z*diff.z;
	}
}

// qsortのための比較関数
int processCompareZvalue(const struct normalVectorAndZvalue *a,const struct normalVectorAndZvalue *b)
{
	if(a->Z<b->Z) return 1;
	else if (a->Z>b->Z) return -1;
	else return 0;
}

// ベジエ曲面上の平面の配列を奥行ソート
void processSortBezierSurfacesByZvalue(void)
{
	// qsort関数により奥行値をクイックソート(比較関数はcompareZvalue)
	qsort(bezierSurface,(usubdivi-1)*(vsubdivi-1),sizeof(struct normalVectorAndZvalue),(int(*)(const void*,const void*))processCompareZvalue);
}

// 最近接の制御点のインデックスを取得
struct index2D getTheNearestCtrlPointIndex(struct vector2D input)
{
	int i,j;
	const float threshold=6.25f; // 最近接であると判断する半径の閾値(の自乗)(2.5f*2.5f)
	unsigned int quadrant; // 引数の象限
	float diffx,diffy; // inputと透視投影後の制御点の差
	float diff; // inputと透視投影後の制御点の距離(の自乗)
	float mindist; // 最小のinputと透視投影後の制御点の距離(の自乗)
	struct index2D output; // 最近接の制御点のインデックス

	// mindistとoutputの初期化
	mindist=threshold; output.k=-1; output.l=-1;

	// 引数の象限の取得
	if(input.x>=0.0f)
	{
		if(input.y>=0.0f) quadrant=1;
		else quadrant=2;
	}
	else
	{
		if(input.y<0.0f) quadrant=3;
		else quadrant=4;
	}

	// 最近接の制御点のインデックスの取得(高速化のため各象限で別処理)
	switch(quadrant)
	{
		case 1:
			for(i=0;i<n+1;i++)
			{
				for(j=0;j<m+1;j++)
				{
					if(ctrlPoint[i][j].p>=0 && ctrlPoint[i][j].q>=0)
					{
						diffx=input.x-ctrlPoint[i][j].p; diffy=input.y-ctrlPoint[i][j].q;
						diff=diffx*diffx+diffy*diffy;
						if(diff<=mindist)
						{
							mindist=diff;
							output.k=i; output.l=j;
						}
					}
				}
			}
			break;
		case 2:
			for(i=0;i<n+1;i++)
			{
				for(j=0;j<m+1;j++)
				{
					if(ctrlPoint[i][j].p>=0 && ctrlPoint[i][j].q<0)
					{
						diffx=input.x-ctrlPoint[i][j].p; diffy=input.y-ctrlPoint[i][j].q;
						diff=diffx*diffx+diffy*diffy;
						if(diff<=mindist)
						{
							mindist=diff;
							output.k=i; output.l=j;
						}
					}
				}
			}
			break;
		case 3:
			for(i=0;i<n+1;i++)
			{
				for(j=0;j<m+1;j++)
				{
					if(ctrlPoint[i][j].p<0 && ctrlPoint[i][j].q<0)
					{
						diffx=input.x-ctrlPoint[i][j].p; diffy=input.y-ctrlPoint[i][j].q;
						diff=diffx*diffx+diffy*diffy;
						if(diff<=mindist)
						{
							mindist=diff;
							output.k=i; output.l=j;
						}
					}
				}
			}
			break;
		case 4:
			for(i=0;i<n+1;i++)
			{
				for(j=0;j<m+1;j++)
				{
					if(ctrlPoint[i][j].p<0 && ctrlPoint[i][j].q>=0)
					{
						diffx=input.x-ctrlPoint[i][j].p; diffy=input.y-ctrlPoint[i][j].q;
						diff=diffx*diffx+diffy*diffy;
						if(diff<=mindist)
						{
							mindist=diff;
							output.k=i; output.l=j;
						}
					}
				}
			}
			break;
	}

	return output;
}

// 初期処理
void myInit(void)
{
	glutSetWindow(main_window);
	glClearColor(1.0f, 1.0f, 1.0f, 1.0f);
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
	// OpenGLでは2次元平行投影として透視投影はソフトで計算
	gluOrtho2D(-100.0f, 200.0f, -150.0f, 100.0f);
	// 制御点の座標を格納する配列の作成
	makeCtrlPointsArray();
	// 制御点の初期化
	processInitCtrlPoints();
	// ベジエ曲面上の点の座標を格納する配列の作成
	makeBezierPointsArray();
	// ベジエ曲面上の点のx, y, z座標の取得
	getBezierPointsCoordinate();
	// ベジエ曲面上の平面の法線ベクトルと奥行値を格納する配列の作成
	makeBezierSurfacesArray();
	// ベジエ曲面上の平面の法線ベクトルの取得
	getBezierSurfacesNormalVector();
}

void myDisplay(void)
{
	if(glutGetWindow()!=main_window)glutSetWindow(main_window);
	// 視点と変換行列の取得
	getEyeCoordinateAndTransformationMatrix();
	// 制御点の透視変換後のx (p), y (q)座標の取得
	getCtrlPointsTransformatedCoordinate();
	// ベジエ曲面上の点の透視変換後のx (p), y (q)座標の取得
	getBezierPointsTransformatedCoordinate();
	// ベジエ曲面上の平面の表裏の識別及び奥行値の取得
	getBezierSurfacesZvalue();
	// ベジエ曲面上の平面の配列を奥行ソート
	processSortBezierSurfacesByZvalue();
	// バッファクリア
	glClear(GL_COLOR_BUFFER_BIT);
	// ベジエ曲面の描画
	drawBezierSurfaces();
	// 制御点の描画
	drawCtrlPoints();
	// 画像の表示
	glutSwapBuffers();
}

void processDisplayWhileCtrlPointSelected(void)
{
	// バッファクリア
	glClear(GL_COLOR_BUFFER_BIT);
	// ベジエ曲面の描画
	drawBezierSurfaces();
	// 制御点の描画
	drawCtrlPoints();
	// 画像の表示
	glutSwapBuffers();
}

// アニメーション
void myIdle(void)
{
	clock_t start, finish;

	if(glutGetWindow()!=main_window)glutSetWindow(main_window);

	counter++;
	start=clock(); // 時間計測の開始

	// 視点と変換行列の取得
	getEyeCoordinateAndTransformationMatrix();
	// 制御点の透視変換後のx (p), y (q)座標の取得
	getCtrlPointsTransformatedCoordinate();
	// ベジエ曲面上の点の透視変換後のx (p), y (q)座標の取得
	getBezierPointsTransformatedCoordinate();
	// ベジエ曲面上の平面の表裏の識別及び奥行値の取得
	getBezierSurfacesZvalue();
	// ベジエ曲面上の平面の配列を奥行ソート
	processSortBezierSurfacesByZvalue();
	// バッファクリア
	glClear(GL_COLOR_BUFFER_BIT);
	// ベジエ曲面の描画
	drawBezierSurfaces();
	// 制御点の描画
	drawCtrlPoints();
	// 回転
	if(animationFlag==1) theta+=1.0f*PI/180.0f;
	else if(animationFlag==2) phi+=1.0f*PI/180.0f;
	else if(animationFlag==3)
	{
		theta+=1.0f*PI/180.0f; phi+=1.0f*PI/180.0f;
	}
	// 画像の表示
	glutSwapBuffers();

	finish=clock(); // 時間計測の終了
	timer+=(float)(finish-start)/CLOCKS_PER_SEC;
	fps=(float)counter/timer;
	// glui2に反映
	glui2->sync_live();
}

void myReshape(int w,int h)
{
	glutSetWindow(main_window);
	glutReshapeWindow(576,480);
	glutPostRedisplay();
}

void myKeyboard(unsigned char key,int x,int y)
{
  switch(key)
  {
	case 27: 
	case 'q':
		exit(0);
		break;
  }
  
  glutPostRedisplay();
}

void myMouse(int button,int state,int x,int y)
{
	struct vector2D pointedCoordinate;

	switch(button)
	{
		case GLUT_LEFT_BUTTON:
			if(state==GLUT_DOWN && showCtrlPointsFlag && animationFlag==0)
			{
				pointedCoordinate.x=(float)x*0.520833333f-100.0f; pointedCoordinate.y=100.0f-(float)y*0.520833333f;
				selectedCtrlPointIndex=getTheNearestCtrlPointIndex(pointedCoordinate);
				if(selectedCtrlPointIndex.k!=-1) ctrlPointSelectedFlag=1;
				processDisplayWhileCtrlPointSelected();
			}
			else if(state==GLUT_UP && ctrlPointSelectedFlag)
			{
				ctrlPointSelectedFlag=0;
				// 制御点の移動に伴う逆透視変換(奥行値を固定)
				processInversePerspectiveTransformation();
				// ベジエ曲面上の点のx, y, z座標の取得
				getBezierPointsCoordinate();
				// ベジエ曲面上の平面の法線ベクトルと奥行値を格納する配列の作成
				makeBezierSurfacesArray();
				// ベジエ曲面上の平面の法線ベクトルの取得
				getBezierSurfacesNormalVector();
				glutPostRedisplay();
			}
			break;
	}
}

void myMotion(int x,int y)
{
	int i,j;
	struct vector2D pointedCoordinate;

	if(ctrlPointSelectedFlag)
	{
		pointedCoordinate.x=(float)x*0.520833333f-100.0f; pointedCoordinate.y=100.0f-(float)y*0.520833333f;
		i=selectedCtrlPointIndex.k; j=selectedCtrlPointIndex.l;
		ctrlPoint[i][j].p=pointedCoordinate.x; ctrlPoint[i][j].q=pointedCoordinate.y;
		processDisplayWhileCtrlPointSelected();
	}
}

// gluiコールバックのためのユーザID
#define N_ID					200
#define M_ID					201
#define USUBDIVI_ID				300
#define VSUBDIVI_ID				301
#define OMEGADEG_ID				400
#define RADIUS_ID				401
#define SHOW_CTRLPOINTS_ID		500
#define SHOW_BEZIERWIRES_ID		501
#define SHOW_BEZIERSURFACES_ID	502
#define HIDDENLINEREMOVAL_ID	600
#define SOURSESHAPE_ID			700
#define ROTATION_ID				800
#define ZOOM_ID					801
#define ANIMATION_ID			900

// myGLUIで使われる変数
GLUI_Checkbox *hiddenLineRemovalCheckbox;
GLUI_Translation *rotationTranslation;
GLUI_Translation *zoomTranslation;
float rot[3],zoom[3];

// コールバック関数
void control_cb(int control)
{
	switch(control)
	{
		case N_ID:
		case M_ID:
			// 制御点の座標を格納する配列の開放
			free(ctrlPoint);
			// 制御点の座標を格納する配列の作成
			makeCtrlPointsArray();
			// 制御点の初期化
			processInitCtrlPoints();
			// ベジエ曲面上の点のx, y, z座標の取得
			getBezierPointsCoordinate();
			// ベジエ曲面上の平面の法線ベクトルと奥行値のポインタを格納する配列の作成
			makeBezierSurfacesArray();
			// ベジエ曲面上の平面の法線ベクトルの取得
			getBezierSurfacesNormalVector();
			// main_windowに反映
			if(glutGetWindow()!=main_window) glutSetWindow(main_window);
			glutPostRedisplay();
			break;
		case USUBDIVI_ID:
		case VSUBDIVI_ID:
			// ベジエ曲面上の点の座標を格納する配列の開放
			free(bezierPoint);
			// ベジエ曲面上の平面の法線ベクトルと奥行値を格納する配列の開放
			free(bezierSurface);
			// ベジエ曲面上の点の座標を格納する配列の作成
			makeBezierPointsArray();
			// ベジエ曲面上の点のx, y, z座標の取得
			getBezierPointsCoordinate();
			// ベジエ曲面上の平面の法線ベクトルと奥行値を格納する配列の作成
			makeBezierSurfacesArray();
			// ベジエ曲面上の平面の法線ベクトルの取得
			getBezierSurfacesNormalVector();
			// main_windowに反映
			if(glutGetWindow()!=main_window) glutSetWindow(main_window);
			glutPostRedisplay();
			break;
		case OMEGADEG_ID:
		case RADIUS_ID:
			omegarad=(float)omegadeg*PI/180.0f; // 画角(弧度法)
			f=radius/(float)tan(omegarad/2.0f); // 視距離
			// glui1に反映
			glui1->sync_live();
			// main_windowに反映
			if(glutGetWindow()!=main_window) glutSetWindow(main_window);
			glutPostRedisplay();
			break;
		case SHOW_CTRLPOINTS_ID:
			counter=0; timer=0.0f; fps=0.0f;
			// glui2に反映
			glui2->sync_live();
			// main_windowに反映
			if(glutGetWindow()!=main_window) glutSetWindow(main_window);
			glutPostRedisplay();
			break;
		case SHOW_BEZIERWIRES_ID:
		case SHOW_BEZIERSURFACES_ID:
			counter=0; timer=0.0f; fps=0.0f;
			// glui2に反映
			glui2->sync_live();
			if(showBezierSurfacesFlag || showBezierWiresFlag!=1) hiddenLineRemovalCheckbox->disable();
			else hiddenLineRemovalCheckbox->enable();
			// main_windowに反映
			if(glutGetWindow()!=main_window) glutSetWindow(main_window);
			glutPostRedisplay();
			break;
		case HIDDENLINEREMOVAL_ID:
			counter=0; timer=0.0f; fps=0.0f;
			// glui2に反映
			glui2->sync_live();
			// main_windowに反映
			if(glutGetWindow()!=main_window) glutSetWindow(main_window);
			glutPostRedisplay();
			break;
		case SOURSESHAPE_ID:
			// 制御点の初期化
			processInitCtrlPoints();
			// ベジエ曲面上の点のx, y, z座標の取得
			getBezierPointsCoordinate();
			// ベジエ曲面上の平面の法線ベクトルと奥行値を格納する配列の作成
			makeBezierSurfacesArray();
			// ベジエ曲面上の平面の法線ベクトルの取得
			getBezierSurfacesNormalVector();
			// main_windowに反映
			if(glutGetWindow()!=main_window) glutSetWindow(main_window);
			glutPostRedisplay();
			break;
		case ROTATION_ID:
			//printf("%f %f\n",rot[0],rot[1]);
			if(rot[0]<0.0f) rot[0]+=360.0f; else if(rot[0]>=360.0f) rot[0]-=360.0f;
			if(rot[1]<0.0f) rot[1]+=360.0f; else if(rot[1]>=360.0f) rot[1]-=360.0f;
			theta=rot[0]*PI/180.0f; phi=rot[1]*PI/180.0f;
			// main_windowに反映
			if(glutGetWindow()!=main_window) glutSetWindow(main_window);
			glutPostRedisplay();
			break;
		case ZOOM_ID:
			//printf("%f\n",zoom[0]);
			d=d0+zoom[0];
			// main_windowに反映
			if(glutGetWindow()!=main_window) glutSetWindow(main_window);
			glutPostRedisplay();
			break;
		case ANIMATION_ID:
			if(animationFlag==0)
			{
				counter=0; timer=0.0f; fps=0.0f;
				// glui2に反映
				glui2->sync_live();
				glutIdleFunc(NULL);
				rotationTranslation->enable(); zoomTranslation->enable();
			}
			else
			{
				glutIdleFunc(myIdle);
				rotationTranslation->disable(); zoomTranslation->disable();
			}
			// main_windowに反映
			if(glutGetWindow()!=main_window) glutSetWindow(main_window);
			glutPostRedisplay();
			break;
	}
}

// GLUIコード
void myGlui()
{
	// glui subwindowの作成
	glui1=GLUI_Master.create_glui_subwindow(main_window,GLUI_SUBWINDOW_RIGHT);
	glui2=GLUI_Master.create_glui_subwindow(main_window,GLUI_SUBWINDOW_BOTTOM);

	// glui1
	GLUI_Panel *degreePanel=glui1->add_panel("Degree");
	GLUI_Spinner *nSpinner=glui1->add_spinner_to_panel
		(degreePanel,"u direction",GLUI_SPINNER_INT,&n,N_ID,control_cb);
	nSpinner->set_int_limits(1, 100);
	nSpinner->set_alignment(GLUI_ALIGN_RIGHT);
	GLUI_Spinner *mSpinner=glui1->add_spinner_to_panel
		(degreePanel,"v direction",GLUI_SPINNER_INT,&m,M_ID,control_cb);
	mSpinner->set_int_limits(1, 100);
	mSpinner->set_alignment(GLUI_ALIGN_RIGHT);

	GLUI_Panel *subdivisionPanel=glui1->add_panel("Subdivision");
	GLUI_Spinner *usubdiviSpinner=glui1->add_spinner_to_panel
		(subdivisionPanel,"u direction",GLUI_SPINNER_INT,&usubdivi,USUBDIVI_ID,control_cb);
	usubdiviSpinner->set_int_limits(2, 300);
	usubdiviSpinner->set_alignment(GLUI_ALIGN_RIGHT);
	GLUI_Spinner *vsubdiviSpinner=glui1->add_spinner_to_panel
		(subdivisionPanel,"v direction",GLUI_SPINNER_INT,&vsubdivi,VSUBDIVI_ID,control_cb);
	vsubdiviSpinner->set_int_limits(2, 300);
	vsubdiviSpinner->set_alignment(GLUI_ALIGN_RIGHT);

	GLUI_Panel *viewpointPanel=glui1->add_panel("Viewpoint");
	GLUI_Spinner *omegadegSpinner=glui1->add_spinner_to_panel
		(viewpointPanel,"angle [deg.]",GLUI_SPINNER_INT,&omegadeg,OMEGADEG_ID,control_cb);
	omegadegSpinner->set_int_limits(1, 180);
	omegadegSpinner->set_alignment(GLUI_ALIGN_RIGHT);
	GLUI_Spinner *radiusSpinner=glui1->add_spinner_to_panel
		(viewpointPanel,"radius",GLUI_SPINNER_FLOAT,&radius,RADIUS_ID,control_cb);
	radiusSpinner->set_float_limits(0.0f, 1000.0f);
	radiusSpinner->set_alignment(GLUI_ALIGN_RIGHT);
	GLUI_Spinner *fSpinner=glui1->add_spinner_to_panel
		(viewpointPanel,"distance",GLUI_SPINNER_FLOAT,&f);
	fSpinner->disable();
	fSpinner->set_alignment(GLUI_ALIGN_RIGHT);

	glui1->add_statictext(""); // hspace

	GLUI_Panel *showWhatPanel=glui1->add_panel("Show what");
	GLUI_Checkbox *showCtrlPointsCheckbox=glui1->add_checkbox_to_panel
		(showWhatPanel,"Control points",&showCtrlPointsFlag,SHOW_CTRLPOINTS_ID,control_cb);
	GLUI_Checkbox *showBezierWiresCheckbox=glui1->add_checkbox_to_panel
		(showWhatPanel,"Bezier wires",&showBezierWiresFlag,SHOW_BEZIERWIRES_ID,control_cb);
	GLUI_Checkbox *showBezierSurfacesCheckbox=glui1->add_checkbox_to_panel
		(showWhatPanel,"Bezier surfaces",&showBezierSurfacesFlag,SHOW_BEZIERSURFACES_ID,control_cb);

	GLUI_Panel *hiddenLineRemovalPanel=glui1->add_panel("Hidden line removal");
	hiddenLineRemovalCheckbox=glui1->add_checkbox_to_panel
	(hiddenLineRemovalPanel,"Hidden line removal",&hiddenLineRemovalFlag,HIDDENLINEREMOVAL_ID,control_cb);
	hiddenLineRemovalCheckbox->disable();

	glui1->add_statictext(""); // hspace

	GLUI_Panel *sourseShapePanel=glui1->add_panel("Source shape");
	GLUI_RadioGroup *showWhatRadiobutton=glui1->add_radiogroup_to_panel
		(sourseShapePanel,&sourceShapeFlag,SOURSESHAPE_ID,control_cb);
	glui1->add_radiobutton_to_group(showWhatRadiobutton,"Convex surface");
	glui1->add_radiobutton_to_group(showWhatRadiobutton,"Sine surface");
	glui1->add_radiobutton_to_group(showWhatRadiobutton,"Plane");

	// glui2
	rotationTranslation=glui2->add_translation("Rotation",GLUI_TRANSLATION_XY,rot,ROTATION_ID,control_cb);
	glui2->add_column(false);
	zoomTranslation=glui2->add_translation("Zoom",GLUI_TRANSLATION_Z,zoom,ZOOM_ID,control_cb);
	glui2->add_column(false);
	GLUI_RadioGroup *animationRadiobutton=glui2->add_radiogroup
		(&animationFlag,ANIMATION_ID,control_cb);
	glui2->add_radiobutton_to_group(animationRadiobutton,"Manual");
	glui2->add_radiobutton_to_group(animationRadiobutton,"theta rotation");
	glui2->add_radiobutton_to_group(animationRadiobutton,"phi rotation");
	glui2->add_radiobutton_to_group(animationRadiobutton,"theta+phi rotation");
	glui2->add_column(false);
	glui2->add_statictext("Benchmark [fps]");
	GLUI_EditText *fpsEdittext=glui2->add_edittext
		("",GLUI_EDITTEXT_FLOAT,&fps);
}

main(int argc, char **argv)
{
	glutInit(&argc,argv);
	
	glutInitDisplayMode(GLUT_DOUBLE|GLUT_RGB);
	glutInitWindowSize(576,480);
	glutInitWindowPosition(100,100);
	main_window=glutCreateWindow("Bezier Surface");

	glutDisplayFunc(myDisplay);
	glutIdleFunc(NULL);
	glutReshapeFunc(myReshape);
	glutKeyboardFunc(myKeyboard);
	glutMouseFunc(myMouse);
	glutMotionFunc(myMotion);

	myGlui();
	myInit();

	glutMainLoop();

	return 0;
}

Ryoichi Mizuno - Scientific Computer Graphics
Supplementary Information: Ryoichi Mizuno - Google+ g+External link