一、目标
鑒于國内網絡能夠搜到關于SVO的文章較少，閱讀國外Paper十分費心費力，我準備将最近閱讀的VCTRenderer代碼中所學到的東西記錄下來，作為自己的知識儲備。
二、VCTRenderer的大體渲染流程
三、陰影貼圖的計算
VCTRenderer中為了達到軟陰影的效果，沒有使用PCF，而是使用了EVSM——Exponential Variance Shadow Map，該算法結合了Exponential Shadow Map和Variance Shadow Map，能夠在達到軟陰影效果的同時有效減少光滲透現象。
1、Exponential Shadow Map
首先給出Shadow Map的數學模型：
該公式中，d(x)代表的是在Fragment Shader中将x變換到光源空間時所得到的深度值，z(p)代表的是從陰影貼圖中紋理坐标為p的位置采樣得到的深度值。f則是映射關系，比如硬陰影的映射關系如下：
if (d(x) < z(p))
return 1;//被遮擋
else
return 0;//未被遮擋
為了得到軟陰影，我所知道的一般有兩種處理方式，第一種是不對陰影貼圖做任何處理，在Fragment Shader中通過PCF多次采樣進行模糊處理，代碼如下：
shadow = 0.0;
vec2 texelSize = 1.0 / textureSize(ShadowMap, 0);
for(int x = -1; x <= 1; ++x)
{
for(int y = -1; y <= 1; ++y)
{
float pcfDepth = texture(ShadowMap, projCoords.xy + vec2(x, y) * texelSize).r;
shadow += currentDepth - bias > pcfDepth ? 1.0 : 0.0;
}
}
shadow /= 9.0;
第二種方式是直接對生成後的陰影貼圖進行模糊處理，例如高斯模糊，然後在Fragment Shader中進行采樣，但是由于進行了模糊處理——或者說卷積處理：
那麼将f(d(x),z(p))替換g(p-q)得到下式：
因為我們隻需要對z(p)進行卷及處理，為了将d(x)與z(p)分離，Exponential Shadow Map使用了指數函數代替了Shadow Map的數學模型：
将上式替換f(d(x),z)，且由于我們隻關心Shadow Map可以認為d(x)是一個常數（因為它在Fragment Shader中對于同一個像素來說，值恒定），可以得到下式：
由此可見，成功将d(x)與z(p)分離，我們接下來隻需要對z(p)進行卷積操作即可。
PS：該節内容主要參考自該部落格。
2、Variance Shadow Map
Variance Shadow Map的基本思路與Exponential Shadow Map相同，這裡直接貼上一篇講的非常好的文章：方差陰影貼圖與切比雪夫不等式。
3、Exponential Variance Shadow Map
Exponential Variance Shadow Map是Exponential Shadow Map與Variance Shadow Map相結合的産物，網絡上的相關資料比較少（沒有找到），這裡直接貼上VCTRenderer中的實作方式：
depth_texture.frag中：
vec2 WarpDepth(float depth)
{
// rescale depth into [-1, 1]
depth = 2.0 * depth - 1.0;
float pos = exp( exponents.x * depth);//exponents的值為vec2(40,5)
float neg = -exp(-exponents.y * depth);
return vec2(pos, neg);
}
vec4 ShadowDepthToEVSM(float depth)
{
vec2 moment1 = WarpDepth(depth);
vec2 moment2 = moment1 * moment1;
return vec4(moment1, moment2);
}
void main()
{
vec4 diffuseColor = texture(diffuseMap, texCoord);
if (diffuseColor.a <= alphaCutoff) { discard; }
outColor = ShadowDepthToEVSM(gl_FragCoord.z);
}
計算光照的Shader中：
float linstep(float low, float high, float value)
{
return clamp((value - low) / (high - low), 0.0f, 1.0f);
}
float ReduceLightBleeding(float pMax, float Amount)
{
return linstep(Amount, 1, pMax);
}
vec2 WarpDepth(float depth)
{
depth = 2.0f * depth - 1.0f;
float pos = exp(exponents.x * depth);
float neg = -exp(-exponents.y * depth);
return vec2(pos, neg);
}
float Chebyshev(vec2 moments, float mean, float minVariance)
{
if(mean <= moments.x)
{
return 1.0f;
}
else
{
float variance = moments.y - (moments.x * moments.x);
variance = max(variance, minVariance);
float d = mean - moments.x;
float lit = variance / (variance + (d * d));
return ReduceLightBleeding(lit, lightBleedingReduction);
}
}
// 指數方差陰影貼圖計算可見度
float Visibility(vec3 position)
{
//變換到光源視角空間
vec4 lsPos = lightViewProjection * vec4(position, 1.0f);
// 避免除以0造成算數錯誤
if(lsPos.w == 0.0f)
return 1.0f;
// 轉換到NDC空間
lsPos /= lsPos.w;
vec4 moments = texture(shadowMap, lsPos.xy);
// 偏移Z值避免發生light acen現象
vec2 wDepth = WarpDepth(lsPos.z - 0.0001f);
// 計算一個可接受的最小方差
vec2 depthScale = 0.0001f * exponents * wDepth;
vec2 minVariance = depthScale * depthScale;
// 分别對正數項和負數項使用切比雪夫不等式計算可見度
float positive = Chebyshev(moments.xz, wDepth.x, minVariance.x);
float negative = Chebyshev(moments.yw, wDepth.y, minVariance.y);
// 選取最小的可見度
return min(positive, negative);
}
四、模型體素化
1、什麼是體素化
如上圖所示，用一個個小立方體将模型表示出來即為體素化過程，比如《我的世界》這款遊戲就可以當做是用體素将一個個模型構造出來。
2、使用正交投影矩陣，以及判斷投影方向
因為體素本質上是一個個軸對齊的立方體，為了将模型的每個三角面用一堆體素表示出來，我們需要使用正交投影，如下圖：
但我們都清楚的是，一個立方體有六個面，難道我們需要進行六次投影嗎？答案是否定的，首先由于是正交投影，是以投影到上面和下面的結果是一樣的，左右面、前後面同理，是以可以首先排除三個面，再有就是裂縫的問題，如下兩張圖：
第一幅圖是投影到右面得到的結果，第二幅圖是投影到上面得到的結果，很明顯的看到第二幅圖體素化的結果中，體素與體素之間有明顯的裂縫，為了解決這個問題，我們要讓模型投影後的面積盡可能大，解決方法很簡單，隻需要判斷頂點的法線的最大分量即可，代碼如下：
int CalculateAxis()
{
vec3 p1 = gl_in[1].gl_Position.xyz - gl_in[0].gl_Position.xyz;
vec3 p2 = gl_in[2].gl_Position.xyz - gl_in[0].gl_Position.xyz;
vec3 faceNormal = cross(p1, p2);
float nDX = abs(faceNormal.x);
float nDY = abs(faceNormal.y);
float nDZ = abs(faceNormal.z);
if( nDX > nDY && nDX > nDZ )
{
return 0;//投影到X軸方向上
}
else if( nDY > nDX && nDY > nDZ )
{
return 1;//投影到Y軸方向上
}
else
{
return 2;//投影到Z軸方向上
}
}
3、解決孔洞問題
選擇了投影方向後，如果我們直接進行投影操作會出現孔洞問題，其原因是我們最終要在Fragment Shader中完成體素化過程，而在執行Fragment Shader之前會有一個光栅化過程，該過程中會判斷模型是否覆寫了像素的中心位置，如果覆寫了則表示該像素需要執行後續的Fragment Shader，否則則不需要，是以模型的一些部分可能會被丢失，為了解決這個問題，需要使用一種稱之為保守光栅化的算法。
如上圖所示，第一步我們需要給三角形面片構造一個包圍盒，由于這是投影後的三角形面片，是以隻需要一個二維AABB包圍盒即可：
vec4 AxisAlignedBoundingBox(vec4 pos[3], vec2 pixelDiagonal)
{
vec4 aabb;
aabb.xy = min(pos[2].xy, min(pos[1].xy, pos[0].xy));
aabb.zw = max(pos[2].xy, max(pos[1].xy, pos[0].xy));
// 讓包圍盒稍微大一點
aabb.xy -= pixelDiagonal;
aabb.zw += pixelDiagonal;
return aabb;
}
第二步，計算三角面片投影後所在的平面方程的系數
vec4 trianglePlane;
trianglePlane.xyz = cross(pos[1].xyz - pos[0].xyz, pos[2].xyz - pos[0].xyz);
trianglePlane.xyz = normalize(trianglePlane.xyz); //計算平面的法線
trianglePlane.w = -dot(pos[0].xyz, trianglePlane.xyz); //計算平面距原點的距離
第三步，計算向外擴充後的三角形面片的頂點坐标
由于是正交投影，我們其實可以先忽略頂點坐标的Z分量，隻關心X、Y分量，并以三角形面片的每條邊的兩個端點以及原點構造出三個平面。
已知平面方程為：
忽略Z分量，并且該平面過原點，可知：
又因為：
化簡後可得：
是以我們首先構造出三個過原點的平面：
vec3 planes[3];
planes[0] = cross(pos[0].xyw - pos[2].xyw, pos[2].xyw);
planes[1] = cross(pos[1].xyw - pos[0].xyw, pos[0].xyw);
planes[2] = cross(pos[2].xyw - pos[1].xyw, pos[1].xyw);
planes[0].z -= dot(halfPixel, abs(planes[0].xy));
planes[1].z -= dot(halfPixel, abs(planes[1].xy));
planes[2].z -= dot(halfPixel, abs(planes[2].xy));
此時planes[i].x即為公式中的a，planes[i].y即為公式中的b，planes[i].z即為公式中的d。
接下來，求出這三個平面相交所得到的交線，并除以W分量得到Xndc與Yndc：
vec3 intersection[3];
intersection[0] = cross(planes[0], planes[1]);
intersection[1] = cross(planes[1], planes[2]);
intersection[2] = cross(planes[2], planes[0]);
intersection[0] /= intersection[0].z;
intersection[1] /= intersection[1].z;
intersection[2] /= intersection[2].z;
最後帶入到第二步求得的三角形面片所在的平面公式中算出Z分量，即為：
其中，intersection[i].x即為公式中的Xndc，intersection[i].y即為公式中的Yndc，trianglePlane[i].x即為公式中的a，trianglePlane[i].y即為公式中的b，trianglePlane[i].z即為公式中的c，trianglePlane[i].w即為公式中的d，Wndc為1。
是以有以下代碼：
float z[3];
z[0] = -(intersection[0].x * trianglePlane.x + intersection[0].y * trianglePlane.y + trianglePlane.w) / trianglePlane.z;
z[1] = -(intersection[1].x * trianglePlane.x + intersection[1].y * trianglePlane.y + trianglePlane.w) / trianglePlane.z;
z[2] = -(intersection[2].x * trianglePlane.x + intersection[2].y * trianglePlane.y + trianglePlane.w) / trianglePlane.z;
pos[0].xyz = vec3(intersection[0].xy, z[0]);
pos[1].xyz = vec3(intersection[1].xy, z[1]);
pos[2].xyz = vec3(intersection[2].xy, z[2]);
4、根據保守光栅化後的頂點位置計算對應的體素坐标
//從ndc空間變換到世界空間
vec4 voxelPos = vec4(pos[i].xyz * pos[i].w, pos[i].w);
voxelPos.xyz = (viewProjectionI * voxelPos).xyz;
//因為imageStore這個3D紋理的大小是volumeDimension^3（比如256*256*256），
//是以要先将世界坐标轉換到模型的AABB包圍盒的坐标系當中
//然後除以voxelSize，獲得該坐标對應的是哪一個體素，voxelSize=voxelScale * volumeDimension
voxelPos.xyz -= worldMinPoint;
voxelPos *= voxelScale;
Out.wsPosition = voxelPos.xyz * volumeDimension;
5、在Fragment Shader中的一些處理
經過上述四個步驟後，已經做到将模型分割成一個個小立方體，由于我們使用3D紋理貼圖來存儲每個體素的資訊，是以需要将每個立方體的資訊進行處理（顔色、法線、金屬度、粗糙度）并通過imageStore函數存儲到3D貼圖當中。
要是用的3D紋理如下（可以添加其他的3D紋理來存儲其他的資訊）：
//volatile保證每次采樣該變量時都要從顯存中重新讀取而不是使用寄存器中的備份
//coherent保證每次采樣資料的時候不會從cache中讀取資料而是從記憶體中讀取，維護cache與記憶體的資料一緻性
//綜上所述volatile coherent保證着色器程式能夠在每次采樣時獲得最新的資料
//題外話memory barrier——記憶體屏障，對于多線程操作，使用該函數可以保證通路記憶體的順序
layout(binding = 0, r32ui) uniform volatile coherent uimage3D voxelAlbedo;
layout(binding = 1, r32ui) uniform volatile coherent uimage3D voxelNormal;
layout(binding = 2, r32ui) uniform volatile coherent uimage3D voxelEmission;
layout(binding = 3, r8) uniform image3D staticVoxelFlag;
首先，由3中的圖檔，我們知道擴充後的三角形面片，有一部分是在AABB包圍盒之外的，對此我們要将其裁減掉：
if( In.position.x < In.triangleAABB.x || In.position.y < In.triangleAABB.y || In.position.x > In.triangleAABB.z || In.position.y > In.triangleAABB.w )
{
discard;
}
接下來是非常重要的一步，因為同一個體素往往對應多個不同的三角形面片，是以我們需要求得平均值：
void imageAtomicRGBA8Avg(layout(r32ui) volatile coherent uimage3D grid, ivec3 coords, vec4 value)
{
value.rgb *= 255.0;
uint newVal = convVec4ToRGBA8(value);
uint prevStoredVal = 0;
uint curStoredVal;
uint numIterations = 0;
//imageAtomicCompSwap：從grid中紋理坐标為coords處采樣，若與prevStoredVal相等則将值替換成newVal，傳回值是未替換之前的紋理坐标所對應的值
//這部分代碼的作用是當發生重疊的時候要計算出一個平均值
while((curStoredVal = imageAtomicCompSwap(grid, coords, prevStoredVal, newVal)) != prevStoredVal && numIterations < 255)
{
prevStoredVal = curStoredVal;
vec4 rval = convRGBA8ToVec4(curStoredVal);
rval.rgb = (rval.rgb * rval.a); // Denormalize
vec4 curValF = rval + value; // Add
curValF.rgb /= curValF.a; // Renormalize
newVal = convVec4ToRGBA8(curValF);
++numIterations;
}
}
以上步驟完成了模型的體素化。
五、通過mipmap來建構體素八叉樹
當我們将模型體素化完畢後，我們需要使用一種資料結構對其進行管理，八叉樹的主要思想是将一個3維空間（立方體）等分成八個次級立方體，并以此不停的疊代下去，直到每個立方體唯一代表一個物體（或其一部分）——葉子結點，是以其應該是一個從根節點進行擴充的樹形結構，但是現在的情況是我們已經完成了體素化——即所有的葉子節點已經被計算出來，那麼通過将相鄰的八個體素合并成一個，并不斷的疊代直到合并成一個和模型的AABB包圍盒同樣大小的體素來建構八叉樹——是不是很像mipmap呢。
1、生成mipmap之前，先計算一下漫反射吧。
VCTRenderer的作者在生成mipmap之前先計算了漫反射和陰影，且分成了直接光和經過一次反射的間接光兩部分的漫反射計算。
第一步計算體素的可見度（即陰影），這裡作者采用了兩種方式，第一種是根據計算得到的陰影貼圖進行計算，其代碼就是三、3：Exponential Variance Shadow Map中所述的代碼。
第二種方式則是使用了光線步進——ray marching的方法來計算體素的可見度：
// 光線追蹤計算陰影
float TraceShadow(vec3 position, vec3 direction, float maxTracingDistance)
{
// 縮放系數
float k = traceShadowHit * traceShadowHit;
// 計算體素在體素空間中的大小
float voxelTexSize = 1.0f / volumeDimension;
// 初始時移動兩倍于體素大小的距離來避免發生自碰撞
float dst = voxelTexSize * 2.0f;
vec3 samplePos = direction * dst + position;
// control variables
float visibility = 0.0f;
// accumulated sample
float traceSample = 0.0f;
while (visibility <= 1.0f && dst <= maxTracingDistance)
{
//判斷邊界條件
if (samplePos.x < 0.0f || samplePos.y < 0.0f || samplePos.z < 0.0f || samplePos.x > 1.0f || samplePos.y > 1.0f || samplePos.z > 1.0f)
{
break;
}
//ceil傳回大于等于X的最小整數值
traceSample = ceil(texture(voxelAlbedo, samplePos).a) * k;
// 如果透明度乘上縮放系數大于1則傳回不可見
if(traceSample > 1.0f - EPSILON)
{
return 0.0f;
}
// 計算可見度
visibility += (1.0f - visibility) * traceSample / dst;
// 按光照方向移動
dst += voxelTexSize;
samplePos = direction * dst + position;
}
return 1.0f - visibility;
}
第二步則是光照的計算，這裡VCTRenderer的作者雖然函數名字寫的是BRDF但是與我所知的BRDF相差甚遠，更像是針對體素的立方體形狀的一種trick算法（個人猜想，BRDF中的D項（法線分布函數）和G項（幾何遮蔽函數）在這裡并沒有很大的作用，因為我們是按照體素進行光照計算，法線分布和幾何遮蔽自然是根據體素的堆疊方式的不同而不同，）：
vec3 BRDF(Light light, vec3 normal, vec3 albedo)
{
float nDotL = 0.0f;
if(normalWeightedLambert == 1)
{
vec3 weight = normal * normal;
// 對X、Y、Z軸三個方向分别計算與光線的夾角
float rDotL = dot(vec3(1.0, 0.0, 0.0), light.direction);
float uDotL = dot(vec3(0.0, 1.0, 0.0), light.direction);
float fDotL = dot(vec3(0.0, 0.0, 1.0), light.direction);
rDotL = normal.x > 0.0 ? max(rDotL, 0.0) : max(-rDotL, 0.0);
uDotL = normal.y > 0.0 ? max(uDotL, 0.0) : max(-uDotL, 0.0);
fDotL = normal.z > 0.0 ? max(fDotL, 0.0) : max(-fDotL, 0.0);
// 按照權重得到夾角大小
nDotL = rDotL * weight.x + uDotL * weight.y + fDotL * weight.z;
}
else
{
nDotL = max(dot(normal, light.direction), 0.0f);
}
return light.diffuse * albedo * nDotL;
}
平行光：
return vec4(BRDF(light, normal, albedo) * visibility, visibility);
點光源：
return vec4(BRDF(light, normal, albedo) * falloff * visibility, visibility);
聚光燈：
return vec4(BRDF(light, normal, albedo) * falloff * spotFalloff * visibility, visibility);
2、間接光的漫反射計算（注意需要先生成Mipmap，我這裡稍微調整了一下順序）
這裡使用了cone trace算法，每個圓錐的示意圖如下：
由上圖可知，圓錐體并不是真正意義上的圓錐體，而是由不同level的3D紋理（或者說體素）并接而成的類錐體的形狀。
每個圓錐體都需要起始點C0、方向Cd、角度
、追蹤的長度t，由這四個變量可以計算出：
由d可以計算出需要在哪個level的mipmap上進行采樣：
Vsize代表的含義是體素在最高的level的mipmap時的大小（這裡一般是128）。
使用單個圓錐體并不能滿足對整個半球空間進行積分運算（光照計算的基礎，每個點都要計算其法線半球内所有光線對該點的影響才能得到最後的光照結果）這個條件，是以需要多個圓錐體進行計算，圓錐體的數量越多、角度
越小結果越趨近于真正的半球積分效果，VCTRenderer的作者這裡使用了四個圓錐體進行光照計算。
vec4 CalculateIndirectLighting(vec3 position, vec3 normal)
{
// 沿着法線前進2/volumeDimension機關距離，因為0級的mipmap并不是256^3而是128^3
position = position + normal * (1.0f / (volumeDimension / 2.0f));
vec4 diffuseTrace = vec4(0.0f);
// 設定上向量——用于計算TBN的輔助向量
const float aperture = 1.0f;
vec3 guide = vec3(0.0f, 1.0f, 0.0f);
if (abs(dot(normal, guide)) == 1.0f)
{
guide = vec3(0.0f, 0.0f, 1.0f);
}
// 計算切線和副切線
vec3 right = normalize(guide - dot(normal, guide) * normal);
vec3 up = cross(right, normal);
for(int i = 0; i < 4; i++)
{
vec3 coneDirection = normal;
coneDirection += propagationDirections[i].x * right + propagationDirections[i].z * up;
coneDirection = normalize(coneDirection); //确定光線追蹤的方向
diffuseTrace += TraceCone(position, coneDirection, aperture) * diffuseConeWeights[i];
}
return clamp(diffuseTrace, 0.0f, 1.0f);
}
vec4 TraceCone(vec3 position, vec3 direction, float aperture)
{
uvec3 visibleFace;
float anisoDimension = volumeDimension / 2.0f;
// 隻對看得見的面的mipmap進行光線追蹤
visibleFace.x = (direction.x < 0.0) ? 0 : 1;
visibleFace.y = (direction.y < 0.0) ? 2 : 3;
visibleFace.z = (direction.z < 0.0) ? 4 : 5;
// 每個軸方向上的權重
vec3 weight = direction * direction;
float voxelSize = 1.0f / anisoDimension;
// 為了避免自身碰撞移動一個體素的距離
float dst = voxelSize;
float diameter = aperture * dst;
vec3 samplePos = position + direction * dst;
float mipLevel = 0.0f;
vec4 coneSample = vec4(0.0f);
vec4 anisoSample = vec4(0.0f);
if(samplePos.x < 0.0f || samplePos.y < 0.0f || samplePos.z < 0.0f || samplePos.x > 1.0f || samplePos.y > 1.0f || samplePos.z > 1.0f)
{
return coneSample;
}
//疊代光線追蹤
while(coneSample.a <= 1.0f && dst <= maxTracingDistanceGlobal)
{
if (checkBoundaries > 0 && (samplePos.x < 0.0f || samplePos.y < 0.0f || samplePos.z < 0.0f || samplePos.x > 1.0f || samplePos.y > 1.0f || samplePos.z > 1.0f))
{
break;
}
// 根據光線前進的距離計算需要采樣的mip等級
mipLevel = log2(diameter * anisoDimension);
mipLevel = max(mipLevel - 1.0f, 0.0f);
// 根據miplevel進行各向異性采樣
anisoSample = weight.x * textureLod(voxelTexMipmap[visibleFace.x], samplePos, mipLevel)
+ weight.y * textureLod(voxelTexMipmap[visibleFace.y], samplePos, mipLevel)
+ weight.z * textureLod(voxelTexMipmap[visibleFace.z], samplePos, mipLevel);
coneSample += (1.0f - coneSample.a) * anisoSample;
// 光線步進
dst += max(diameter, voxelSize);
diameter = dst * aperture;
samplePos = direction * dst + position;
}
return coneSample;
}
3、生成mipmap
VCTRenderer的作者将其分成了兩個步驟，第一步是從原始的3D紋理中生成第0級的mipmap（并沒有把原始的3D紋理作為最終的體素的mipmap紋理）。
第二步則是根據第一步生成的Level0的3Dmipmap紋理生成其餘level的mipmap，shader十分簡單：
#version 430
layout (local_size_x = 8, local_size_y = 8, local_size_z = 8) in;
layout(binding = 0, rgba8) uniform image3D voxelMipmapDst[6];
layout(binding = 5) uniform sampler3D voxelMipmapSrc[6];
uniform vec3 mipDimension;
uniform int mipLevel;
const ivec3 anisoOffsets[] = ivec3[8]
(
ivec3(1, 1, 1),
ivec3(1, 1, 0),
ivec3(1, 0, 1),
ivec3(1, 0, 0),
ivec3(0, 1, 1),
ivec3(0, 1, 0),
ivec3(0, 0, 1),
ivec3(0, 0, 0)
);
void FetchTexels(ivec3 pos, int dir, inout vec4 val[8])
{
for(int i = 0; i < 8; i++)
{
val[i] = texelFetch(voxelMipmapSrc[dir], pos + anisoOffsets[i], mipLevel);
}
}
void main()
{
if(gl_GlobalInvocationID.x >= mipDimension.x || gl_GlobalInvocationID.y >= mipDimension.y || gl_GlobalInvocationID.z >= mipDimension.z)
{
return;
}
ivec3 writePos = ivec3(gl_GlobalInvocationID);
ivec3 sourcePos = writePos * 2;
// fetch values
vec4 values[8];
// x -
FetchTexels(sourcePos, 0, values);
imageStore(voxelMipmapDst[0], writePos,
(
values[0] + values[4] * (1 - values[0].a) +
values[1] + values[5] * (1 - values[1].a) +
values[2] + values[6] * (1 - values[2].a) +
values[3] + values[7] * (1 - values[3].a)) * 0.25f
);
// x +
FetchTexels(sourcePos, 1, values);
imageStore(voxelMipmapDst[1], writePos,
(
values[4] + values[0] * (1 - values[4].a) +
values[5] + values[1] * (1 - values[5].a) +
values[6] + values[2] * (1 - values[6].a) +
values[7] + values[3] * (1 - values[7].a)) * 0.25f
);
// y -
FetchTexels(sourcePos, 2, values);
imageStore(voxelMipmapDst[2], writePos,
(
values[0] + values[2] * (1 - values[0].a) +
values[1] + values[3] * (1 - values[1].a) +
values[5] + values[7] * (1 - values[5].a) +
values[4] + values[6] * (1 - values[4].a)) * 0.25f
);
// y +
FetchTexels(sourcePos, 3, values);
imageStore(voxelMipmapDst[3], writePos,
(
values[2] + values[0] * (1 - values[2].a) +
values[3] + values[1] * (1 - values[3].a) +
values[7] + values[5] * (1 - values[7].a) +
values[6] + values[4] * (1 - values[6].a)) * 0.25f
);
// z -
FetchTexels(sourcePos, 4, values);
imageStore(voxelMipmapDst[4], writePos,
(
values[0] + values[1] * (1 - values[0].a) +
values[2] + values[3] * (1 - values[2].a) +
values[4] + values[5] * (1 - values[4].a) +
values[6] + values[7] * (1 - values[6].a)) * 0.25f
);
// z +
FetchTexels(sourcePos, 5, values);
imageStore(voxelMipmapDst[5], writePos,
(
values[1] + values[0] * (1 - values[1].a) +
values[3] + values[2] * (1 - values[3].a) +
values[5] + values[4] * (1 - values[5].a) +
values[7] + values[6] * (1 - values[7].a)) * 0.25f
);
}
六、全局光照
這裡的全局光照計算其算法與五中計算直接漫反射和間接漫反射的算法差别不大，這裡主要說一下幾個沒叙述過的點。
1、通過深度和紋理坐标計算出世界坐标
因為紋理坐标覆寫的是整個視窗，是以本質上此時紋理坐标等價于頂點轉換到ndc空間後*2-1，是以可由下述函數計算出世界坐标：
// 是以texCoord可以變換到[-1，1]範圍内來當作ndc空間内該片段所對應的點的x、y坐标
// 随後通過采樣深度貼圖獲得深度值，再乘上逆矩陣便可以得到坐标
vec3 PositionFromDepth()
{
float z = texture(gDepth, texCoord).x * 2.0f - 1.0f;
vec4 projected = vec4(texCoord * 2.0f - 1.0f, z, 1.0f);
projected = inverseProjectionView * projected;
return projected.xyz / projected.w;
}
2、射線與AABB求交
設射線：
Origin代表射線的起始點，dir是射線方向。
設平面：
normal是平面的法線，X是平面上一點，d是平面到原點的距離。
可得：
由此求出最大的t并保證交點在AABB上，即可判定射線與AABB的哪個面相交。
bool IntersectRayWithWorldAABB(vec3 ro, vec3 rd, out float enter, out float leave)
{
//因為AABB的法線都是(0,1,0),(1,0,0),(0,0,1)
//是以直接對應分量相除即可
vec3 tempMin = (worldMinPoint - ro) / rd;
vec3 tempMax = (worldMaxPoint - ro) / rd;
vec3 v3Max = max (tempMax, tempMin);
vec3 v3Min = min (tempMax, tempMin);
leave = min (v3Max.x, min (v3Max.y, v3Max.z));
enter = max (max (v3Min.x, 0.0), max (v3Min.y, v3Min.z));
return leave > enter;
}
3、直接光光照計算
vec3 BRDF(Light light, vec3 N, vec3 X, vec3 ka, vec4 ks)
{
vec3 L = light.direction;
vec3 V = normalize(cameraPosition - X);
vec3 H = normalize(V + L);
float dotNL = max(dot(N, L), 0.0f);
float dotNH = max(dot(N, H), 0.0f);
float dotLH = max(dot(L, H), 0.0f);
// 解碼高光功率
float spec = exp2(11.0f * ks.a + 1.0f);
// 菲涅爾效應
vec3 fresnel = ks.rgb + (1.0f - ks.rgb) * pow(1.0f - dotLH, 5.0f);
// 高光因子
float blinnPhong = pow(dotNH, spec);
// 能量守恒
blinnPhong *= spec * 0.0397f + 0.3183f;
// 高光
vec3 specular = ks.rgb * light.specular * blinnPhong * fresnel;
// 漫反射
vec3 diffuse = ka.rgb * light.diffuse;
return (diffuse + specular) * dotNL;
}
平行光：
return BRDF(light, normal, position, albedo, specular) * visibility;
點光源：
return BRDF(light, normal, position, albedo, specular) * falloff * visibility;
聚光燈：
return BRDF(light, normal, position, albedo, specular) * falloff * spotFalloff * visibility;
4、間接光光照計算
vec4 CalculateIndirectLighting(vec3 position, vec3 normal, vec3 albedo, vec4 specular, bool ambientOcclusion)
{
vec4 specularTrace = vec4(0.0f);
vec4 diffuseTrace = vec4(0.0f);
vec3 coneDirection = vec3(0.0f);
if(any(greaterThan(specular.rgb, specularTrace.rgb)))
{
vec3 viewDirection = normalize(cameraPosition - position);
vec3 coneDirection = reflect(-viewDirection, normal);
coneDirection = normalize(coneDirection);
// 鏡面反射錐體角度大小設定，最低數值為1機關梯度，再小的話會嚴重影響性能
float aperture = clamp(tan(HALF_PI * (1.0f - specular.a)), 0.0174533f, PI);
specularTrace = TraceCone(position, normal, coneDirection, aperture, false);
specularTrace.rgb *= specular.rgb;
}
if(any(greaterThan(albedo, diffuseTrace.rgb)))
{
// 這裡的漫反射調整了角度的大小以及圓錐體的數量
const float aperture = 0.57735f;
vec3 guide = vec3(0.0f, 1.0f, 0.0f);
if (abs(dot(normal,guide)) == 1.0f)
{
guide = vec3(0.0f, 0.0f, 1.0f);
}
// Find a tangent and a bitangent
vec3 right = normalize(guide - dot(normal, guide) * normal);
vec3 up = cross(right, normal);
for(int i = 0; i < 6; i++)
{
coneDirection = normal;
coneDirection += diffuseConeDirections[i].x * right + diffuseConeDirections[i].z * up;
coneDirection = normalize(coneDirection);
// cumulative result
diffuseTrace += TraceCone(position, normal, coneDirection, aperture, ambientOcclusion) * diffuseConeWeights[i];
}
diffuseTrace.rgb *= albedo;
}
vec3 result = bounceStrength * (diffuseTrace.rgb + specularTrace.rgb);
return vec4(result, ambientOcclusion ? clamp(1.0f - diffuseTrace.a + aoAlpha, 0.0f, 1.0f) : 1.0f);
}
七、結束
以上就是對VCTRenderer的大體解析，對于光照部分一些細節處理，我并沒有完全看懂，作者是按照Interactive Indirect Illumination Using Voxel Cone Tracing這篇論文寫得cone trace的代碼，如果有疑問可以去看一下這篇論文。