體繪製RayCasting（光線投影演算法）C++/OpenGL原始碼

RayCasting（光線投影演算法）C++原始碼

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include <GL/glut.h>
#define EPSILON 0.000001
#define WIDTH 400
#define HEIGTH 500
float Image[WIDTH*HEIGTH*4];

void GenerateVolume(int *Data,int* Dim);//生成體資料
void GenCube(int x,int y,int z,int side,int density, int *Data,int *Dim);//生成正方體資料
void GenSphere(int x,int y,int z,int radius,int density,int *Data,int *Dim);//生成球體資料
void Classify(float* CData,int *Data,int *Dim);//資料分類
void RotationMatrix(float *R,float *eye,float *center,float *up);//求取影象空間到物體空間變換的旋轉矩陣
void Composite(float *rgba,int x0,int y0,float *CData,int *Dim,float *R,float *T);//合成畫素顏色值
bool Intersection(float *startpos, float *pos, float *dir,int *Dim);//求光線與包圍盒交點座標
void TrInterpolation(float *rgba,float *pos,float *CData,int *Dim);//三線性插值
bool CheckinBox(float *point,int *Dim);//判斷點是否在包圍盒內
void MatrixmulVec(float *c,float *a,float *b);//矩陣向量乘積
void CrossProd(float *c,float *a,float *b);//向量叉乘
void Normalize(float *norm,float *a);//向量歸一化
void Mydisplay();//顯示影象

int main(int argc,char **argv)
{
	int Dim[3]={200,200,200};//體資料大小
	int *Data=(int *)malloc(sizeof(int)*Dim[0]*Dim[1]*Dim[2]);
	float *CData=(float*)malloc(sizeof(float)*Dim[0]*Dim[1]*Dim[2]*4);
	float R[9];//旋轉矩陣
	float T[3]={0,0,450};//平移向量（要根據R調整，以保證獲得體資料全貌）
	float eye[3]={0.5,0.5,1};//視點位置
	float center[3]={0,0,0};//物體參考點位置
	float up[3]={0,1,0};//相機朝上的方向
	RotationMatrix(R,eye,center,up);//獲得旋轉矩陣
	GenerateVolume(Data,Dim);//生成原始體資料
	Classify(CData,Data,Dim);//對體資料分類
	float *LinePS=Image;
	for(int j=0;j<HEIGTH;j++)//逐個合成畫素值
	{
		for(int i=0;i<WIDTH;i++)
		{
			Composite(LinePS,i,j,CData,Dim,R,T);
			LinePS+=4;
		}
	}
	free(Data);
	free(CData);
	//使用OpenGL顯示此二維影象
	glutInit(&argc,argv);
	glutInitDisplayMode(GLUT_SINGLE | GLUT_RGBA);
	glutInitWindowSize(500,500);
	glutInitWindowPosition(200,200);
	glutCreateWindow("Ray-Casting");
	glClearColor (1,1,1,1);//背景設為白色
	glutDisplayFunc(Mydisplay);//顯示影象
	glutMainLoop();
}

//生成體資料
//********************************************************************//
//生成一個大正方體，內部包含一個球體，球體中間又包含一個小正方體
//Data:體積資料
//Dim:體資料大小
//********************************************************************//
void GenerateVolume(int *Data,int *Dim)
{
	GenCube(0,0,0,200,100,Data,Dim);//大正方體
	GenSphere(100,100,100,80,200,Data,Dim);//球體
	GenCube(70,70,70,60,300,Data,Dim);//小正方體
}

//生成正方體資料
//********************************************************************//
//x,y,z:正方體左下角座標
//side:正方體邊長
//density:正方體對應的標量值
//Data:體積資料
//Dim:體資料大小
//********************************************************************//
void GenCube(int x,int y,int z,int side,int density, int *Data,int *Dim)
{
	int max_x=x+side,max_y=y+side,max_z=z+side;
	int Dimxy=Dim[0]*Dim[1];
	for(int k=z;k<max_z;k++)
	{
		for(int j=y;j<max_y;j++)
		{
			for(int i=x;i<max_x;i++)
			{
				Data[k*Dimxy+j*Dim[0]+i]=density;
			}
		}
	}
}

//生成球體資料
//********************************************************************//
//x,y,z:球心座標
//radius:球半徑
//density:球體對應的標量值
//Data:體積資料
//Dim:體資料大小
//********************************************************************//
void GenSphere(int x,int y,int z,int radius,int density,int *Data,int *Dim)
{
	int radius2=radius*radius;
	int Dimxy=Dim[0]*Dim[1];
	for(int k=0;k<Dim[2];k++)
	{
		for(int j=0;j<Dim[1];j++)
		{
			for(int i=0;i<Dim[0];i++)
			{
				if((i-x)*(i-x)+(j-y)*(j-y)+(k-z)*(k-z)<=radius2)
				{
					Data[k*Dimxy+j*Dim[0]+i]=density;
				}			
			}
		}
	}
}

//資料分類
//********************************************************************//
//將原始體資料的標量值對映為顏色和不透明度 
//這裡處理的比較簡單，直接將之前生成的資料分三類：大正方體白色、球體紅色、小正方體黃色
//CData:分類後體資料
//Data:原始體資料
//Dim:體資料大小
//********************************************************************//
void Classify(float* CData,int *Data,int *Dim)
{
	int *LinePS=Data;
	float *LinePD=CData;
	for(int k=0;k<Dim[2];k++)
	{
		for(int j=0;j<Dim[1];j++)
		{
			for(int i=0;i<Dim[0];i++)
			{
				if(LinePS[0]<=100)
				{
					//白色
					LinePD[0]=1.0;
					LinePD[1]=1.0;
					LinePD[2]=1.0;
					LinePD[3]=0.005;
				}
				else if(LinePS[0]<=200)
				{
					//紅色
					LinePD[0]=1.0;
					LinePD[1]=0.0;
					LinePD[2]=0.0;
					LinePD[3]=0.015;
				}
				else
				{
					//黃色
					LinePD[0]=1.0;
					LinePD[1]=1.0;
					LinePD[2]=0.0;
					LinePD[3]=0.02;
				}
				LinePS++;
				LinePD+=4;
			}
		}
	}
}

//求取從影象空間到物體空間變換的旋轉矩陣
//********************************************************************//
//功能類似於OpenGL中的gluLookAt函式
//參考：http://blog.csdn.net/popy007/article/details/5120158
//R:旋轉矩陣
//eye:視點位置
//center:物體參考點位置
//up:相機朝上的方向
//********************************************************************//
void RotationMatrix(float *R,float *eye,float *center,float *up)
{
	float XX[3],YY[3],ZZ[3];//影象空間的基向量
	ZZ[0]=eye[0]-center[0];
	ZZ[1]=eye[1]-center[1];
	ZZ[2]=eye[2]-center[2];
	CrossProd(XX,up,ZZ);
	CrossProd(YY,ZZ,XX);
	Normalize(XX,XX);
	Normalize(YY,YY);
	Normalize(ZZ,ZZ);
	//由影象空間基向量構成旋轉矩陣
	R[0]=XX[0];R[1]=YY[0];R[2]=ZZ[0];
	R[3]=XX[1];R[4]=YY[1];R[5]=ZZ[1];
	R[6]=XX[2];R[7]=YY[2];R[8]=ZZ[2];
}

//合成畫素值
//********************************************************************//
//rgba:合成顏色值
//x0,y0:二維影象畫素座標
//CData:分類後體資料
//Dim:體資料大小
//R:旋轉矩陣（換用於影象空間到物體空間的轉換）
//T:平移向量（同上）
//********************************************************************//
void Composite(float *rgba,int x0,int y0,float *CData,int *Dim,float *R,float *T)
{ 
	int stepsize=1;//取樣步長
	float cumcolor[4];//累計顏色值
	cumcolor[0]=cumcolor[1]=cumcolor[2]=cumcolor[3]=0.0;
	float pos[3],dir[3];//投射光線起點、方向
	float startpos[3];//光線與包圍盒近視點處的交點座標
	float samplepos[3];//取樣點座標
	float samplecolor[4];//取樣點顏色
	//採用平行投影，故在影象空間中投射光線的方向(0,0,-1),起點(x0,y0,0)
	pos[0]=x0;pos[1]=y0;pos[2]=0;
	//將光線描述轉換到物體空間
	//*********************************//
	dir[0]=-R[2];dir[1]=-R[5];dir[2]=-R[8];//光線方向在物體空間的表達
	MatrixmulVec(pos,R,pos);//旋轉
	pos[0]+=T[0];//平移
	pos[1]+=T[1];
	pos[2]+=T[2];
	//*********************************//
	if(Intersection(startpos,pos,dir,Dim))//判斷光線與包圍盒是否相交
	{
		samplepos[0]=startpos[0];
		samplepos[1]=startpos[1];
		samplepos[2]=startpos[2];
		while(CheckinBox(samplepos,Dim)&&cumcolor[3]<1)//當光線射出包圍盒或累計不透明度超過1時中止合成
		{
			TrInterpolation(samplecolor,samplepos,CData,Dim);//三線性插值獲得取樣點處的顏色及不透明度
			//合成顏色及不透明度,採用的是從前到後的合成公式
			cumcolor[0] +=samplecolor[0]*samplecolor[3]*(1-cumcolor[3]);//R
			cumcolor[1] +=samplecolor[1]*samplecolor[3]*(1-cumcolor[3]);//G
			cumcolor[2] +=samplecolor[2]*samplecolor[3]*(1-cumcolor[3]);//B
			cumcolor[3] +=samplecolor[3]*(1-cumcolor[3]);				//A
			//下一個取樣點
			samplepos[0]+=dir[0]*stepsize;
			samplepos[1]+=dir[1]*stepsize;
			samplepos[2]+=dir[2]*stepsize;
		}
		rgba[0]=cumcolor[0];
		rgba[1]=cumcolor[1];
		rgba[2]=cumcolor[2];
		rgba[3]=cumcolor[3];
		return;
	}
	rgba[0]=rgba[1]=rgba[2]=rgba[3]=1.0;//若光線與包圍盒不相交，賦白色
}

//判斷投射光線與包圍盒是否相交（若相交，求靠近視點處的交點座標）
//********************************************************************//
//思路：將包圍盒6個面無限擴充套件，並分成3組，即平行於XOY,YOZ,ZOX平面的各有2個；
//求光線與6個平面的交點，從每組的2個交點中選出距離視點較近者，這樣得到3個候
//選交點；從這3個候選交點中選出距離視點最遠的那個。最後判斷這個點是否落在包
//圍盒內，若在，即為光線與包圍盒的靠近視點處的交點。
//stratpos:靠近視點處的交點座標
//pos:光線起點座標
//dir:光線方向向量
//Dim:包圍盒右上角座標（左下角座標為（0,0,0））
//********************************************************************//
bool Intersection(float *startpos, float *pos, float *dir,int *Dim)
{
	float nearscale = -1000000;
	float scale1, scale2;
	//光線與包圍盒平行於YOZ的2個平面交點
	if ((dir[0] <=-EPSILON)||(dir[0] >=EPSILON))
	{
		scale1 = (0- pos[0]) / dir[0];
		scale2 = (Dim[0] -1 - pos[0]) / dir[0];
		//選出靠近視點的交點，並與當前候選點比較，保留較遠者
		if (scale1 < scale2) 
		{
			if (scale1 > nearscale) 
				nearscale = scale1;
		}
		else
		{
			if (scale2 > nearscale)
				nearscale = scale2;
		}
	}
	//光線與包圍盒平行於ZOX的2個平面交點
	if ((dir[1] <=-EPSILON)||(dir[1] >=EPSILON))
	{
		scale1 = (0 - pos[1]) / dir[1];
		scale2 = (Dim[1] -1 - pos[1]) / dir[1];
		//選出靠近視點的交點，並與當前候選點比較，保留較遠者
		if (scale1 < scale2) 
		{
			if (scale1 > nearscale) 
				nearscale = scale1;
		}
		else
		{
			if (scale2 > nearscale) 
				nearscale = scale2;
		}
	}
	//光線與包圍盒平行於XOY的2個平面交點
	if ((dir[2] <=-EPSILON)||(dir[2] >=EPSILON))
	{
		scale1 = (0 - pos[2]) / dir[2];
		scale2 = (Dim[2] -1 - pos[2]) / dir[2];
		//選出靠近視點的交點，並與當前候選點比較，保留較遠者
		if (scale1 < scale2) 
		{
			if (scale1 > nearscale) 
				nearscale = scale1;
		}
		else
		{
			if (scale2 > nearscale) 
				nearscale = scale2;
		}
	}
	startpos[0] = pos[0] + nearscale * dir[0] ;
	startpos[1] = pos[1] + nearscale * dir[1] ;
	startpos[2] = pos[2] + nearscale * dir[2] ;
	return CheckinBox(startpos,Dim);  //判斷該點是否在包圍盒內
}

//三線性插值
//********************************************************************//
//rgba:插值結果
//pos:取樣點座標
//CData:分類後體資料
//Dim:體資料大小
//********************************************************************//
void TrInterpolation(float *rgba ,float *pos,float *CData,int *Dim)
{
	int x0,y0,z0,x1,y1,z1;
	float fx,fy,fz;
	float v0,v1,v2,v3,v4,v5,v6;
	int Slicesize=Dim[0]*Dim[1]*4;
	int Stepsize=Dim[0]*4;
	x0=(int)pos[0];//整數部分
	y0=(int)pos[1];
	z0=(int)pos[2];
	fx=pos[0]-x0;//小數部分
	fy=pos[1]-y0;
	fz=pos[2]-z0;
	x1=x0+1;
	y1=y0+1;
	z1=z0+1;
	if(x1>=Dim[0])x1=Dim[0]-1;//防止越界
	if(y1>=Dim[1])y1=Dim[1]-1;
	if(z1>=Dim[2])z1=Dim[2]-1;
	for(int i=0;i<4;i++)
	{
		//取樣點處的值由鄰近的8個點插值獲得
		v0=CData[z0*Slicesize+y0*Stepsize+4*x0+i]*(1-fx)+CData[z0*Slicesize+y0*Stepsize+4*x1+i]*fx;
		v1=CData[z0*Slicesize+y1*Stepsize+4*x0+i]*(1-fx)+CData[z0*Slicesize+y1*Stepsize+4*x1+i]*fx;
		v2=CData[z1*Slicesize+y0*Stepsize+4*x0+i]*(1-fx)+CData[z1*Slicesize+y0*Stepsize+4*x1+i]*fx;	
		v3=CData[z1*Slicesize+y1*Stepsize+4*x0+i]*(1-fx)+CData[z1*Slicesize+y1*Stepsize+4*x1+i]*fx;
		v4=v0*(1-fy)+v1*fy;
		v5=v2*(1-fy)+v3*fy;
		v6=v4*(1-fz)+v5*fz;
		if(v6>1)v6=1;//防止越界
		rgba[i]=v6;
	}
}

//判斷點是否在包圍盒內
//********************************************************************//
//point:點座標
//Dim:包圍盒右上角座標（左下角座標為（0,0,0））
//********************************************************************//
bool CheckinBox(float *point,int *Dim)
{
	if (point[0] < 0||point[0] >= Dim[0]||point[1] < 0||point[1] >= Dim[1]||point[2] < 0||point[2] >= Dim[2]) 
		return false;
	else
		return true;
}

//矩陣與向量乘積
//********************************************************************//
//c=a*b
//c:輸出向量
//a:輸入矩陣
//b:輸入向量
//********************************************************************//
void MatrixmulVec(float *c,float *a,float *b)
{
	float x,y,z;
	x=a[0]*b[0]+a[1]*b[1]+a[2]*b[2];
	y=a[3]*b[0]+a[4]*b[1]+a[5]*b[2];
	z=a[6]*b[0]+a[7]*b[1]+a[8]*b[2];
	c[0]=x;
	c[1]=y;
	c[2]=z;
}

//向量叉乘
//********************************************************************//
//c=a x b
//c:輸出向量
//a:輸入向量
//b:輸入向量
//********************************************************************//
void CrossProd(float *c,float *a,float *b)
{
	float x,y,z;
	x=a[1]*b[2]-b[1]*a[2];
	y=a[2]*b[0]-b[2]*a[0];
	z=a[0]*b[1]-b[0]*a[1];
	c[0]=x;
	c[1]=y;
	c[2]=z;
}

//向量歸一化
//********************************************************************//
//norm:歸一化結果
//a:輸入向量
//********************************************************************//
void Normalize(float *norm,float *a)
{
	float len=sqrt(a[0]*a[0]+a[1]*a[1]+a[2]*a[2]);
	norm[0]=a[0]/len;
	norm[1]=a[1]/len;
	norm[2]=a[2]/len;
}

//顯示函式
void Mydisplay()
{
	glClear(GL_COLOR_BUFFER_BIT);
	glDrawPixels(WIDTH,HEIGTH,GL_RGBA,GL_FLOAT,Image);//使用OpenGL的繪圖函式
	glFlush();
}