Houdini 中Volume Retiming 放慢速度的方法

閱讀這篇博客讀者可能需要對體積的概念掌握程度比較高，如果因概念不太熟悉而造成閱讀中身體不適的，還請回頭自行補課。

設計這個解決方案的原因主要是有些鏡頭中的煙火效果可能形態很好但就是覺得速度需要放慢點。而恰恰Houdini自帶的retiming方法做插值的話是會出現類似鬼影的問題，常常結果是不能用的。昨天通過挪威同事指點后，自己也做了一個方案，這里講講我的思路。

首先是效果：

 

從圖中可以看出，雖然用這個方法仍然不能避免抖動的存在，但是整體過渡確實緩和流暢多了，傳統方法效果是整個煙霧不論速度快的還是慢的都會出現一定程度的抖動，而改進后的方法抖動的地方都是速度過快導致誤差增大而產生的。所以如果煙效的速度在一定程度的話，或者density比較厚而不太透明的話，這個插值方法是能比較准確表達出過渡的。

如下圖，我們假定當前要計算的voxel是P，我們需要計算這個voxel里面density的值從n到n+1的過渡，假設我們的substep設置成0.1，那么相當於放慢了十倍，中間需要插九幀。

我們先考慮第一個substep的情況，也就是第一個0.1的情況。因為是從n到n+1的幀的過渡，所以當前幀的基礎幀我們會用到第n幀，采樣出他的速度為vel1，通過速度乘以TimeInc和substep之后的向量尺度則相當為一個substep的跨度（這里我們只能假設幀與幀之間是線性變化的，即便肯定不是線性的）。雖然vel1在理解上肯定是向前的，但在計算voxel的時候我們一直是在討論一個固定點上采樣他的某個值，這里是density，所以我們不能用sop vop里的常用思維來思考物理位置的變化。我們先考慮n.1步的時候點P的density值最有可能是哪些位置的density轉移過來的，我們只能從兩幀中來做評斷，在n幀里面我們把vel1計算出的向量向后推一步得到A1。在判定n+1幀里面哪個位置的density最可能組成n.1的density前，我們先需要得到次P點在n+1幀的速度，假設為Vel2。根據這個速度我們可以猜出（只能是猜）n幀在點P的density很可能在n+1幀中在點PNext的位置，所以n.1的density在n+1幀中很可能就在B1的位置。同理，在n.9步中最可靠的值在n幀的C2和n+1幀的D2中。是不是很繞啊，沒辦法就通過一個速度來判定一個固定地方某個值的前世今生確實比較耗腦力。

在如何混合A1和B1這兩個density的問題上，我們還得清楚以點P為參考的計算，越是離點P近的點的值越是准確，因為速度越大，跨度就越大，這樣density在速度的積累下就越散，所以就越不精確。所以在計算n.1步時，如果是使用線性混合的方法的話，那么比較准確的方式為density = A1*0.9+B1*0.1。同理在n.9步中density = C2*0.1 + D2*0.9。

到目前為止我們的插值模型就算完成了，接着就是怎么實現的了，我簡單截一下用volume vop連的方法和直接用Vex寫的方法兩種，大同小異，計算時間也相差不大。

在鏈接volume vop之前需要得到timeshift之后的n幀和n+1幀，鏈接如下：

 

在volume vop里面的鏈接為：

 

這是Vex的方法，是通過vop的連接為原型直接轉譯寫出來的，沒什么太多好講的。

float subProcess = ch("subProcess");

float currentVelX, currentVelY, currentVelZ;
float nextVelX, nextVelY,nextVelZ;

int curAttribValueX, curAttribValueY, curAttribValueZ;
int nextAttribValueX, nextAttribValueY, nextAttribValueZ;

//find the primitive index for each input
//find the right way to get the velocity values
curAttribValueX = findattribval( @OpInput1, "primitive", "name", "vel.x");
curAttribValueY = findattribval( @OpInput1, "primitive", "name", "vel.y");
curAttribValueZ = findattribval( @OpInput1, "primitive", "name", "vel.z");

nextAttribValueX = findattribval( @OpInput2, "primitive", "name", "vel.x");
nextAttribValueY = findattribval( @OpInput2, "primitive", "name", "vel.y");
nextAttribValueZ = findattribval( @OpInput2, "primitive", "name", "vel.z");

//get the velocity of frame n
currentVelX = volumesample(@OpInput1, curAttribValueX, @P);
currentVelY = volumesample(@OpInput1, curAttribValueY, @P);
currentVelZ = volumesample(@OpInput1, curAttribValueZ, @P);
vector currentVel = set(currentVelX, currentVelY, currentVelZ);

//get the velocity of n+1 frame
nextVelX = volumesample(@OpInput2, nextAttribValueX, @P);
nextVelY = volumesample(@OpInput2, nextAttribValueY, @P);
nextVelZ = volumesample(@OpInput2, nextAttribValueZ, @P);
vector nextVel = set(nextVelX, nextVelY, nextVelZ);

//get the density sample position from frame n
vector forbackCurrentDirection, currentCheckPos;
forbackCurrentDirection = -1 * currentVel * @TimeInc * subProcess;
currentCheckPos = @P + forbackCurrentDirection;

//get the density sample position from frame n+1
vector forwardNextDirection, nextCheckPos;
forwardNextDirection = nextVel * @TimeInc * (1 - subProcess);
nextCheckPos = @P + forwardNextDirection;

//sample the density
float forbackCurrentDensity, forwardNextDensity;
int curAttribValueDensity, nextAttribValueDensity;
curAttribValueDensity  = findattribval( @OpInput1, "primitive", "name", "density");
nextAttribValueDensity = findattribval( @OpInput2, "primitive", "name", "density");

forbackCurrentDensity = volumesample(@OpInput1, curAttribValueDensity, currentCheckPos);
forwardNextDensity = volumesample(@OpInput2, nextAttribValueDensity, nextCheckPos);

//mix the density by the distance from the current sample position
float finalDensity;
finalDensity = subProcess * forbackCurrentDensity + (1 - subProcess) * forwardNextDensity;
@density = finalDensity;