視線シミュレーションにより球の画像を生成するプログラム。出力にlibpngを使用。

gazou_seisei.c

#include <stdio.h>
#include <math.h>
#include <png.h>

#define EPS (1e-7)

/* ---------- 反射の計算関係の関数 ---------- */

/*
 * (s[0],s[1])を(d[0],d[1])に原点を中心とする回転運動で移す行列を求める
 * @param outmat 出力を書き出す行列
 * @param s 移される元の点
 * @param d 移す先の点
 */
void getRotateMatrix(double outmat[2][2],double s[2],double d[2]) {
	double a,b,delta;
	delta=s[0]*s[0]+s[1]*s[1];
	if(delta==0) {
		/* error */
		a=1;
		b=0;
	} else {
		a=(s[0]*d[0]+s[1]*d[1])/delta;
		b=(-s[1]*d[0]+s[0]*d[1])/delta;
	}
	outmat[0][0]=a;outmat[0][1]=-b;
	outmat[1][0]=b;outmat[1][1]=a;
}

/*
 * mat*inを求め、outに格納する
 * @param out 乗算結果の点の出力先
 * @param mat 行列
 * @param in 乗算する点
 */
void mulMatrix(double out[2],double mat[2][2],double in[2]) {
	double outbuf[2];
	outbuf[0]=mat[0][0]*in[0]+mat[0][1]*in[1];
	outbuf[1]=mat[1][0]*in[0]+mat[1][1]*in[1];
	out[0]=outbuf[0];
	out[1]=outbuf[1];
}

/*
 * ベクトルを正規化(スカラー倍して長さを1に)する
 * @param out 結果の出力先
 * @param in 入力
 */
void vectorNormalize(double out[3],double in[3]) {
	double len=sqrt(in[0]*in[0]+in[1]*in[1]+in[2]*in[2]);
	if(len==0) {
		out[0]=0;
		out[1]=0;
		out[2]=0;
	} else {
		out[0]=in[0]/len;
		out[1]=in[1]/len;
		out[2]=in[2]/len;
	}
}

/*
 * inが法線nの平面で鏡面反射したベクトルを求める
 * @param out 出力
 * @param in 入力の反射するベクトル
 * @param n 入力の法線ベクトル
 */
void getVectorMillor(double out[3],double in[3],double n[3]) {
	double mat1[2][2]; /* (x,y,z) -> (x,sqrt(y*y+z*z),0) */
	double mat2[2][2]; /* (x,sqrt(y*y+z*z),0) -> (sqrt(x*x+y*y+z*z),0,0) */
	double imat1[2][2]; /* (x,sqrt(y*y+z*z),0) -> (x,y,z) */
	double imat2[2][2]; /* (sqrt(x*x+y*y+z*z),0,0) -> (x,sqrt(y*y+z*z),0) */
	double n_tvec1_s[2]={n[1],n[2]};
	double n_tvec1_d[2]={sqrt(n[1]*n[1]+n[2]*n[2]),0};
	double n_tvec2_s[2]={n[0],n_tvec1_d[0]};
	double n_tvec2_d[2]={sqrt(n[0]*n[0]+n[1]*n[1]+n[2]*n[2]),0};
	double in_tvec1_s[2]={in[1],in[2]};
	double in_tvec1_d[2];
	double in_tvec2_s[2];
	double in_tvec2_d[2];
	double out_tvec1_s[2];
	double out_tvec1_d[2];
	double out_tvec2_s[2];
	double out_tvec2_d[2];
	/* nがx軸に平行になるように回転 */
	getRotateMatrix(mat1,n_tvec1_s,n_tvec1_d);
	getRotateMatrix(imat1,n_tvec1_d,n_tvec1_s);
	getRotateMatrix(mat2,n_tvec2_s,n_tvec2_d);
	getRotateMatrix(imat2,n_tvec2_d,n_tvec2_s);
	mulMatrix(in_tvec1_d,mat1,in_tvec1_s);
	in_tvec2_s[0]=in[0];in_tvec2_s[1]=in_tvec1_d[0];
	mulMatrix(in_tvec2_d,mat2,in_tvec2_s);
	/* 対象のベクトルのx成分を反転(-1倍) */
	out_tvec2_s[0]=-in_tvec2_d[0];
	out_tvec2_s[1]=in_tvec2_d[1];
	/* 逆回転 */
	mulMatrix(out_tvec2_d,imat2,out_tvec2_s);
	out_tvec1_s[0]=out_tvec2_d[1];
	out_tvec1_s[1]=in_tvec1_d[1];
	mulMatrix(out_tvec1_d,imat1,out_tvec1_s);
	/* 結果を書き出す */
	out[0]=out_tvec2_d[0];
	out[1]=out_tvec1_d[0];
	out[2]=out_tvec1_d[1];
}

/* ---------- 直線と球の交点を求める処理 ---------- */

/*
 * 直線P=ins+t*indと球面の交点を求める
 * @param ot 結果のtを格納する配列
 * @param ins 入力直線の基準となる点
 * @param ind 入力直線の方向ベクトル
 * @param incc 入力球の中心
 * @param incr 入力球の半径
 * @return 交点の個数
 */
int getLineSphereHitPoint(double ot[2],
		double ins[3],double ind[3],double incc[3],double incr) {
	double xt,yt,zt;
	double a,b,c;
	double d;
	double s;
	xt=ins[0]-incc[0];
	yt=ins[1]-incc[1];
	zt=ins[2]-incc[2];
	a=ind[0]*ind[0]+ind[1]*ind[1]+ind[2]*ind[2];
	b=2.0*(ind[0]*xt+ind[1]*yt+ind[2]*zt);
	c=xt*xt+yt*yt+zt*zt-incr*incr;
	d=b*b-4.0*a*c;
	if(d<0)return 0;
	else if(-EPS<d && d<EPS) {
		ot[0]=-b/(2.0*a);
		return 1;
	}
	s=sqrt(d);
	/* a>0なので、ot[0]<ot[1] */
	ot[0]=(-b-s)/(2.0*a);
	ot[1]=(-b+s)/(2.0*a);
	return 2;
}

/* ---------- メイン処理 ---------- */

int mixColor(int x,int y,double xRatio) {
	int r1,g1,b1,r2,g2,b2;
	int r3,g3,b3;
	r1=(x>>16)&0xFF;g1=(x>>8)&0xFF;b1=x&0xFF;
	r2=(y>>16)&0xFF;g2=(y>>8)&0xFF;b2=y&0xFF;
	r3=(int)(r1*xRatio+r2*(1.0-xRatio));
	g3=(int)(g1*xRatio+g2*(1.0-xRatio));
	b3=(int)(b1*xRatio+b2*(1.0-xRatio));
	if(r3<0)r3=0;
	if(r3>255)r3=255;
	if(g3<0)g3=0;
	if(g3>255)g3=255;
	if(b3<0)b3=0;
	if(b3>255)b3=255;
	return (r3<<16)|(g3<<8)|b3;
}

/*
x : 右
y : 奥
z : 上
*/

int getColorOfOnePoint(double sp[3],double d[3]) {
	double sphereCenter[3]={0,50,20};
	double sphereRadius=15;
	int sphereColor=0xFFD700;
	double sphereAlpha=0.4;

	double lightPos[3]={30,0,50};
	int lightColor=0xFFFFFF;
	double lightAlpha=0.4;

	int backGroundColor=0x8080FF;

	int doCalculateLight=0;
	double colPoint[3];
	double colN[3];
	int col;
	double t[2];
	int tnum;
	tnum=getLineSphereHitPoint(t,sp,d,sphereCenter,sphereRadius);
	if(tnum==2 && (t[0]>EPS || t[1]>EPS)) {
		double np[3];
		double n[3],nn[3];
		double nd[3];
		if(t[0]>EPS) {
			np[0]=sp[0]+t[0]*d[0];
			np[1]=sp[1]+t[0]*d[1];
			np[2]=sp[2]+t[0]*d[2];
		} else {
			np[0]=sp[0]+t[1]*d[0];
			np[1]=sp[1]+t[1]*d[1];
			np[2]=sp[2]+t[1]*d[2];
		}
		n[0]=np[0]-sphereCenter[0];
		n[1]=np[1]-sphereCenter[1];
		n[2]=np[2]-sphereCenter[2];
		vectorNormalize(nn,n);
		getVectorMillor(nd,d,nn);
		col=mixColor(sphereColor,getColorOfOnePoint(np,nd),sphereAlpha);
		doCalculateLight=1;
		colPoint[0]=np[0];
		colPoint[1]=np[1];
		colPoint[2]=np[2];
		colN[0]=nn[0];
		colN[1]=nn[1];
		colN[2]=nn[2];
	} else if(sp[2]*d[2]<0) {
		double t=-sp[2]/d[2];
		double x=sp[0]+t*d[0];
		double y=sp[1]+t*d[1];
		int xx=(int)floor(x);
		int yy=(int)floor(y);
		int p=0;
		if(xx<0) {
			xx=-xx;
			if(xx%20<10)p++;
		} else {
			if(xx%20>=10)p++;
		}
		if(yy<0) {
			yy=-yy;
			if(yy%20<10)p++;
		} else {
			if(yy%20>=10)p++;
		}
		col=p%2?0xFFFFFF:0x000000;
		doCalculateLight=1;
		colPoint[0]=x;
		colPoint[1]=y;
		colPoint[2]=sp[2]+t*d[2];
		colN[0]=0;
		colN[1]=0;
		colN[2]=1;
	} else {
		col=backGroundColor;
	}
	if(doCalculateLight) {
		double toLight[3]={
			lightPos[0]-colPoint[0],lightPos[1]-colPoint[1],lightPos[2]-colPoint[2]
		};
		double hwv[3];
		double sot[2];
		double lightValue;
		int sotnum;
		sotnum=getLineSphereHitPoint(sot,colPoint,toLight,sphereCenter,sphereRadius);
		if(sotnum==2 && ((EPS<sot[0] && sot[0]+EPS<1.0) || (EPS<sot[1] && sot[1]+EPS<1.0))) {
			lightValue=0;
		} else {
			vectorNormalize(hwv,d);
			hwv[0]=-hwv[0]+toLight[0];
			hwv[1]=-hwv[1]+toLight[1];
			hwv[2]=-hwv[2]+toLight[2];
			vectorNormalize(hwv,hwv);
			lightValue=hwv[0]*colN[0]+hwv[1]*colN[1]+hwv[2]*colN[2];
			if(lightValue<EPS) {
				lightValue=0;
			} else {
				lightValue/=sqrt(hwv[0]*hwv[0]+hwv[1]*hwv[1]+hwv[2]*hwv[2]);
				lightValue/=sqrt(colN[0]*colN[0]+colN[1]*colN[1]+colN[2]*colN[2]);
				lightValue=pow(lightValue,1.2);
				if(lightValue<0)lightValue=0;
				if(lightValue>1)lightValue=1;
			}
		}
		col=mixColor(mixColor(lightColor,0,lightValue),col,lightAlpha);
	}
	return col;
}

int main(int argc,char* argv[]) {
	double cameraDir[3]={0,1,0};
	double cameraX[3]={1,0,0};
	double cameraY[3]={0,0,1};
	double cameraPos[3]={0,0,25};
	int width=512;
	int height=512;
	int max;
	int x,y;
	char fileName[128];
	FILE* fp=NULL;
	png_structp pngPtr;
	png_infop pngInfoPtr;
	png_bytepp pngImageData;
	if(argc>=2)sscanf(argv[1],"%d",&width);
	if(argc>=3)sscanf(argv[2],"%d",&height);
	if(width<=0 || height<=0) {
		puts("invalid size");
		return 1;
	}
	max=(width>=height?width:height);
	pngPtr=png_create_write_struct(PNG_LIBPNG_VER_STRING,NULL,NULL,NULL);
	if(!pngPtr) {
		puts("png_create_write_struct failed");
		return 1;
	}
	pngInfoPtr=png_create_info_struct(pngPtr);
	if(!pngInfoPtr) {
		puts("png_create_info_struct failed");
		png_destroy_write_struct(&pngPtr,NULL);
		return 1;
	}
	sprintf(fileName,"gazou_seisei_out_%d_%d.png",width,height);
	fp=fopen(fileName,"wb");
	if(setjmp(png_jmpbuf(pngPtr))) {
		png_destroy_write_struct(&pngPtr,&pngInfoPtr);
		fclose(fp);
		return 1;
	}
	png_init_io(pngPtr,fp);
	png_set_IHDR(pngPtr,pngInfoPtr,width,height,8,PNG_COLOR_TYPE_RGB,
		PNG_INTERLACE_NONE,PNG_COMPRESSION_TYPE_BASE,PNG_FILTER_TYPE_BASE);

	puts("generating image data...");
	pngImageData=png_malloc(pngPtr,height*sizeof(png_bytep));
	for(y=0;y<height;y++)pngImageData[y]=png_malloc(pngPtr,width*3);

	for(y=0;y<height;y++) {
		for(x=0;x<width;x++) {
			double xoff=(double)(x-width/2)/(max/2);
			double yoff=(double)(height/2-y)/(max/2);
			double vec[3]={
				cameraDir[0]+cameraX[0]*xoff+cameraY[0]*yoff,
				cameraDir[1]+cameraX[1]*xoff+cameraY[1]*yoff,
				cameraDir[2]+cameraX[2]*xoff+cameraY[2]*yoff
			};
			int color=getColorOfOnePoint(cameraPos,vec);
			pngImageData[y][x*3+0]=(color>>16)&0xFF;
			pngImageData[y][x*3+1]=(color>>8)&0xFF;
			pngImageData[y][x*3+2]=color&0xFF;
		}
	}
	png_set_rows(pngPtr,pngInfoPtr,pngImageData);

	png_write_png(pngPtr,pngInfoPtr,PNG_TRANSFORM_IDENTITY,NULL);
	for(y=0;y<height;y++)png_free(pngPtr,pngImageData[y]);
	png_free(pngPtr,pngImageData);
	png_destroy_write_struct(&pngPtr,&pngInfoPtr);
	fclose(fp);
	puts("done");
	return 0;
}