PPClock 8 技术专题之二·改进的水波效果 
早在上个世纪(……)我就看到过一个DOS下C的效果讲解。当时个人水平有限,依样画了个葫芦,结果失败了,完全得不到想要的效果。

这个最早的概念总的来说是一个四点递推的过程(方案A),当前点下一时刻的振幅除了取决于该点当前时刻的振幅以外,还取决于它的上下左右四个点的振幅。我昨天准备照此重搭炉灶的时候本也准备采用这个原理,但是很不幸的,结果跟N年前没有区别,速度奇慢不说,而且效果极其难看。照说流传这么广泛的文章应该没有什么大问题,可是我直到现在仍不得不怀疑它的有效性。

后来我在网上又看到一篇文章,采用的原理是一个八点递推的过程(方案B),而且这个原理源自微分方程,同时也有数学上的表述,看上去比较可靠。可惜已经无法确定是谁写的了,先在这里表示一下敬意。

下面是这篇文章的原文:

这个水波模拟程序和其它常见的算法不太一样。这是用虎克定理来近似模拟水波波动的。具体的算法如下。

我们知道,当水面是平静的是时候。也就是说水面各点处于平衡位置时候。水面上的任何点的受力都是平衡的。但是有任何的扰动发生,就会有点偏离平衡位置。受力也就会不平衡,即会产生加速度。由于这些点受力是遵守虎克定理的,如果知道某一刻水面的状态,理论上,我们可以精确的计算出任何时候的水面状态(外干扰为0)。

当然在我们的程序里,计算是不需要精确的。否则我们得建立复杂的微分方程……。我们在这里做如下的假设:一)水面上的点是离散的。二)一个点在dt的时间里受力是相同的。三)一个点受力只和周围的8个点有关。有了这样的三个假设,我们就可以建立我们的物理近似模型了。



首先,我们来计算一个点的受力和他们周围点的位置的关系。设某个时刻S点的位移为S,速度为Vs

S点的受力 Fs = k * b * [(B1+B2+B3+B4) - 4 * S ] + k * a * [(A1+A2+A3+A4)-4*S]。k为虎克系数。a,b分别是两种距离点上的权系数。(距离不一样,效果当然不一样了。)

所以在这一刻,S 点的加速度就是 a = Fs * dm (dm是水的质量。)。 下个时刻S点的速度就是Vs2 = Vs - a * dt (dt为时间间隔)。下个时刻的惟一为 S2 = S + Vs2 * dt (dt为时间间隔)。

一切都准备好了。现在该是我们实现程序的时候了。水面上的点的位移可以由一个数组来保存起来。我们需要两个数组来保存位移,一个保存当前的位移,一个用来计算下一刻的位移(计算的时候要用到当前的位移)。同时还要有一个数组来保存当前的速度。在程序里,我们计算加速度的时候,可以把dt假设成1。这样可以提高速度。同时可以把dm*k当成一个数就可以了。注意的是,如果这些常数取的不好,就得不到你想要的结果。在我的程序里,dm*k=1。b=1。a=3。dt=1。其次我们还要对水波进行衰减,不然你的水波就停不下来了:-)。水波的扩散具体的代码如下:

void RippleSpread()
{
for(int y=1;y<H-1;y++)
for(int x=1;x<W-1;x++)
{
int a =( 16*pBack[y][x]-
(3*pBack[y-1][x]+
3*pBack[y+1][x]+
3*pBack[y][x-1]+
3*pBack[y][x+1]+
1*pBack[y-1][x-1]+
1*pBack[y-1][x+1]+
1*pBack[y+1][x+1]+
1*pBack[y+1][x-1])
)/16;
pV[y][x] -= a;//计算下一刻的速度
pV[y][x] -= pV[y][x]>>6;//速度的自动衰减
pFront[y][x] =pBack[y][x] + pV[y][x];//下一刻的位移
pFront[y][x]-= pFront[y][x]>>8;//位移的衰减
}
//交换缓冲区
int (*pTemp)[W]=pFront;
pFront=pBack;
pBack=pTemp;
}

最后,我们要根据水面位移的数据来渲染水波。我们知道水面波动后,是有坡度的,我们真正好可以用坡度扭曲图象来产生水波的效果。做法如下。对水面上的任何一个点(X,Y)。算出offsetX和offsetY。然后把图象中(X+offsetX,Y+offsetY)里的数据拷贝到该位置就可以了。这里给出了渲染的代码:

void RenderRipple()
{
int xoff,yoff,nx,ny;
COLOR cl;
DRAWSTRUCT ds;
DRAWSTRUCT dst;
BeginDraw(&dst,lpDDS_Texture);
BeginDraw(&ds,lpDDS_Back);
for(int y=1;y<H-1;y++)
{
for(int x=1;x<W-1;x++)
{
xoff=pFront[y][x+1]-pFront[y][x-1];
yoff=pFront[y-1][x]-pFront[y+1][x];
nx=x+xoff/10;
ny=y+yoff/10;
if(nx>=W) continue;
if(nx<0) continue;
if(ny>=H) continue;
if(ny<0) continue;
GetColor(dst,nx,ny,&cl);
SetColor(ds,x,y,&cl);
}
}
EndDraw(lpDDS_Texture);
EndDraw(lpDDS_Back);
}


上面的这个方法是才正确的,能够生成非常逼真的水波。但是遗憾的一点,因为采用逐点运算,所以运算量非常大,速度仍然跟不上。且因为要递归,所以水波的扩散速度是非常慢的,每个周期才往外扩散一个像素,看上去有点像慢动作……

针对这些不足,我放弃了逐点计算振幅的方法,而采用了一个二维网格(Mesh)来代替相应的像素,计算出每个网格结点的振幅以及折射偏移量,对于网格内的点就使用线性插值。这听起来有点像3D贴图方面的算法,其实确实也差不多。



网格大小的选取是有讲究的,因为水波的扩散速度是每个周期一个网格,所以你可以自行选取网格的大小来取得喜欢的效果,一般取4到16之间。这样一来,计算振幅的开销就会下降到原来的1/16到1/256。

计算网格节点振幅和偏移量跟上面提到的方案B完全一样,顶多修改一下衰减的特性而已。下面我再罗嗦一下插值方面的东西。

已知量是屏幕某一点的相对坐标(X, Y)及其所属网格四个角的折射偏移量,需要从这四个偏移量求出当前点的折射偏移量。线性插值的方法估计大家都能够做得出来,我就不多说了,它实际上是目标缓冲区上的正四边形到源缓冲区内的一个不规则四边形的映射。已知目标四边形中的坐标(X, Y)经过映射后取得源坐标(X', Y'),然后把源缓冲区的像素信息复制到目标缓冲区就可以了。

我在实现这个效果的时候使用了静态的网格大小,目的是使得编译器能够对其进行优化。其实如果代码优化得好的话,在插值的过程中将不会出现乘法运算,只有加法和移位,因而也能轻易取得相当不错的速度。

相关的程序可以在这里下载。源代码我会在发布 PPClock 8 的时候一起放出。最后再哆嗦一点,本程序在底层图形方面选用的是Graphics32库。感觉非常棒,值得推荐!

下面是核心代码,因为博客不支持缩进,所以大家看起来可能会累一点^_^bbb

首先是计算振幅和折射偏移量
procedure TGlobeUI.UpdateWaveMesh;
var
tmp: PWaveMesh;
i, j, a, s: Integer;
MESH_SIZE, MESH_WIDTH, MESH_HEIGHT: Integer;
begin
MESH_SIZE := FMeshSize;
MESH_WIDTH := MAP_WIDTH div MESH_SIZE + 1;
MESH_HEIGHT := MAP_HEIGHT div MESH_SIZE + 1;

case MESH_SIZE of
2: UpdateWaveMesh2;
4: UpdateWaveMesh4;
8: UpdateWaveMesh8;
16: UpdateWaveMesh16;
else
begin
tmp := FWaveBuf;
FWaveBuf := FWaveMesh;
FWaveMesh := tmp;
FHasWave := False;

for j := 1 to MESH_HEIGHT - 2 do
begin
for i := 1 to MESH_WIDTH - 2 do
begin
a := (16 * FWaveBuf[i, j] -
3 * (FWaveBuf[i - 1, j] + FWaveBuf[i + 1, j] + FWaveBuf[i, j - 1] + FWaveBuf[i, j + 1]) -
(FWaveBuf[i - 1, j - 1] + FWaveBuf[i - 1, j + 1] + FWaveBuf[i + 1, j - 1] + FWaveBuf[i + 1, j + 1])
) div 16;

s := FSpeed[i, j] - a;
FSpeed[i, j] := s;

s := FWaveBuf[i, j] + FSpeed[i, j];
FWaveMesh[i, j] := (s * 32 - s) div 32;
end;
end;

for j := 1 to MESH_HEIGHT - 2 do
for i := 1 to MESH_WIDTH - 2 do
begin
FWaveDiffX[i, j] := (FWaveMesh[i + 1, j] - FWaveMesh[i - 1, j]) div 8;
FWaveDiffY[i, j] := (FWaveMesh[i, j + 1] - FWaveMesh[i, j - 1]) div 8;
FHasWave := FHasWave or ((FWaveDiffY[i, j] <> 0) or (FWaveDiffY[i, j] <> 0));
end;
end;
end;
end;

然后是内插并绘制像素
procedure TGlobeUI.DrawWave;
var
i, j, dx, dy, x, y, px, py: Integer;
p1, p2, p3, p4, p_row_start, p_row_end, p: TPoint;
row_start_inc, row_end_inc, p_inc: TPoint;
//m2 : Integer;
MESH_SIZE, MESH_WIDTH, MESH_HEIGHT: Integer;
begin
MESH_SIZE := FMeshSize;
MESH_WIDTH := MAP_WIDTH div MESH_SIZE + 1;
MESH_HEIGHT := MAP_HEIGHT div MESH_SIZE + 1;

case MESH_SIZE of
2: DrawWave2;
4: DrawWave4;
8: DrawWave8;
16: DrawWave16;
else
begin
//m2 := MESH_SIZE * 3;
FTemp.Width := FGlobe.Width;
FTemp.Height := FGlobe.Height;
//FTemp.Assign(FGlobe);

for j := 1 to MESH_HEIGHT - 1 do
begin
for i := 1 to MESH_WIDTH - 1 do
begin
p1 := Point(FWaveDiffX[i - 1, j - 1], FWaveDiffY[i - 1, j - 1]);
p2 := Point(FWaveDiffX[i, j - 1], FWaveDiffY[i, j - 1]);
p3 := Point(FWaveDiffX[i - 1, j], FWaveDiffY[i - 1, j]);
p4 := Point(FWaveDiffX[i, j], FWaveDiffY[i, j]);

//--------------------------------------
p_row_start := Point(p1.X * MESH_SIZE, p1.Y * MESH_SIZE);
row_start_inc := Point(p3.X - p1.X, p3.Y - p1.Y); //scaled by MESH_SIZE times
p_row_end := Point(p2.X * MESH_SIZE, p2.Y * MESH_SIZE);
row_end_inc := Point(p4.X - p2.X, p4.Y - p2.Y); //scaled by MESH_SIZE times

py := (j - 1) * MESH_SIZE;
for y := 0 to MESH_SIZE - 1 do
begin
p := Point(p_row_start.X, p_row_start.Y);
p_inc := Point((p_row_end.X - p_row_start.X) div MESH_SIZE,
(p_row_end.Y - p_row_start.Y) div MESH_SIZE); //also scaled by MESH_SIZE times

px := (i - 1) * MESH_SIZE;
for x := 0 to MESH_SIZE - 1 do
begin
dx := px + p.X div MESH_SIZE;
dy := py + p.Y div MESH_SIZE;

if (dx >= 0) and (dy >= 0) and (dx < MAP_WIDTH) and (dy < MAP_HEIGHT) then
begin
//if (p.X < - m2) or (p.Y < - m2) or (p.X > m2) or (p.Y > m2) then
//FTemp.PixelPtr[px, py]^ := clBlack32 // completely reflected
//else
FTemp.PixelPtr[px, py]^ := FGlobe.PixelPtr[dx, dy]^;
end
else FTemp.PixelPtr[px, py]^ := clBlack32;

Inc(p.X, p_inc.X);
Inc(p.Y, p_inc.Y);
Inc(px);
end;

Inc(p_row_start.X, row_start_inc.X);
Inc(p_row_start.Y, row_start_inc.Y);
Inc(p_row_end.X, row_end_inc.X);
Inc(p_row_end.Y, row_end_inc.Y);
Inc(py);
end;
//--------------------------------------
end;
end;
end;
end;
end;


虽然里面的代码是针对Graphics32的,但并不依赖于Graphics32提供的特性,同样可以移植到TBitmap上。

[ 查看全文 ] ( 34 评论 / 8146 次浏览 )   |  永久链接  |  相关链接  |   ( 3 / 1930 )