water Caustics

　　水下环境的模拟在现阶段的计算能力下，Jos Stam(https://www.opengl.org/archives/resources/code/samples/mjktips/caustics/) 提出的方法是预先计算多张caustics图，根据时间选取不同的图片，然后与水下模型进行混合。在underWater.c中，

glEnable(GL_BLEND); 

glBlendFunc(GL_ZERO, GL_SRC_COLOR);

上面这两行代表绘制时与现有模型进行颜色混合。

GLfloat sPlane[4] = { 0.05, 0.03, 0.0, 0.0 };

GLfloat tPlane[4] = { 0.0, 0.03, 0.05, 0.0 };

glTexGeni(GL_S, GL_TEXTURE_GEN_MODE, GL_OBJECT_LINEAR);

glTexGeni(GL_T, GL_TEXTURE_GEN_MODE, GL_OBJECT_LINEAR);

glTexGenfv(GL_S, GL_OBJECT_PLANE, sPlane);

glTexGenfv(GL_T, GL_OBJECT_PLANE, tPlane);

glEnable(GL_TEXTURE_GEN_S);

glEnable(GL_TEXTURE_GEN_T);

glEnable(GL_TEXTURE_2D);

代表着caustics纹理坐标的计算，可以看到sPlane和tPlane的Y值并不为0，首先我们设置Y0时cube的效果，这时候cube的竖直面没有ripple效果，这是因为纹理坐标在Y值没有偏移时，X,Z一致值取得的caustics纹理坐标完全一样。左图是Y值为0的情况，右图是Y值不为0的情况。

　　在GpuGem中，文章用一种简单的构造折射向量的方式来计算纹理的映射，用来映射的纹理就是 sun map，方法分为三步：

1.构造入射向量。

2.根据波当前处理的顶点位置和顶点法线来计算折射向量。

3.根据折射向量来计算映射的纹理坐标。

在计算光照折射效果之前，需要利用波公式来构建波，这里所用的wave function 如下：

下面就是该函数的实现部分：

float func(float x,float z)
{
float y=0;

float factor=1;
float d=sqrt(x*x+z*z);
d=d/40;
if (d>1.5) d=1.5;
if (d<0) d=0;
for (int i=0;i<octaves;i++)
    {
    y-=    (factor)*VTXSIZE*d*cos(((float)timer*SPEED)+(1/factor)*x*z*WAVESIZE)+
            (factor)*VTXSIZE*d*cos(((float)timer*SPEED)+(1/factor)*(x*z)*WAVESIZE)
            ;
    //y-=factor*VTXSIZE*d*noise((float)(x*(1/factor))/5,SPEED*timer,(float)z*(1/factor)/5);
    factor=factor/1.7;
    }

return y;
}

而在波纹理的构造过程中，主要使用下面的公式进行纹理坐标的计算。

double plane::testline(point r,point vec,point &pres)
// return value is distance along the line to the collision point
{
// a*x +b*y+c*z+d =0为平面构造公式，r为波的顶点坐标，vec为顶点的法线，该方法主要
//是求pres，也就是平面与顶点和法线构造的直线相交点。
double fres;
double t;

//fres = |vec|*|n|cos(Q)
fres=vec*n;
pres=r;
if ((fres)!=0)
    {
// r到平面的距离公式= |a*r.x + b*r.y+c*r.z +d|/(sqrt(a*a+b*b+c*c)) = |n*r+d|
// cos(Q) = (nr+d)/t = fres/(|vec|*|n|)
    t=-(n*r + d)/fres;
//pres相交点。
    pres=r+vec*t;
    }
// hack... arreglar
else t=0;
return t;
}

最后计算的是折射光线的效果，原始的snell’s law计算方法不太适用于计算机计算，Foley等在1996提出了下面一个便于计算的公式。

其中n代表面的法线方向，E为入射向量， h _1/ h _{2  是折射率，计算的函数如下。}

void sample_caustic(float xi,float zi,float &u,float &v)
// here‘s the main stuff for the caustic
{
point p(xi,0,zi);p.y=func(xi,zi);
point q(xi+QUADSIZE,0,zi);q.y=func(xi+1,zi);
point r(xi,0,zi+QUADSIZE);r.y=func(xi,zi+1);

point e1=q-p; point e2=r-p;    point n=e1^e2;    // the normal above the sampling point
u=0; v=0;
for (int i=0;i<numrays;i++)
    {
    float alpha=6.28*(float)i/numrays;
    float rad=(float)i/numrays;
    float xf=rad*cos(alpha)/2.0;
    float zf=rad*sin(alpha)/2.0;

    //构造入射向量
    point incident(xf,1,zf);

    float c=incident*n;
    float sq=1+(IOR*IOR*(c*c-1));
    if (sq>0) sq=sqrt(sq); else sq=0;

    //折射向量的计算。
    point transmitted=incident*IOR + n*(IOR*c-sq);
    transmitted.normalize();

    u+=transmitted.x;
    v+=transmitted.z;
    }
u=4*(u/numrays) + 0.5;
v=4*(v/numrays) + 0.5;
}

　　最后比较几张wave function 中octaves参数从1到3的情况。