Vulkan_基于分形布朗运动的体积云实现

一、分形布朗运动FBM

分形布朗运动（Fractal Brownian Motion）也就是fbm，它不是噪声，但是他可以让噪声有更多的细节。可以看成把不同比例位置的一张噪声合并在一起。

通过添加不同的噪声迭代（octaves），我们以规则的步长（lacunarity）连续增加频率并降低噪声的幅度（增益），我们可以获得更精细的噪声粒度并获得更精细的细节。这种技术称为“分形布朗运动”（FBM），或简称为“分形噪声”。

具体原理可参照thebookofshaders一书详细简介。

二、raymarching 算法

在讲述此算法之前我们首先科普区分下一下相近名词

• Ray tracing： 这其实是个框架，而不是个方法。符合这个框架的都叫ray tracing。这个框架就是从视点发射ray，与物体相交就根据规则反射、折射或吸收。遇到光源或者走太远就停住。一般来说运算量不小。
• Ray casting： 其实这个和volumetric可以脱钩。它就是ray tracing的第一步，发射光线，与物体相交。这个可以做的很快，在Doom 1里用它来做遮挡。
• Path tracing： 是ray tracing +蒙特卡洛法。在相交后会选一个随机方向继续跟踪，并根据BRDF计算颜色。运算量也不小。还有一些小分类，比如Bidirectional path tracing。
• Ray Marching： 顾名思义，是一根ray一步一步向前走（marching），知道与物体相交。基本只用于volumetric，或可以当作volumetric处理的情况。

2.1 算法思想

raymarching 算法思想很直观：首先有一个3D的体纹理，然后从相机发射n条射线，射线有一个采样的步长。当射线处在体纹理中时，每个步长采一次样，获取纹理值（实际上表示该点的密度值），计算光照，然后和该条射线当前累积的颜色值进行混合。

为什么这样就可以渲染出正确的图案呢？因为光路是可逆的，从光源射出的光线经过散射，最终进入摄像机的效果等同于从摄像机发出的射线进行着色和采样，这个raytracing的道理是一样的。

这种算法很适合在GPU上实现，因为每条射线的计算都是独立并行的，GPU在大量并行计算上有先天的优势。为了在GPU上实现，我们需要解决的问题主要有2个：

1. 哪些片段需要raymarching。
2. raymarching的方向和终点在哪里。

2.2 确定raymarching的片段

体绘制首先需要一个载体(proxy geometry)，也就是为了确定屏幕上的哪些像素是属于某个体纹理的。这个问题很容易就让人联想到包围盒，问题也就引刃而解。

我们只需将体纹理的包围盒绘制出来，那么包围盒在屏幕上覆盖的片段自然就是需要进行raymarching的了。如下图所示：

随后只需要执行raymarching的片断着色器即可。

2.3 raymarching的方向和终止点

在使用包围盒作为体绘制的载体时，起/终点就是每根ray进出包围盒时的两个交点。关于如何得到这两个点的坐标，有一种2个pass的算法：

two pass算法：

1. 绘制包围盒的背面，即将OpenGL背面剔除设置为GL_FRONT，并将每个片段的世界坐标保存在纹理缓存中。
2. 绘制正面，将每个片段的坐标和上一个pass中的每个片段的坐标相减，即可的到ray的方向和长度，然后进行raymarching算法，达到长度终止即可（采样时要转换为纹理坐标）。

one pass算法：

1. 获得片段的世界坐标，然后减去视点位置得到ray的方向。然后每次步进时都判断当前的纹理坐标是否超出了包围盒的边界，一旦超出，就停止算法。

本案例中使用one pass算法实现：

三、体积云实现

3.1 CPU端

CPU端用到的是一个128³ 的3D材质提供基本形状和一个32³的3D材质提供细节。128的材质里面用了4个通道：1个Perlin和3个越来越密集的Worely；32的材质里面用了3个通道：3个越来越密集的Worely。这两个三维噪声材质都是在CPU上生成的，因为不需要在实时的情况下对这些噪声进行修改所以没有用GPU每一帧再生成的必要，生成完了之后存储为两个raw file之后每一次运行程序的时候直接尝试读取，读不到就再重新生成。

3.2 着色器实现

着色器我们主要用的是fragment shader，画一个屏幕大小的quad做ray marching。我这里的方法还是先和装云的盒子相交然后在盒子里面ray marching。两个hemisphere的方法应该更好但是我想不清楚坐标的对应关系。在盒子里面放就直接可以转到（0，1）的本地坐标。。

在ray marching的每一步都对两个噪声进行采样然后对不同通道加不同的权重叠出最终效果。

顶点着色器:

#version 430            
layout(location = 0) in vec3 pos_attrib;

out vec2 tex_coord;  

void main(void)
{
	gl_Position = vec4(pos_attrib, 1.0);
	tex_coord = 0.5*pos_attrib.xy + vec2(0.5);
}

片段着色器:

#version 430
layout(location = 0) uniform float time;
layout(location = 1) uniform vec2 windowSize;
layout(location = 2) uniform mat4 _CameraToWorld;
layout(location = 3) uniform mat4 _CameraInverseProjection;
layout(location = 6) uniform vec4 _shape;
layout(location = 7) uniform sampler3D _CloudShape;
layout(location = 8) uniform sampler3D _CloudDetail;
layout(location = 9) uniform vec3 _detail;
layout(location = 10) uniform vec3 _offset;

out vec4 fragcolor;           
in vec2 tex_coord;

vec2 uv;	// 当前像素在屏幕上位置
const int RAY_MARCHING_STEPS = 64;
const float PI = 3.1415926;
const vec4 skyColor = vec4(0.2,0.3,0.6,1);
const vec3 lightColor= vec3(1.0,1.0,1.0);
const vec3 cloudColor= vec3(0.78,0.7,0.65);

const vec3 boxScale= vec3(5,2,4.5);
const vec3 boxPosition= vec3(0,1,0.4);
//射线Info
struct Ray
{
    vec3 origin;
    vec3 direction;
};

Ray CreateRay(vec3 origin, vec3 direction)
{
    Ray ray;
    ray.origin = origin;
    ray.direction = direction;
    return ray;
}

//交点Info
struct RayHit{
	vec3 position;
	float hitDist;
	float alpha;
	float entryPoint;
	float exitPoint;
};

RayHit CreateRayHit()
{
	RayHit hit;
	hit.position = vec3(0,0,0);
	hit.hitDist = 99999;
	hit.alpha = 0;
	hit.entryPoint = 0;
	hit.exitPoint = 99999;
	return hit;
}

//根据相机位置创建射线
Ray CreateCameraRay(vec2 uv)
{
    // 将相机原点转换为世界空间
    vec3 origin = (_CameraToWorld * vec4(0.0f, 0.0f, 0.0f, 1.0f)).xyz;
    // 反转视图空间位置的透视投影
    vec3 direction = (_CameraInverseProjection * vec4(uv, 0.0f, 1.0f)).xyz;
    // 转换方向从相机到世界空间，并归一化
    direction = normalize((_CameraToWorld * vec4(direction, 0.0f)).xyz);
    return CreateRay(origin, direction);
}

//云体包围盒Info
struct Box
{
	vec3 position;
	vec3 scale;
}_box;

//创建云体包围盒
Box CreateBox(vec3 position, vec3 scale)
{
	Box box;
	box.position = position;
	box.scale = scale;
	return box;
}

//包围盒求交并记录
void IntersectBox(Ray ray, inout RayHit bestHit, Box box)
{
	vec3 minBound = box.position - box.scale;
	vec3 maxBound = box.position + box.scale;

	vec3 t0 = (minBound - ray.origin) / ray.direction;
	vec3 t1 = (maxBound - ray.origin) / ray.direction;

	vec3 tsmaller = min(t0, t1);
	vec3 tbigger = max(t0, t1);

	float tmin = max(tsmaller[0], max(tsmaller[1], tsmaller[2]));
    float tmax = min(tbigger[0], min(tbigger[1], tbigger[2]));

	if(tmin > tmax) return;

	// Hit a box!
	if(tmax > 0 && tmin < bestHit.hitDist)
	{
		if(tmin < 0) tmin = 0;
		bestHit.hitDist = tmin;
		bestHit.position = ray.origin + bestHit.hitDist * ray.direction;
		bestHit.alpha = 1;
		// For volumetric rendering 对于体积渲染
		bestHit.entryPoint = tmin;
		bestHit.exitPoint = tmax;
	}
}

float Remap(float v, float l0, float h0, float ln, float hn)
{
	return ln + ((v - l0) * (hn - ln)) / (h0 - l0);
}

float SAT(float v)
{
	if(v > 1) return 1;
	if(v < 0) return 0;
	return v;
}

vec3 GetLocPos(vec3 pos)
{
	// 采样4D纹理
	pos = ((pos -_box.position)/(_box.scale)+vec3(1,1+abs(sin(time*0.1)*sin(time*0.05)*cos(time*0.08)),1))*0.5;
	return pos;
}

// 融合天空颜色
vec4 GetSkyColor(float y)
{
	 y = y*0.5 + 0.5;
	return mix(skyColor,vec4(lightColor/2,1),y);
}

//获取云层密度参数
float SampleDensity(vec3 p)
{	
	// p为世界系下坐标点位置
	vec3 shapePose = min(vec3(0.95),max(vec3(0.05),fract(GetLocPos(p)*_offset.x)));
	vec4 cloudShape =  texture(_CloudShape, shapePose);
	
	// 计算基本形状密度
	float boxBottom = _box.position.y - _box.scale.y;
	float heightPercent = (p.y - boxBottom) / (_box.scale.y);
	float heightGradient = SAT(Remap(heightPercent, 0.0, 0.2,0,1)) * SAT(Remap(heightPercent, 1, 0.7, 0,1));
	float shapeFBM = dot(cloudShape,_shape)*heightGradient;
	
	if(shapeFBM > 0)
	{
		vec3 detailPos = min(vec3(0.95),max(vec3(0.05),fract(GetLocPos(p)*4*_offset.y)));
		vec4 cloudDetail =  texture(_CloudDetail, detailPos);
		// 从基本形状中减去细节噪声(通过逆密度加权，使边缘比中心更容易被侵蚀)
		float detailErodeWeight = pow((1 - shapeFBM), 3);
		float detailFBM = dot(cloudDetail.xyz, _detail)*0.6;
		float density = shapeFBM - (1 - detailFBM)*detailErodeWeight*100;
		if(density < 0) return 0;
		else return density*10;
	}

	return -1;
}

vec3 sunDir = vec3(0,1,0);
const int LIGHT_MARCH_NUM = 8;
//RayMarch算法
float lightMarch(vec3 p)
{
	float totalDensity = 0; 
	vec3 lightDir = normalize(sunDir);
	Ray lightRay = CreateRay(p, lightDir);
	RayHit lightHit = CreateRayHit();
	IntersectBox(lightRay, lightHit,_box);

	// 云体包围盒内的距离
	float distInBox = abs(lightHit.exitPoint - lightHit.entryPoint);

	float stepSize = distInBox / float(LIGHT_MARCH_NUM);
	for(int i = 0; i < LIGHT_MARCH_NUM; i++)
	{
		p += lightDir*stepSize;
		totalDensity += max(0.0, SampleDensity(p)/8*stepSize);
	}
	float transmittance = exp(-totalDensity*2);
	float darknessThreshold = 0.6;
	return darknessThreshold + transmittance*(1-darknessThreshold);
}

//体积云光照融合
void Volume(Ray ray, RayHit hit, inout vec4 result)
{
	vec3 start = hit.position;
	vec3 end = hit.position + ray.direction*abs(hit.exitPoint-hit.entryPoint);
	float len = distance(start,end);
	float stepSize = len / float(RAY_MARCHING_STEPS);

	vec3 eachStep =  stepSize * normalize(end - start);
	vec3 currentPos = start;

	// 云体反射光量
	float lightEnergy = 0;

	//若透过率为1，表示总透过率，即天空的颜色。
	float transmittance = 1;
	for(int i = 0; i < RAY_MARCHING_STEPS; i++)
	{
	float density = SampleDensity(currentPos);
	if(density > 0)
	{
		//对太阳方向的云进行采样(平行光)
		float lightTransmittance = lightMarch(currentPos);
		lightEnergy += density*stepSize*transmittance*lightTransmittance;
		//云层密度越大，透光率越小
		transmittance *= exp(-density*stepSize*0.643);

		if(transmittance < 0.01 || lightEnergy > 2)
			break;

	}
	currentPos += eachStep;	 	
	}

	 vec3 merCloudColor = lightEnergy * lightColor * cloudColor;
	 result.xyz = result.xyz*transmittance + merCloudColor;
}

//光线投射
RayHit CastRay(Ray ray)
{
	RayHit bestHit = CreateRayHit();
	IntersectBox(ray, bestHit, _box);
	return bestHit;
}

void main(void)
{   
	//根据屏幕大小确定屏幕像素位置
	uv = (vec2(gl_FragCoord.xy) + vec2(0.5, 0.5))/windowSize.xy * 2.0 - 1.0;
	// 创建射线
    Ray ray = CreateCameraRay(uv);
	vec4 result = GetSkyColor(-ray.direction.y);
	_box = CreateBox(boxPosition*5, boxScale);
	RayHit hit = CastRay(ray);
	if(hit.alpha != 0)
	{
		Volume(ray,hit,result);
	}
	fragcolor = result;
}

具体实现见着色器注释吧，每步骤都是依托于上述两种算法原理进行实现。你也可以自行调整上述可变参数进行设置效果，上述着色器具体实现如下图：

如果为了实现更多的天空效果也可加上雾化等效果。