monotonic cubic interpolation
Author
Hi,
I would like to implement the Monotonic cubic interpolation scheme found in "visual simulation of smoke" by Fedkiw and Stam.
I am currently using trilinear interpolation but would like to use the cubic one as it shall give visually better results.
But I have a problem understanding it. My fluid field (grid) has for example the scalar values (density) stored at the vertices. I am using them to interpolate the values to the current position within the voxel.
Ok, but the cubic interpolation scheme talks about taking dk = (fk+1 - fk-1) in the appendix. Especially the fk-1 I don´t understand. I thouhgt fk would be my first vertex of the voxel and fk+1 would be the one on the opposite site. So like this
| | /
| |/
o----o
fk fk+1
Does this now mean that fk-1 is the vertex to the left of "fk"?
Here is the hermite cubic Formula in written text. This is for one-dimensional interpolation, where fk is a list of values (vertices) defined at the location k=0...N. A value at a point t(tk,tk+1).
//////////////////////////////////////////////
//////// monotonic cubic interpolation ///////
//////////////////////////////////////////////
//ip = a3(t-tk)³ + a2(t-tk)² + a1(t-tk) + a0
//
//a3 = dk + dk+1 - Dk
//a2 = 3Dk - 2dk - dk+1
//a1 = dk
//a0 = fk
//
//dk = (fk+1 - fk-1)*0.5
//Dk = (fk+1 - fk)
//////////////////////////////////////////////
Thanks in advance for any help!
Best
Sam
[Edited by - Katachi on August 11, 2006 5:25:14 PM]
Author
Here is some code. I calculated this for x,y and z axis:
But what must I do now with these 3 values that I stored in i[0..2]? Do I need to interpolate them somehow to get the final interpolated value?
Thanks
Sam
//(t-tk)float vx = v(0) - floorf(v(0));//x axisk[0] = ((*this) ((int) (v(0) + 1), (int) (v(1)), (int) (v(2)))) - ((*this) ((int) (v(0)), (int) (v(1)),(int) (v(2)))); //fk = a0k[1] = (((*this) ((int) (v(0) + 1), (int) (v(1)), (int) (v(2)))) - ((*this) ((int) (v(0) - 1), (int) (v(1)),(int) (v(2)))))*0.5; //dk = a1if(!(Sgn(k[1])== Sgn(k[0]))) k[1] = 0.0;k[2] = (((*this) ((int) (v(0) + 2), (int) (v(1)), (int) (v(2)))) - ((*this) ((int) (v(0)), (int) (v(1)),(int) (v(2)))))*0.5; //dk+1if(!(Sgn(k[2])== Sgn(k[0]))) k[2] = 0.0;k[3] = k[0]*3 - k[1]*2 - k[2];//a2if(!(Sgn(k[3])== Sgn(k[0]))) k[3] = 0.0;k[4] = k[1] + k[2] - k[0];//a3if(!(Sgn(k[4])== Sgn(k[0]))) k[4] = 0.0;//Is this now the x-axis value?i[0] = k[4]*pow(vx,3) + k[3]*pow(vx,2) + k[1]*vx + k[0];//---------------------vx = v(1) - floorf(v(1));k[0] = ((*this) ((int) (v(0)), (int) (v(1) + 1), (int) (v(2)))) - ((*this) ((int) (v(0)), (int) (v(1)),(int) (v(2)))); //fk = a0k[1] = (((*this) ((int) (v(0)), (int) (v(1) + 1), (int) (v(2)))) - ((*this) ((int) (v(0)), (int) (v(1) - 1),(int) (v(2)))))*0.5; //dk = a1if(!(Sgn(k[1])== Sgn(k[0]))) k[1] = 0.0;k[2] = (((*this) ((int) (v(0)), (int) (v(1) + 2), (int) (v(2)))) - ((*this) ((int) (v(0)), (int) (v(1)),(int) (v(2)))))*0.5; //dk+1if(!(Sgn(k[2])== Sgn(k[0]))) k[2] = 0.0;k[3] = k[0]*3 - k[1]*2 - k[2];//a2if(!(Sgn(k[3])== Sgn(k[0]))) k[3] = 0.0;k[4] = k[1] + k[2] - k[0];//a3if(!(Sgn(k[4])== Sgn(k[0]))) k[4] = 0.0;//Is this now the y-axis value?i[1] = k[4]*pow(vx,3) + k[3]*pow(vx,2) + k[1]*vx + k[0];//---------------------vx = v(2) - floorf(v(2));k[0] = ((*this) ((int) (v(0)), (int) (v(1) + 1), (int) (v(2)))) - ((*this) ((int) (v(0)), (int) (v(1)),(int) (v(2)))); //Dk = a0k[1] = (((*this) ((int) (v(0)), (int) (v(1) + 1), (int) (v(2)))) - ((*this) ((int) (v(0)), (int) (v(1) - 1),(int) (v(2)))))*0.5; //dk = a1if(!(Sgn(k[1])== Sgn(k[0]))) k[1] = 0.0;k[2] = (((*this) ((int) (v(0)), (int) (v(1) + 2), (int) (v(2)))) - ((*this) ((int) (v(0)), (int) (v(1)),(int) (v(2)))))*0.5; //dk+1if(!(Sgn(k[2])== Sgn(k[0]))) k[2] = 0.0;k[3] = k[0]*3 - k[1]*2 - k[2];//a2if(!(Sgn(k[3])== Sgn(k[0]))) k[3] = 0.0;k[4] = k[1] + k[2] - k[0];//a3if(!(Sgn(k[4])== Sgn(k[0]))) k[4] = 0.0;//Is this now the z-axis value?i[2] = k[4]*pow(vx,3) + k[3]*pow(vx,2) + k[1]*vx + k[0];But what must I do now with these 3 values that I stored in i[0..2]? Do I need to interpolate them somehow to get the final interpolated value?
Thanks
Sam
Author
hmm, noone? I would appreciate moving this thread to maths and physics if it´s more likely to get an answer there.
Hi... i might try to answer, since no one else tried.
I checked the papers' appendix and wasn't able to follow it very well... but since they stated it was an hermite interpolation i checked wikipedia cubic hermite.
Working out the equations and placing things in a polynomial form, i got (usig the same notation in wikipedia):
P(t) = a*t^3 + b*t^2 + c*t + d
where
a = 2*P0 + m0 - 2*P1 + m1
b = -3*P0 - 2*m0 + 3*P1 - m1
c = m1
d = P0
In wikipedia, it is assumed that t varies between 0 and 1. P0 and P1 are the values at the extreme of the interval (at t=0 and t=1). The thing here are the m's (m0 and m1) that are worrying you. They are both tangents to the points P0 and P1.
I think you are correct at your initial guess. Since m0 is the tangent at P0, it is calculated taking into account P1 and P-1. With m1, it uses P0 and P2. In the paper it is the dk calculated by taking into account f_(k+1) and f_(k-1) and the same reasoning for dk+1.
So, for the given points P0 and P1, you still need 2 extra points for the cubic interpolation: P-1 and P2...
I checked the papers' appendix and wasn't able to follow it very well... but since they stated it was an hermite interpolation i checked wikipedia cubic hermite.
Working out the equations and placing things in a polynomial form, i got (usig the same notation in wikipedia):
P(t) = a*t^3 + b*t^2 + c*t + d
where
a = 2*P0 + m0 - 2*P1 + m1
b = -3*P0 - 2*m0 + 3*P1 - m1
c = m1
d = P0
In wikipedia, it is assumed that t varies between 0 and 1. P0 and P1 are the values at the extreme of the interval (at t=0 and t=1). The thing here are the m's (m0 and m1) that are worrying you. They are both tangents to the points P0 and P1.
I think you are correct at your initial guess. Since m0 is the tangent at P0, it is calculated taking into account P1 and P-1. With m1, it uses P0 and P2. In the paper it is the dk calculated by taking into account f_(k+1) and f_(k-1) and the same reasoning for dk+1.
So, for the given points P0 and P1, you still need 2 extra points for the cubic interpolation: P-1 and P2...
I've also checked some books i have here, and it follows the previous conclusions. Using either a nearest-neighbor interpolation or a trilinear interpolation inside a cell, all you need is the eight values at the 8 grid points. No need to check values that dont belong to the cell.
Using any kind of cubic interpolation, you need to pass the cells extent and go over to the adjacent cells for getting the adittional values. In this case, the values at the 8 grid points are not enough. Thats why you need the f_{k-1}, f_{k}, f_{k+1} and f_{k+2}.
Using any kind of cubic interpolation, you need to pass the cells extent and go over to the adjacent cells for getting the adittional values. In this case, the values at the 8 grid points are not enough. Thats why you need the f_{k-1}, f_{k}, f_{k+1} and f_{k+2}.
Author
Hi to the both of you,
thank you for re-checking and confirming my thoughts and all the nice explanations. That´s already good to know and then my code seems to be correct. great! :)
What I still need is, how do I get the final value for the point t(x,y,z)? I can get the values for each dimension with my code but considering that t(x,y,z) is a scalar value, it´s necessary to somehow combine these 3 values.
Any idea how? Interpolation? Sum up? anything else? What makes sense? (I am not very familiar with interpolation techniques, so this is all a bit confusing...hehe, as everything new ;) )
Thank you guys!
Best
Sam
thank you for re-checking and confirming my thoughts and all the nice explanations. That´s already good to know and then my code seems to be correct. great! :)
What I still need is, how do I get the final value for the point t(x,y,z)? I can get the values for each dimension with my code but considering that t(x,y,z) is a scalar value, it´s necessary to somehow combine these 3 values.
Any idea how? Interpolation? Sum up? anything else? What makes sense? (I am not very familiar with interpolation techniques, so this is all a bit confusing...hehe, as everything new ;) )
Thank you guys!
Best
Sam
First, there is no both of you... both answers were from the same person :D
Also, i made a mistake in one the indexes in the formulas i've worked up from wikipedia when i posted here. The correct formula is:
P(t) = a*t^3 + b*t^2 + c*t + d
where
a = 2*P0 + m0 - 2*P1 + m1
b = -3*P0 - 2*m0 + 3*P1 - m1
c = m0
d = P0
Since i've never worked with hermite interpolation in 3d data, i don't know how to exactly make a correct reconstruction in the general case for 3d from the 1d case. Perhapas someone that is more confortable with this math can answer this.... However, there seems to be some errors in your papers appendix formulas:
http://www-static.cc.gatech.edu/~sarah/smoke.htm
Also, the usual paper that analizes 3d reconstruction with 3d data is "An Evaluation of Reconstruction Filters for Volume Rendering". If you google for it, you'll find it. They analizes several reconstruction filters, some of them being cubic and separable (i.e. h(x,y,z)=h(x)*h(y)*h(z)) that simplify the calculations. You can always switch to one of them if no one else helps.
Also, i made a mistake in one the indexes in the formulas i've worked up from wikipedia when i posted here. The correct formula is:
P(t) = a*t^3 + b*t^2 + c*t + d
where
a = 2*P0 + m0 - 2*P1 + m1
b = -3*P0 - 2*m0 + 3*P1 - m1
c = m0
d = P0
Since i've never worked with hermite interpolation in 3d data, i don't know how to exactly make a correct reconstruction in the general case for 3d from the 1d case. Perhapas someone that is more confortable with this math can answer this.... However, there seems to be some errors in your papers appendix formulas:
http://www-static.cc.gatech.edu/~sarah/smoke.htm
Also, the usual paper that analizes 3d reconstruction with 3d data is "An Evaluation of Reconstruction Filters for Volume Rendering". If you google for it, you'll find it. They analizes several reconstruction filters, some of them being cubic and separable (i.e. h(x,y,z)=h(x)*h(y)*h(z)) that simplify the calculations. You can always switch to one of them if no one else helps.
Author
arg, yes, it was only wolverine. :) Sorry.
Thanks for the fix info, got it in there now. Thanks for the links too! Will check out that paper.
Sam
Thanks for the fix info, got it in there now. Thanks for the links too! Will check out that paper.
Sam
This topic is closed to new replies.
Advertisement
Popular Topics
Advertisement