Global Illumination_Screen-Space Ray Tracing(SSR)

本次介绍一种把原本应该在三维空间的Ray Marching 转换到平面空间中采用DDA画线算法的方式来进行纹理采样的屏幕空间局部光追的算法，通过此种转换可以使得在屏幕空间的采样非常均匀，并且减少普通SSR重复采样的问题。

一、DDA画线算法

DDA算法是计算机图形学中最简单的绘制直线算法。其主要思想是由直线公式y = kx + b推导出来的。

本文为了均匀采样，将三维空间的取点转换到屏幕空间用。具体算法原理可以参照DDA算法和Bresenham算法。

DDA画线算法是在屏幕空间最简单的一种画线算法（首先要明白在显示器上画线本质是画一系列的点）。如上图所示，如果需要画一条A到B的线，就需要确定AB线段经过了那些像素，然后把那些像素涂黑即可。

那么DDA画线算法的步骤是怎样的呢？其实非常简单，如上图所示，已知点 A = (x,y) ，那么下一个需要画线的点就是 P0 = A + (△x,△y) ，依次类推再下一个点就是 P1 = P0 + (△x,△y) ，直到画到点 B 。

可以看出DDA画线算法的关键就在于选取 △x,△y ，对于斜率小于1的通常的选取为：

△x = 1 只的是一个像素。对于斜率大于1，为了代码的简洁则可以交换 x,y 。这样对于 x 而言每次像素增加1，对于 y 而言每次最多增加1。DDA画线算法伪代码如下：

if(abs(xB - xA) < abs(yB -yA))
    swap(A,B);
float deltaX = 1.0f;
float deltaY = (yB - yA) / (xB - xA);
Point P = A;
for A to B
    P += float2(deltaX, deltaY);
    DrawPixel(int(P.x),int(P.y));

DDA非常适合GPU，因为GPU的浮点数计算能力更强，而对于CPU中画线，则可以采用浮点计算更少的Bresenham画线算法。

二、屏幕空间的光线反射

2.1 三维空间的Ray Marching

如果想要做出镜面的效果在原理上非常简单，如上图所示，只需要根据入射光线 I 和法线 n 求出反射光线 R 然后计算出与物体的交点 P 即可，那么点 x 的颜色就是 P 点的颜色。那么在具体执行的时候求出交点 P 非常的困难（三角形数量可能巨大，进行射线与大量三角形交点判断代价太高），因此常见的做法就是采用Ray Marching的方式来判断交点，步骤如下（注意以下步骤都是在世界空间中）：

• 对于着色点 x 计算出反射方向 R ;
• 然后以着色点 x 为起点，每次沿着 R 方向行进一个定长得到一个新的点 xi = x + i△p* ;
• 将这个新的点投影到屏幕空间中，得到其屏幕空间中的 uv 坐标。
• 在有了 uv 坐标之后，再从深度缓冲中取得 xi 点所对应的深度 d （由于上述步骤是在屏幕空间，因此还需要将其变换到世界空间），然后比较点 xi 的 z 值与深度 d 即可。

但上面的步骤会带来一个问题，如下图所示

从上图可以看到，非常多的点 xi 对应的是同一个像素，这就导致了一些区域过采样。

又如上图所示，点 xi 之间所采样的深度 d 可能有比较大的间隔，这就导致了一些区域低采样。

2.2 屏幕空间的Ray Marching

为了解决采样不均匀问题，可以把处理全部都转化到屏幕空间中，即对于线段起始点 Q0 和结束点 Q1 （ Q0，Q1 是世界空间的点），首先将其转换到屏幕空间中得到点 H0，H1 （屏幕空间中的点是二维点只有 x，y 分量）。如此一来对于点 Q0，Q1 就可以借助DDA画线算法的思想来进行采样（原本的画点操作就变成了采样操作）。这样的好处就在于绝不会重复采样，也保证会连续采样。

但是由于在屏幕空间进行Ray Marching丢了 Z 值信息，无法进行深度比较，因此还需要单独保存当在屏幕空间从 H0 变化到 H1 的 Z 值。

有人可能会说，当的时候，同样也让[公式] 以相同的步数变化到 [公式]不就好了吗？因为屏幕空间和世界空间并不呈一个线性关系，所以不能从 H0 变化到 H1 等比成从 Q0 变化到 Q1 。

当 H0 到 H1 是均匀变化的情况下， Q0 到 Q1 不是均匀变化。原因就是透视除法，因此为了避免这种情况，还需要存储以下变量：

上式的顺序不能乱。当从 H0 变化到 H1 的时候，不在变化 Q0，Q1 而是变化 V0，V1 （此时经历透视除法的 V 已经和屏幕空间的点呈线性相关了），同时为了将点 V0，V1 还原到 Q0，Q1 还需要将从 k0 变化到 k0 。通过 V 还原到 Q ：

三、代码实现

3.1 G-Buffer生成

顶点着色器

#version 430 core
layout(location = 0) in vec3 _Position;
layout(location = 1) in vec3 _Normal;
layout(location = 2) in vec2 _TexCoord;
layout(std140, binding = 0) uniform u_Matrices4ProjectionWorld
{
	mat4 u_ProjectionMatrix;
	mat4 u_ViewMatrix;
};

uniform mat4 u_ModelMatrix;
out vec2 v2f_TexCoords;
out vec3 test;
void main()
{
	vec4 FragPosInViewSpace = u_ViewMatrix * u_ModelMatrix * vec4(_Position, 1.0f);
	gl_Position = u_ProjectionMatrix * FragPosInViewSpace;
	v2f_TexCoords = _TexCoord;
	test = FragPosInViewSpace.xyz;
}

片元着色器

#version 430 core

in  vec2 v2f_TexCoords;
in  vec3 test;
layout (location = 0) out vec4 AlbedoAndMetallic_;

uniform mat4 u_TransposeInverseViewModelMatrix;

uniform sampler2D u_DiffuseTexture;
uniform sampler2D u_NormalTexture;
uniform float u_Near = 0.1;
uniform float u_Far = 1000.0f;

void main()
{
	float gamma = 2.2;
	vec3 diffuseColor = pow(texture(u_DiffuseTexture, v2f_TexCoords).rgb, vec3(gamma));	
	vec3 mapped = vec3(1.0) - exp(-diffuseColor);
	mapped = pow(mapped, vec3(1.0f / 2.2f));
	AlbedoAndMetallic_ = vec4(mapped,1.0);
	float alpha = textureLod(u_DiffuseTexture, v2f_TexCoords,0).a;
	if(alpha != 1.0f)
		discard;
}

3.2 SSR计算

顶点着色器

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

layout(std140, binding = 0) uniform u_Matrices4ProjectionWorld
{
	mat4 u_ProjectionMatrix;
	mat4 u_ViewMatrix;
};

uniform mat4 u_ModelMatrix;

out vec3 v2f_FragPosInWorldSpace;
void main()
{
	vec4 FragPosInWorldSpace = u_ModelMatrix * vec4(_Position, 1.0f);
	gl_Position = u_ProjectionMatrix * u_ViewMatrix * FragPosInWorldSpace;
	v2f_FragPosInWorldSpace = vec3(FragPosInWorldSpace);
}

片元着色器

#version 430 core

in vec3 v2f_FragPosInWorldSpace;
layout(location = 0) out vec4 Color_;

uniform vec3 u_DiffuseColor;
uniform sampler2D u_DepthTexture;
uniform sampler2D u_AlbedoTexture;
uniform float u_Near = 0.1;
uniform float u_Far = 1000.0f;
uniform float u_WindowWidth;
uniform float u_WindowHeight;
uniform vec3  u_CameraPosInWorldSpace;
uniform float u_RayLength = 10000;

layout(std140, binding = 0) uniform u_Matrices4ProjectionWorld
{
	mat4 u_ProjectionMatrix;
	mat4 u_ViewMatrix;
};

float LinearizeDepth(float vDepth)
{
    float z = vDepth * 2.0 - 1.0; 
    return (2.0 * u_Near * u_Far) / (u_Far + u_Near - z * (u_Far - u_Near));    
}

struct Ray
{
    vec3 Origin;
    vec3 Direction;
};

struct Result
{
    bool IsHit;

    vec2 UV;
    vec3 Position;

    int IterationCount;
};

vec4 projectToScreenSpace(vec3 vPoint)
{
	return u_ProjectionMatrix * vec4(vPoint,1);
}

vec3 projectToViewSpace(vec3 vPointInViewSpace)
{
	return vec3(u_ViewMatrix * vec4(vPointInViewSpace,1));
}

float distanceSquared(vec2 A, vec2 B) 
{
    A -= B;
    return dot(A, A);
}
bool Query(vec2 z, vec2 uv)
{
    float depths = -LinearizeDepth(texture(u_DepthTexture, uv / vec2(u_WindowWidth,u_WindowHeight)).r);
    return z.y < depths && z.x > depths;
}

Result RayMarching(Ray vRay)
{
	Result result;

	vec3 Begin = vRay.Origin;
	vec3 End = vRay.Origin + vRay.Direction * u_RayLength;

	vec3 V0 = projectToViewSpace(Begin);
	vec3 V1 = projectToViewSpace(End);

	vec4 H0 = projectToScreenSpace(V0);
	vec4 H1 = projectToScreenSpace(V1);

	float k0 = 1.0 / H0.w;
    float k1 = 1.0 / H1.w;

	vec3 Q0 = V0 * k0; 
    vec3 Q1 = V1 * k1;

	// NDC-space
    vec2 P0 = H0.xy * k0;
    vec2 P1 = H1.xy * k1;
	vec2 Size = vec2(u_WindowWidth,u_WindowHeight);
	//Screen Space
	P0 = (P0 + 1) / 2 * Size;
	P1 = (P1 + 1) / 2 * Size;

	P1 += vec2((distanceSquared(P0, P1) < 0.0001) ? 0.01 : 0.0);

	vec2 Delta = P1 - P0;

	bool Permute = false;
    if (abs(Delta.x) < abs(Delta.y)) 
	{ 
        Permute = true;
        Delta = Delta.yx; P0 = P0.yx; P1 = P1.yx; 
    }
	float StepDir = sign(Delta.x);
    float Invdx = StepDir / Delta.x;
	vec3  dQ = (Q1 - Q0) * Invdx;
    float dk = (k1 - k0) * Invdx;
    vec2  dP = vec2(StepDir, Delta.y * Invdx);
	float Stride = 1.0f;

    dP *= Stride; dQ *= Stride; dk *= Stride;

    P0 += dP; Q0 += dQ; k0 += dk;
	
	int Step = 0;
	int MaxStep = 5000;
	float k = k0;
	float EndX = P1.x * StepDir;
	vec3 Q = Q0;
	float prevZMaxEstimate = V0.z;

	for(vec2 P = P0;  Step < MaxStep;Step++,P += dP, Q.z += dQ.z, k += dk)
	{
		result.UV = Permute ? P.yx : P;
		vec2 Depths;
		Depths.x = prevZMaxEstimate;
		Depths.y = (dQ.z * 0.5 + Q.z) / (dk * 0.5 + k);
		prevZMaxEstimate = Depths.y;
		if(Depths.x < Depths.y)
			Depths.xy = Depths.yx;
		if(result.UV.x > u_WindowWidth || result.UV.x < 0 || result.UV.y > u_WindowHeight || result.UV.y < 0)
			break;
		result.IsHit = Query(Depths, result.UV);
		if (result.IsHit)
			break;
	}

	return result;
}

void main()
{
	vec3 OrginPoint = v2f_FragPosInWorldSpace;
	vec3 ViewDir = normalize(u_CameraPosInWorldSpace - OrginPoint);
	vec3 Normal = vec3(0,1,0);
	vec3 ReflectDir = normalize(reflect(-ViewDir,Normal));
	Ray ray;
	ray.Origin = OrginPoint;
	ray.Direction = ReflectDir;
	Result result = RayMarching(ray);
	vec3 End = OrginPoint;
	vec3 V1 = projectToViewSpace(End);
	vec4 H1 = projectToScreenSpace(V1);

    float k1 = 1.0 / H1.w;

    vec3 Q1 = V1 * k1;
	if(result.IsHit)
	{
		Color_ = vec4(texture(u_AlbedoTexture,result.UV / vec2(u_WindowWidth,u_WindowHeight)).xyz, 1);
	}
	else
	{
		vec4 PointInScreen = u_ProjectionMatrix * u_ViewMatrix * vec4(OrginPoint,1);
		PointInScreen.xy /= PointInScreen.w;
		PointInScreen.xy = PointInScreen.xy * 0.5 + 0.5;
		Color_ = vec4(texture(u_AlbedoTexture,PointInScreen.xy).xyz, 1);
	}
}

3.3 显示

顶点着色器

#version 430 core

layout (location = 0) in vec2 _Position;
layout (location = 1) in vec2 _TexCoords;

out vec2 v2f_TexCoords;

void main()
{
	gl_Position = vec4(_Position, 0.0, 1.0);
	v2f_TexCoords = _TexCoords;
}

片元着色器

#version 430 core

in  vec2 v2f_TexCoords;
out vec4 Color_;

uniform sampler2D u_Texture;
void main()
{
	Color_ = vec4(texture(u_Texture, v2f_TexCoords).rgb, 1.0f);
}