水动态模拟与渲染
总结自己学习的流体模拟和渲染概念和算法。
一、流体模拟概述
1.流体模拟方法
水这种流体的模拟有基于经验、物理和数据驱动的三类方法。
（1）基于经验水模拟方法将水面简化为高度场，通过对高度场修改实现水面波浪的效果。目前，对水波面动态波浪的模拟方法主要分为这三类，周期函数法、噪声法、傅里叶频谱法。
A. 周期函数法是在高度场上叠加近似波浪的周期函数，Fournier 和Reeves[1]首次提出使用Gerstner方程构建水波模型。《GPU Gems》书中有具体讲解这里。Gerstner方程:
Q控制波陡度的参数，Q=0，就是通常的正弦波的叠加；A<1，形成尖锐的波峰；A>1，波峰处又交叉。这个方程是正弦函数的变形，用余弦函数将坐标点（x,y）向内挤压，再给出高度的sin值。函数图像如下（作者见水印），这样的波形更符合实际水波形状。
B.噪声法
噪声法是在高度场叠加噪声实现类似随机水面波纹的效果。Johanson和 Lejdfors[2]采用Perlin噪声方法.
（2）基于物理水模拟方法
基于物理的方法是依据流体力学方程The Navier-Stokes Equations (NSE)模拟运动，主要分为两种：基于网格的欧拉方法和基于粒子的拉格朗日方法。第一种方法需生成一系列的网格，然后将N-S方程离散到网格上并研究流体上物理状态（速度、压强、密度等）随时间的变化状态。第二种方法将流体离散成一系列具有速度、密度等物理属性的粒子，通过追踪粒子的运动获得整体的流体运动效果。N-S方程有许多不同的变形，主要模块有下面两个方程组成。
2.流体模拟步骤
每个方法具体步骤都不太一致，这篇主要放上噪声法的实现步骤，下一篇是基于物理SPH方法。
首先，需要对水建模，可以利用一张网格模拟水面，或者立方体网格模拟水体，以及利用粒子建模复杂形态。
（1）LOD（level of detail）网格
网格矩形排列方式是采用靠近摄像头的划分等级多，远离摄像头的等级少，这种方式也符合我们对海洋的观察，靠近视线的海洋细节多，远离视线的海洋细节相对少。下面是我的一些注释。
using System.Collections.Generic;
using UnityEngine;
[RequireComponent(typeof(MeshRenderer), typeof(MeshFilter))]
public class OceanLOD : MonoBehaviour {
// 一片区域网格横纵数量
public int width = 32;
public int height = 32;
// 区域的数量和大小
public int tiles = 2;
public Vector3 size = new Vector3(150f, 1f, 150f);
// 材质
public Material material;
// 组成网格横纵的线条数量
int g_height;
int g_width;
// 网格相关
Vector3[] vertices; //顶点
Vector3[] normals; //法线
Vector4[] tangents; //三角
// Mesh baseMesh;
// LOD，越在靠后List的Mesh，网格越少
public int maxLOD = 4;
List<List<Mesh>> tiles_LOD;
// Use this for initialization
void Start ()
{
// 计算线条数量
g_height = height + 1;
g_width = width + 1;
// LOD，Mesh所在的List的LOD List编号越小，Mesh的网格越多
tiles_LOD = new List<List<Mesh>>();
for(int LOD = 0; LOD < maxLOD; LOD++)
{
tiles_LOD.Add(new List<Mesh>());
}
int index = 0;
for(int y = 0; y < tiles; ++y)
{
for(int x = 0; x < tiles; ++x)
{
Debug.Log("创建了一片水");
float cy = y - Mathf.Floor(tiles * 0.5f);
float cx = x - Mathf.Floor(tiles * 0.5f);
// 创建一片水
GameObject tile = new GameObject("WaterTile");
// 坐标以当前节点为中心
tile.transform.parent = transform;
tile.transform.localPosition = new Vector3(cx * size.x, 0f, cy * size.z);
//tile.transform.localScale = new Vector3(150f, 1f, 150f);
// 添加Mesh渲染组件
tile.AddComponent<MeshFilter>();
tile.AddComponent<MeshRenderer>().material = material;
tile.layer = LayerMask.NameToLayer("Water");
tiles_LOD [index++].Add (tile.GetComponent<MeshFilter> ().mesh);
//tiles_LOD [0].Add (tile.GetComponent<MeshFilter> ().mesh);
}
}
GenerateHeightmap();
// OnDrawGizmos ();
//Mesh tile.GetComponentin<MeshFilter> ().mesh=tiles_LOD;
Debug.Log ("mesh绘制完毕");
}
void Update()
{
//Mesh tile.GetComponent<MeshFilter> ().mesh=meshLOD;
}
// 初始化Mesh信息
void GenerateHeightmap()
{
Mesh mesh = new Mesh();
int y = 0;
int x = 0;
// 创建顶点和uv坐标
Vector3[] vertices = new Vector3[g_height * g_width];
Vector4[] tangents = new Vector4[g_height * g_width];
Vector3[] normals = new Vector3[g_height*g_width];
Vector2[] uv = new Vector2[g_height * g_width];
// uv和顶点坐标的缩放值（如果要创建width*height的网格）
Vector2 uvScale = new Vector2(1.0f / (g_width - 1f), 1.0f / (g_height - 1f));
Vector3 sizeScale = new Vector3(size.x / (g_width - 1f), size.y, size.z / (g_height - 1f));
// 顶点和uv坐标一个一个排列过去，在之前MC创建方块的时候没用使用这样的方法，是每个面就对应四个顶点，很多顶点都重复了
for(y = 0; y < g_height; ++y)
{
for(x = 0; x < g_width; ++x)
{
vertices[y * g_width + x] = Vector3.Scale(new Vector3(x, 0f, y), sizeScale);
uv[y * g_width + x] = Vector2.Scale(new Vector2(x, y), uvScale);
}
}
mesh.vertices = vertices;
mesh.uv = uv;
// 设置切线
for(y = 0; y < g_height; ++y)
{
for(x = 0; x < g_width; ++x)
{
tangents[y * g_width + x] = new Vector4(1f, 0f, 0f, -1f);
}
}
mesh.tangents = tangents;
//baseMesh = mesh;
//set normal
for(y = 0; y < g_height; ++y)
{
for(x = 0; x < g_width; ++x)
{
normals[y * g_width + x] = new Vector3(0f, 1f, 0f);
}
}
mesh.normals = normals;
// 生成LOD对应的网格，数组越靠后，网格越大、数量越少
for (int LOD = 0; LOD < maxLOD; ++LOD)
{
Vector3[] verticesLOD = new Vector3[(int)(height / System.Math.Pow(2, LOD) + 1) * (int)(width / System.Math.Pow(2, LOD) + 1)];
Vector2[] uvLOD = new Vector2[(int)(height / System.Math.Pow(2, LOD) + 1) * (int)(width / System.Math.Pow(2, LOD) + 1)];
Vector3[] normalsLOD=new Vector3[(int)(height / System.Math.Pow(2, LOD) + 1) * (int)(width / System.Math.Pow(2, LOD) + 1)];
int idx = 0;
for(y = 0; y < g_height; y += (int)System.Math.Pow(2,LOD))
{
for(x = 0; x < g_width; x += (int)System.Math.Pow(2, LOD))
{
verticesLOD[idx] = vertices[g_width * y + x];
normalsLOD[idx]=normals[g_width * y + x];
uvLOD[idx++] = uv[g_width * y + x];
}
}
// tiles_LOD中的网格都替换成为LOD优化过的网格
// Debug.Log("tiles_LOD[LOD].Count"+tiles_LOD[LOD].Count);
// for(int k = 0; k < tiles_LOD[LOD].Count; ++k)
// {
// Mesh meshLOD = tiles_LOD[LOD][k];
// meshLOD.vertices = verticesLOD;
// meshLOD.normals = normalsLOD;
// meshLOD.uv = uvLOD;
// }
Mesh meshLOD = tiles_LOD[LOD][0];
meshLOD.vertices = verticesLOD;
meshLOD.normals = normalsLOD;
meshLOD.uv = uvLOD;
}
// 三角顶点信息，一个方块对应两个三角、对应六个顶点
for(int LOD = 0; LOD < maxLOD; ++LOD)
{
int index = 0;
int width_LOD = (int)(width / System.Math.Pow(2, LOD) + 1);
int[] triangles = new int[(int)(height / System.Math.Pow(2, LOD) * width / System.Math.Pow(2, LOD)) * 6];
for(y = 0; y < (int)(height / System.Math.Pow(2, LOD)); ++y)
{
for(x = 0; x < (int)(width / System.Math.Pow(2, LOD)); ++x)
{
// 这边逆时针绘制了，按照以前的测试要顺时针才能看见，可能跟切线法线有关
triangles[index++] = (y * width_LOD) + x;
triangles[index++] = ((y + 1) * width_LOD) + x;
triangles[index++] = (y * width_LOD) + x + 1;
triangles[index++] = ((y + 1) * width_LOD) + x;
triangles[index++] = ((y + 1) * width_LOD) + x + 1;
triangles[index++] = (y * width_LOD) + x + 1;
}
}
// 三角替换
for (int k = 0; k < tiles_LOD[LOD].Count; ++k)
{
Mesh meshLOD = tiles_LOD[LOD][k];
meshLOD.triangles = triangles;
//return meshLOD;
}
}
}（2）shader
将噪声纹理类型选为Normal map，用其不断修改法线方向。冯乐乐著《Unity Shader入门精要》水那一章有这一部分详细实现。
（3）渲染
在得到法线信息后，实现反射和折射以及菲涅尔反射，使用立方体纹理作为环境纹理模拟反射，使用GrabPass获取当前屏幕渲染纹理，使用切线空间下法线方向随像素屏幕坐标进行偏移，再使用该坐标对渲染纹理进行屏幕采样，从而模拟折射效果。
最后噪声法的水面模拟和渲染结果图
参考文献
[1]Fournier A,Reeves WT.A simple model of ocean waves[J].ACM Siggraph Computer Graphics,1986,20(4):75-84.
[2]Claes Johanson,(March 2004),Real-time water rendering - Introducing the projected grid concept,Lund University.