Houdini 中 Gray Scott Reaction-Diffusion算法的實現

這篇文章是吧很久以前學的一個神奇算法歸一下檔，在公交車上突然想起來了，覺得還是很有必要再仔細梳理一遍，對以后也許有用。

先看圖再說話：

Gray Scott Reaction-Diffusion算法， 在模擬微觀細胞的運動或者類似的效果是非常神奇。

理論鏈接：http://www.karlsims.com/rd.html

原理：模擬兩種物質之間在平面（暫時是平面）上的相互作用，動作分為反應與擴散。

公式：

右手邊第一項為擴散，第二項為反應，第三項是供給，DA與DB為分散率，▽2A或B為相鄰位置擴散陣列（下圖），f和k是兩個非常重要的外參數，分別叫做feed和kill。可以簡單理解成A是環境中的實物，兩個B和一個A就能變成三個B，過程中feed會變出更多的A，k也會去掉一些B。非常想現實的生長環境。

|0.05| 0.2 |0.05|

| 0.2 | -1  | 0.2 | 

|0.05| 0.2 |0.05|    

上面是拉普拉斯方程用到的陣列，因為外圍總和等於中間數的反數，所以參數化后便可以隨便調

 

在houdini里面，使用sop solver來迭代A和B的物質變化，在進入solver前，先定義好A和B的值，一般是A在整個平面是1.0，B隨機散布開，這樣出來的效果也好,如下圖是我做的B隨機散布示意圖：

同時用noise方法隨機分布feed和kill的值，你當然可以使用固定的f和k，不過有隨機變化的話會讓生長更加自然有趣。下面是detail面板里面的簡單示意：

值得注意的是feed和kill的值是非常敏感的，我自己測試時，也就幾個非常窄的區間內會有生長效果，要不然cb會非常快的全死掉。

在solver里面給preframe 加一個新的vex node，代碼如下：

#pragma label resx "Res X"
#pragma label resy "Res Y"

#pragma label diff_a "Diffusion A"
#pragma label diff_b "Diffusion B"

#pragma label weight "Weight"
#pragma range weight 0.125 0.25

#pragma label feed "Feed"
#pragma range feed 0 1
#pragma label kill "Kill"
#pragma range kill 0 1
#pragma label time_step "Sub Step"
#pragma range time_step 0 20

#pragma hint ca invisible
#pragma hint cb invisible


//this fuction will diffuse 
float laplacian(int ptnum; int resx; int resy; string attrib; float weight){

        // | w_s | w | w_s |
        //-----------------------
        // | w |-1 | w |
        //-----------------------
        // | w_s | w | w_s | weight

        // | pt+w-1 | pt+w | pt+w+1 |
        //-----------------------------
        // | pt-1 | pt | pt+1 |
        //-----------------------------
        // | pt-w-1 | pt-w | pt-w+1 | location

        int vb = resx*2, hb = 1;

        // detects boundaries
        if((ptnum % resx) == 0 || (ptnum % resx) == (resx - 1)){
                hb *= -1;  // horizon
        }

        if(ptnum < resx || ptnum > (resy-1)*resx){
                vb *= -1;  // vetical
        }



        float n[];
        resize(n,9);

        float weight_second = 0.25 - weight;

        n[0] = point(0, attrib, ptnum - vb - 1) * weight_second;  // left bottom
        n[1] = point(0, attrib, ptnum - vb    ) * weight;         // bottom
        n[2] = point(0, attrib, ptnum - vb + 1) * weight_second;  // right bottom
        n[3] = point(0, attrib, ptnum      - 1) * weight;         // left
        n[4] = point(0, attrib, ptnum         ) * (-1.0);           // center
        n[5] = point(0, attrib, ptnum      + 1) * weight;         // right
        n[6] = point(0, attrib, ptnum + vb - 1) * weight_second;  // left top
        n[7] = point(0, attrib, ptnum + vb    ) * weight;         // top
        n[8] = point(0, attrib, ptnum + vb + 1) * weight_second;  // right top

        float sum=0.0;
        foreach(float i; n){
                sum += i;
        }

        return sum;
}

sop
gray_scott(
        int resx = 512;
        int resy = 512;
        float diff_a = 1.0;
        float diff_b = 1.0;
        float weight = 0.15;
        export float feed = 0.05;
        export float kill = 0.05;
        float time_step = 1.0;
        export float ca = 0.0;
        export float cb = 0.0;
        )
{
                float react_rate = ca * pow(cb,2);
                float feed_set = feed * (1.0 - ca);
                float kill_set = (kill + feed) * cb;
                //float substep = 1.0/(float)time_step;

        ca += (diff_a * laplacian(ptnum, resx, resy, "ca", weight) - react_rate + feed_set) * time_step;
        cb += (diff_b * laplacian(ptnum, resx, resy, "cb", weight) + react_rate - kill_set) * time_step;
}

 最后看一看不同的feed和kill值產生的不同效果：

 feed = 0.051  kill = 0.062

 feed = 0.0367 kill = 0.069

 feed = 0.0561 kill = 0.0636