# trilinear interpolation

This topic is 2844 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

If you intended to correct an error in the post then please contact us.

## Recommended Posts

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?

#### Share this post

##### Share on other sites
Quote:
 Original post by tompishhow 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 ;))

#### Share this post

##### Share on other sites
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.

#### Share this post

##### Share on other sites
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]

#### Share this post

##### Share on other sites
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]

#### Share this post

##### Share on other sites
Quote:
 Original post by tompishOh 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.

#### Share this post

##### Share on other sites
Does this look better you think?

#### Share this post

##### Share on other sites
Quote:
 Original post by tompishDoes 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).

#### Share this post

##### Share on other sites
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).

#### Share this post

##### Share on other sites
I don't know what the problem is, but I get only empty posts from you, tompish.

#### Share this post

##### Share on other sites
Dont know why but i cant post a reply from any computer so I am assuming that i am writing a too long message so here is the explenation of what i am working on and the next reply will be the code with comments.

I am working on stable fluids introduced by Jos Stam. Trying to implement into 3d. The representation of the velocity field is not correct in my program and I think that this part of code is not doing the job correctly. The reason I use pointer is because i have a class defined for my physics solver cube so i find it nescesary to use pointers.

#### Share this post

##### Share on other sites
link to my code cause i cant seem to post my code for some damn reason.

http://pastebin.com/vQAxpnhZ

#### Share this post

##### Share on other sites
Quote:
 Orignal post by tompishThe reason I use pointer is because i have a class defined for my physics solver cube so i find it nescesary to use pointers.
At the end the decision is up to you, but I don't see the reason. Somewhere the routine Semiadvect is invoked, and that invocer should be able to read the values from Cube and put them as parameters for Semiadvect. Inside Semiadvect there is no reason to change those values; you actually don't change them there, but there is no documentation about this.

Unfortunately, I don't know anything about Jos' paper. But the main questions are still there: Have you thought about the 1..N and N+2 size problems? And have you tried to drive the tri-linear interpolation isolated, to make sure it works for itself?

Do you use this paper?

#### Share this post

##### Share on other sites
Nope have not tried to isolate the interpolation, that was a good thing to do ill try that later on. And yes I have thought about the grid size.. the actual grid is 1..N while there is an extra layer of cells that account for the boundary conditions 0 and N+1.

tx for the help man.

#### Share this post

##### Share on other sites

This topic is 2844 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

If you intended to correct an error in the post then please contact us.

## Create an account or sign in to comment

You need to be a member in order to leave a comment

## Create an account

Sign up for a new account in our community. It's easy!

Register a new account

## Sign in

Already have an account? Sign in here.

Sign In Now

• ### Forum Statistics

• Total Topics
628698
• Total Posts
2984275

• 20
• 10
• 13
• 13
• 11