Skip to content

Instantly share code, notes, and snippets.

@mikecat
Created March 11, 2014 14:33
Show Gist options
  • Save mikecat/9486977 to your computer and use it in GitHub Desktop.
Save mikecat/9486977 to your computer and use it in GitHub Desktop.
視線シミュレーションにより球の画像を生成するプログラム。出力にlibpngを使用。
#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;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment