用Houdini的grid做畫板學習GLSL的fragment shader原理

最近一直在學習使用GLSL，國外有兩個非常好的資源網站 shadertoy 和 glslsandbox。里面有很多關於glsl的案例，網站維護也都是一些圈里的大神在做着，如果有想法學習比較底層的shader可以看看這兩個網站。

這里主要講一講通過在Houdini里面用grid做畫板，通過編寫VEX代碼來實現GLSL中fragment shader的原理。所謂的fragment其實可以直接理解為像素化。從vertex shader出來后所有的頂點數據（位置，顏色等等）都通過fragment shader來進行插值一個像素一個像素點的畫滿整個屏幕。為了理解這一過程，我直接在最熟悉的vex環境里面做了一些練習。其實vex從某種程度上來講也算houdini的shader語言，一方面他本來就是point或者primitive、vertex遍歷的一種語言，除開pcopen帶來的與自身鄰居點產生互動外，其每一個point或primitive都是相互獨立進行運算的，另一方面你要是牛逼點的，mantra的shader也可以直接使用vex開寫。而且vex的處理過程像極了GPU的並行運算方法。在下面的三個練習過程中，也發現了H14與之前的版本相比，在VOP上有一個比較好的改進的改進，這個我留在后面講。

 

練習一 - 在grid上畫一個帶陰影的圓球，其中光源能跟着攝像機走

用fragment shader的思路來闡述這個效果怎么實現，首先我們把grid上面的每一個點看成屏幕上的一個像素點。在shader進行着色的過程中，是並行的對每一個像素點進行着色。也就是說我隨機找grid上的一個點，那么我計算這個點的顏色時，我是假設這個世界上就只有我這一個點了，我的鄰居點在哪我不知道，我也不用知道。就好比我要做一千個餃子，每個餃子都是用一種方法做出來的，就是拿起餃子皮，用筷子戳一塊餡放在餃子皮中間然后卷起餃子皮包起來。GPU牛逼的地方就是你來一千個餃子的訂單他就出一千個工人同時每人包一個餃子然后出貨，一萬個就一萬個工人同時每人包一個，每次一次性能出多少餃子只取決於GPU有多少工人（管線）。雖然上面寫的有點啰嗦，但是這也是最直白解釋並行運算的原理了。

回到點上，我們開始包這一個餃子了。首先是圓球，那么就有球心。定義好球心和半徑，這個球基本上就出來了。雖然grid上只有x和y這兩個軸，但是我們通過sqrt（radius*radius - x*x - y*y）這個方法直接模擬計算出z軸的正方向值。這樣我們在這個點上以x，y和計算出來的z三個數確定了它在模擬的三維空間中的位置，與圓心相減那就是這個位置在模擬的三維空間中應有的面法線向量surfaceNormal。這樣空間有了面法相也有了，光源light在攝像機位置也算是已知的。通過dot（light，surfaceNormal）能夠簡單的計算出光影明暗關系來。如果把這個值乘以一個指數，高光和陰影的過度能過更加的平滑。下面是代碼：

vector light = point(1, "P", 0);
light = normalize(light);

vector bmax, bmin;
getbbox(geoself(), bmin, bmax);

float posX = (@P.x - bmin.x) / (bmax.x - bmin.x);
float posY = (@P.y - bmin.y) / (bmax.y - bmin.y);

float ratio = (bmax.y - bmin.y) / (bmax.x - bmin.x);

posY *= ratio;
vector pos = set(posX, posY, 0);
vector center = set(0.5, 0.5 * ratio, 0);

float radius = 0.2;
float posZ = 0;

float distance = length(pos - center);
if(distance < radius){
        posZ = sqrt(radius * radius - pow((center.x - posX),2) - pow((center.y - posY),2));
        vector relPos = normalize(set(posX - center.x, posY - center.y, posZ));
        float bright = pow(dot(light, relPos),2);
        @Cd = set(bright, bright, bright);
}else{
        @Cd = set(1,1,1);
}

這里把點的x,y的位置壓縮到了一x軸總長為1的矩形，ratio是高寬比。這里主要是模擬glsl里面的：

vec3 pos = (gl_FragCoord.xy / resolution.xy , 0);

float ratio = resolution.x / resolution.y;

 

練習2 - 生成幾個發光的小亮點，並讓小亮點畫出三葉線。

這個比上一個練習有稍微上升了一個級別，因為我們在這用到了for循環。

同樣是每次只看一個點，for中循環每一次相當於我又要增加一個小亮點了，但是for主體里面不是在創建小亮點，在計算這個點的顏色值之前，我們已經假設所有小亮點在屏幕上什么位置和大小關系都是知道的，我們只是要求出所有亮點的顏色對這個點產生了什么影響，有點像倒推的理解。

#define PI 3.14159265359

float time = chf("time");

vector light = point(1, "P", 0);
light = normalize(light);

vector bmax, bmin;
getbbox(geoself(), bmin, bmax);

float posX = (@P.x - bmin.x) / (bmax.x - bmin.x);
float posY = (@P.y - bmin.y) / (bmax.y - bmin.y);

float ratio = (bmax.y - bmin.y) / (bmax.x - bmin.x);

posY *= ratio;
vector pos = set(posX, posY, 0);
vector center = set(0.5, 0.5 * ratio, 0);

int count = 10;
float radius = 4;
float size = 0.1;
vector color = set(0,0,0);

for(int i = 0; i < count; i++){
        vector starPos = set(sin(float(i)/float(count) * 2 * PI + time) * sin(3 * float(i)/float(count) * 2 * PI  + time) *radius,
                             cos(float(i)/float(count) * 2 * PI + time) * sin(3 * float(i)/float(count) * 2 * PI  + time) *radius,
                             0);
        float distance = size / max(length(@P - starPos), 0.001);
        color += pow(min(distance, set(1,1,1)), 3);
}

@Cd = min(color, set(1,1,1));

這里講一講hou14的一個默不作聲的升級，也就是之前提到的VOP(wrangle)比較好的改進。就是從之前可選單線程到八線程改為了有多少線程就多少線程。按理說vop在cpu中進行多線程並行運算應該是不難也合情合理的事，正如我在最上面講的vex每次是針對一個point或primitive進行計算的，同理vop也是這樣的。在模擬過程我明顯的感覺到了多線程在vop上的優勢。

測試過程我是用公司的機器做的，H14能夠同時使用最多16個線程一起運算，grid樣本一共有400*640=256,000個點。下面的這個表格是我的練習二在H13和H14中wrangle部分的計算表現：

循環次數	H14	H13
10	52.57ms	239.30ms
100	184.81ms	9.170s
1000	1.378s	1m30.06s

可以看到但循環次數成指數增長時，H14的vop中並行運算速度發揮了驚人的作用。其實這個數據也是我為什么這么想要學glsl的一個原因之一，本人非常希望能夠把feature movie中的一些特效能夠拿到實時交互中去和人產生互動，形式我先不多說，單單計算能力就已經非常吃力，而glsl給了我們一個直接能夠把GPU浮點並行運算的能力運用起來的一個途徑。之前很多對於成千上萬個點同時進行計算操作的任務，如果用傳統的方法一一排隊扔給CPU的話，效果是很難達到實時的要求的，之前關於彈力和斥力的博文就談到了這一點。

 

練習3 - 畫出一個太空穿梭效果

這里為了看看fragment shader到底有多神奇，本人直接拿着glslsandbox中的例子轉譯過來試了一下，沒想到效果確實非常好。簡單改了改一些參數直接就轉起來了，真是棒棒噠。

值得一提的是這個例子因為個人原因斷斷續續搞的，中間分別用h14和h13實現了一次，發現vex里面的數據類型在h14里面也有了一些小更新。新增了vector2 和 matrix2這兩個新數據類型，一般人可能用不上他，不過在處理平面數據上這兩個類型起始還是蠻有用的。 :P

說實話，整段代碼本人也沒搞太明白，只是隨便調了一點參數，然后把matrix2需要的計算直接寫成了手算的方法。正所謂前人栽樹后人乘涼，有時間再正經琢磨一下代碼。

#define PI 3.14159265359

float time = chf("time");

vector bmax, bmin;
getbbox(geoself(), bmin, bmax);

float posX = (@P.x - bmin.x) / (bmax.x - bmin.x) - 0.5;
float posY = (@P.y - bmin.y) / (bmax.y - bmin.y);

float ratio = (bmax.y - bmin.y) / (bmax.x - bmin.x);

posY = posY -  0.8*ratio;

posY *= ratio;
vector pos = set(posX, posY, 0);

time = time * 0.5 + ((.25+.05*sin(time*.1))/(length(pos)+.07))* 2.2;
float si = cos(time);
float co = tan(time);
//matrix2 (a b) (c d) -> (co, si) (-si, co)

float c = 1.0;
float v1 = 0.0;
float v2 = 0.0;

for(int i = 0; i < 100; i++){
        float s = float(i) * 0.035;
        vector PP = s * pos;
        PP = set(PP.x * co + PP.y * si, PP.x * (-si) + PP.y * co, 0);
        PP += set(0.22, 0.9, s - 1.9 - sin(time * 1.13) * 0.1);
        for(int j = 0; j < 8; j++){
                PP = abs(PP) / dot(PP, PP) - 0.659;
        }
        v1 += dot(PP,PP) *.0015 * (1.8+sin(length(pos*13.0)+.5-time*.2));
        v2 += dot(PP,PP) *.0015 * (1.5+sin(length(pos*13.5)+2.2-time*.3));
        c = length(PP * 0.5) * 0.35;
}

float len = length(pos);
v1 += smooth(01, 0.0, len);
v2 += smooth(1, 0.0, len);

float re = clamp(c, 0.0, 1.0);
float gr = clamp((v1 + c) * 0.35, 0.0, 1.0);
float bl = clamp(v2, 0.0, 1.0);
vector co1 = set(re, gr, bl) + smooth(0.6, 0.0, len) * 0.9;

float gray = dot(co1, set(0.499, 0.487, 0.014));
@Cd = set(gray, gray, gray);