trilinear interpolation

Started by
12 comments, last by tompish 14 years, 1 month ago
If i have an array like this: #define IX(i,j,z) ((i) + (j) * (*cube.N) + (z) * (*cube.N) * (*cube.N)) And three variables u[IX(i,j,k)] v[IX(i,j,k)] w[IX(i,j,k)] how do i implement trilinear interpolation on these three components?
Advertisement
Quote:Original post by tompish
how do i implement trilinear interpolation on these three components
Not at all!

You need to have 8 elements for tri-linear interpolation, as well as 3 weights. A possibility would be: You may have 3 co-ordinates { x, y, z } each one a floating point value. Then you can compute indices on an integer lattice e.g. as
{ i, j, k } := floor( { x, y, z } )
as well as weights
{ u, v, w } := { x, y, z } - { i, j, k }

You can then fetch the said 8 values from your field by using all combinations of { i, j, k } and { i+1, j+1, k+1 }. Then perform the 4 linear interpolations in the first dimension with the 8 field values, then the 2 linear interpolations in the 2nd dimension with the 4 results from the previous interpolation, and then the 1 linear interpolation in the 3rd dimension with the 2 results of the previous interpolation. That's it.

(Well, I had also hinted at this wikipedia page, but that would be too easy ;))
I belive that this piece of code does what you said... but something is wrong here because my program is not working korrectly

dt0=*dt*(*N);

for(k = 1; k <= *N; k++)
{
for(j = 1; j <= *N; j++)
{
for(i = 1; i <= *N; i++)
{
x = i-dt0*u[IX(i,j,k)];
y = j-dt0*v[IX(i,j,k)];
z = k-dt0*w[IX(i,j,k)];

if (x<0.5f)
x=0.5f;

if (x>*N+0.5f)
x=*N+0.5f;

i0=(int)x;
i1=i0+1;

if (y<0.5f)
y=0.5f;
if (y>*N+0.5f)
y=*N+0.5f;

j0=(int)y;
j1=j0+1;

if (z<0.5f)
z=0.5f;
if (z>*N+0.5f)
z=*N+0.5f;

k0=(int)z;
k1=k0+1;


s1 = x-i0;
s0 = 1.0f-s1;
t1 = y-j0;
t0 = 1.0f-t1;
u1 = z-k0;
u0 = 1.0f-u1;

d[IX(i, j, k)] = s0 * ( t0 * (u0 * d0[IX(i0, j0, k0)]
plus u1 * d0[IX(i0, j0, k1)])
plus ( t1 * (u0 * d0[IX(i0, j1, k0)]
plus u1 * d0[IX(i0, j1, k1)])))
plus s1 * ( t0 * (u0 * d0[IX(i1, j0, k0)]
plus u1 * d0[IX(i1, j0, k1)])
plus ( t1 * (u0 * d0[IX(i1, j1, k0)]
plus u1 * d0[IX(i1, j1, k1)])));

}
}
}

PS. Sory for the horrible code dont realy know how to do nice code in the messege field.
Quote:PS. Sory for the horrible code dont realy know how to do nice code in the messege field.
You can enclose it in source tags, like this:

[source]...code goes here...[/source]
Oh thanks a lot.. here is the code for my interpolation:

void Semiadvect(int *N, int b, float *d, float *d0, float *u, float *v, float *w, float *dt){  int i0, j0, i1, j1, k0, k1, i ,k, j;  float dt0, x, y, z, s1, t1, u1, s0, t0, u0;  dt0=*dt*(*N);  for(k = 1; k <= *N; k++)  {    for(j = 1; j <= *N; j++)    {       for(i = 1; i <= *N; i++)       {	 x = i-dt0*u[IX(i,j,k)];	 y = j-dt0*v[IX(i,j,k)];	 z = k-dt0*w[IX(i,j,k)];	 if (x<0.5f)	   x=0.5f;         if (x>*N+0.5f)	   x=*N+0.5f;         i0=(int)x; i1=i0+1;         if (y<0.5f)	   y=0.5f;	 if (y>*N+0.5f)	   y=*N+0.5f;         j0=(int)y; j1=j0+1;         if (z<0.5f)	   z=0.5f;         if (z>*N+0.5f)	   z=*N+0.5f;         k0=(int)z; k1=k0+1;	 s1 = x-i0;	 s0 = 1.0f-s1;	 t1 = y-j0;	 t0 = 1.0f-t1;	 u1 = z-k0;	 u0 = 1.0f-u1;         d[IX(i, j, k)] =                     s0 * ( t0 * (u0 * d0[IX(i0, j0, k0)]                                +u1 * d0[IX(i0, j0, k1)])                        +( t1 * (u0 * d0[IX(i0, j1, k0)]                                +u1 * d0[IX(i0, j1, k1)])))                   +s1 * ( t0 * (u0 * d0[IX(i1, j0, k0)]                                +u1 * d0[IX(i1, j0, k1)])                        +( t1 * (u0 * d0[IX(i1, j1, k0)]                                +u1 * d0[IX(i1, j1, k1)])));          }       }    }    bounds (N, b, d);} 


Am i missing somthing or doing somthing wrong?

[Edited by - tompish on February 28, 2010 12:46:46 PM]
Quote:Original post by tompish

Oh thanks a lot.. here is the code for my interpolation:

*** Source Snippet Removed ***

Am i missing somthing or doing somthing wrong?
For me at least, that code would be difficult to proofread due to the poor formatting.

I'm not sure exactly how the forum works in this respect, but I've noticed that when posting code, using soft tabs (e.g. spaces) usually gives the best results overall. (I think maybe mixing sort and hard tabs can give the worst results, since depending on how the tabs are interpreted you may get inconsistent indentations.)

What you might try is converting all of your tabs to spaces, and then editing the code further (if necessary) to make sure that it's easy to read and that expressions are relatively short and/or grouped in a clear way (for example, that long expression that computes the result would be a real chore to dig through as is).

Anyway, that's just me - maybe someone else will be able (or willing) to look at the code as is and give you some feedback.
Does this look better you think?
Quote:Original post by tompish
Does this look better you think?
Yes, much :)

However (and I'm not just trying to be difficult here), I'm still having trouble making sense of the code. I could probably sort through it with enough time, but it's usually easier to answer questions like this one when more context is provided.

It looks to me like there's more going on here than just a simple trilinear interpolation. What exactly is the problem you're trying to solve? What do all those cryptically named variables represent? What is the mysterious cube.N? Why do your indices go from 1 to N rather than from 0 to N - 1? You said it's not working correctly - what does that mean exactly? Does it crash? Is the output wrong? If the latter, what's the input, what's the expected output, and what's the actual output? What does the debugger tell you?

From looking at your code, I think a good first step would be to do some refactoring. That would probably make it easier for us to understand it, and you might even find that in the process of refactoring the code, you end up solving the problem (whatever it might be).
First of all, I second jyk's suggestion to do some refactoring.

AFAIS, the Cube.N (or just N as it is named now) is the amount of sample points along each dimension of the cubic lattice with which the value fields are sampled. The opening post let me assume this due to the way Cube.N is incorporated into the array index linearization code. However, why the counter { i, j, k } range from 1 to N is incomprehensible to me, too.

I also don't understand why simple data types that are not altered by the function are specified by pointer (I mean *N and *dt). If, for any reason I can't see ATM, it is necessary to use pointers, then use at least the const qualifier if the language allows it.

I have no real clue what { u, v, w } and dt are good for. They seem to store sampled scalar fields from which the real co-ordinates into the original field are computed. I can't say something about this, but I recommend to drop that for the reason to verify the interpolation code in the first instance.

I also don't understand the limitation of the co-ordinates to the interval [ 0.5, *N+0.5 ]. Let's use x and i0, i1 as example. If x is limited to [ 0.5, *N+0.5 ], then i0 is limited to [ 0, *N ] and i1 is limited to [ 1, *N+1 ], aren't they? But that implies that *N+2 samples along each dimension are available, what doesn't seem to fit the other code.

It would also be nice to have some comments, and variable names that express a meaning! That would make reading and understanding the code much easier. Please bear in mind that also you, as the programmer, will have trouble to understand your own code after some few monthes of progression.

Besides all that, the actual interpolation code seems me correct (after a pure visual inspection).
I don't know what the problem is, but I get only empty posts from you, tompish.

This topic is closed to new replies.

Advertisement