Rotation of the coordinates of a 3D function

[Amentia]

Hello,

I come across a problem in programming and I do not find a lot of help on the internet, so I hope I can find an answer here. I have a 3D array representing a function, say f(i,j,k) and a basis function u(i,j,k). I would like to perform a general rotation of the basis function u so that I can project my function f to any axis I want.

I do not really know how to do that, I am afraid the coordinates will become non integers and I don't know how to deal with this. I was told that it is possible to interpolate with spline functions but I am not familiar with this. Is it a good solution? Are there better ways to tackle this problem?

Best regards!

 

[jedishrfu]

Perhaps a google search will help here:

https://www.google.com/search?q=Rot...rome..69i57.3437j0j8&sourceid=chrome&ie=UTF-8

 

[Amentia]

Hello thank you for your answer but I do not find which links address my problem. To be more precise, I mean that for example mathematically I take a function f(x, y, z) and I make some rotation of x, y, z giving new coordinates x', y', z'. I want to create a new function g(x', y', z') =f(x, y, z). But I could have x'= 1/sqrt(2)*(x+y) and indexes in arrays are integers. So I want to associate a value to g(0.5,1.5,10.3) and I don't know how to do this.

 

[Stephen Tashi]

So I want to associate a value to g(0.5,1.5,10.3) and I don't know how to do this.

Your goal isn't clear.

Is the function $f$ a function of the coordinate indices such as $f(i,j,k) = 3i + 2j - k$?

Or is the function $f$ an array so its value is given by information like $f[ i][j][k] = 23.944$ .

Or is the function $f$ a function of data in other arrays like $f(i,j,k) = 3.6* m[ i] + 2.8*w[j][k]$ ?

 

[Amentia]

It is given like f[ i] [j] [k] =23.944.

 

[Stephen Tashi]

It is given like f[ i] [j] [k] =23.944.

Since you are talking about rotating the indices, do the indices i,j,k represent a position is space? If so, how do you determine what position (i,j,k) represents? For example, is (2,3,5) the location (2.0, 3.0, 5.0) in a cartesian coordinate system. Or is (2,3,5) the location ( X[2],Y[3],Z[5]) where you have arrays that tell you the coordinates?

(The forum's software interprets an "i" in square brackets as a command to display text in italics, so when you want to use "i" to denote a bracketed array index, you must leave a space in front of the "i" - like "[ i]".)
Mentor note: I fixed the problem described above.

 

[Amentia]

Thank you for your answer. In my case I, j, k are indeed a position, but if I could create some arrays like that to get real coordinates, perhaps it would help to rotate the indices? I understand I, j, k as coordinates in a cubic grid like Cartesian coordinates but in a discrete version. So position I, j, k is a point in space and we associate a value to this point which is stored in the array.

 

[.Scott]

Converting from x,y,z, to x',y',z', you have:
x' = ax+by+cz+d
y' = ex+fy+gz+h
z' = ix+jy+kz+l

So you have a matrix Ur:
a,b,c
e,f,g
h,i,j

And Ud = d, h, l

So x',y',z' = [x,y,z]*Ur+Ud

Invert the Ur matrix to get U'r:
a',b',c'
e',f',g'
h',i',j'

And compute U'd with a matrix multiply: U'd = -Ud*Ur

Now you have x,y,z = [x',y',z']*U'r+U'd

Now your g(x',y',z') = f(x,y,z) becomes f([x',y',z']*U'r+U'd)
But, of course, your x,y,z will not be integers - and they may not be positive or within the array bounds of f.. So you will need to clip and interpolate within the f array to get the entries in your g array.

Let's say g(3,5,4) ends up as f(2.6,-4.1,3.3) - then you will use some default value, since you are outside the range of f.
Let's say g(4,5,7) ends up as f(12.9,1.3,2.8), but 12.9 is greater that the first array dimension in f - then again, you will use some default value, since you are outside the range of f.

Let's say g(2,3,4) ends up as f(5.1,3.7,4.9) - that point lands within a cube in f. So interpolate the value of the 8 corners of the cube to get you g () value.

 

[Stephen Tashi]

but if I could create some arrays like that to get real coordinates, perhaps it would help to rotate the indices?

That depends on whether you want g to be a function of integer values or whether g is a function of 3 floating point variables.

If g[2][3][5] is to represent the value of data exactly at (2,3,5) then you have the problem that (2,3,5) can be a point that came by way of rotation from a point that is not necessary given by whole numbers. You can un-rotate (2,3,5) back to a point in original coordinate system like (1.3, 2.8,...). Then you must find a way to interpolate the value of f at (1.3, 2.8,..). There are various ways to interpolate data. If you use this method you will perform the inverse rotation of the rotation you want to perform on the coordinate system used by f.

 

[Amentia]

Hi .Scott, thank you for you answer.

Now your g(x',y',z') = f(x,y,z) becomes f([x',y',z']*U'r+U'd)
But, of course, your x,y,z will not be integers - and they may not be positive or within the array bounds of f.. So you will need to clip and interpolate within the f array to get the entries in your g array.

Let's say g(3,5,4) ends up as f(2.6,-4.1,3.3) - then you will use some default value, since you are outside the range of f.
Let's say g(4,5,7) ends up as f(12.9,1.3,2.8), but 12.9 is greater that the first array dimension in f - then again, you will use some default value, since you are outside the range of f.

Let's say g(2,3,4) ends up as f(5.1,3.7,4.9) - that point lands within a cube in f. So interpolate the value of the 8 corners of the cube to get you g () value.

This is exactly what I have never done. Do you mean that I will lose information after doing the transformation by giving some default values? And when you talk about interpolating, is it by use of some advanced numerical methods? Or in this case some kind of average? But I do not think I will have definite values at integer positions in the new grid.

That depends on whether you want g to be a function of integer values or whether g is a function of 3 floating point variables.

If g[2][3][5] is to represent the value of data exactly at (2,3,5) then you have the problem that (2,3,5) can be a point that came by way of rotation from a point that is not necessary given by whole numbers. You can un-rotate (2,3,5) back to a point in original coordinate system like (1.3, 2.8,...). Then you must find a way to interpolate the value of f at (1.3, 2.8,..). There are various ways to interpolate data. If you use this method you will perform the inverse rotation of the rotation you want to perform on the coordinate system used by f.

Mathematically I would like g to be a function of floating point variables. But I need to store it in an array, so it looks impossible. Could you elaborate on the various ways to interpolate data that could work in my case? Or even in simpler cases like having a 2D grid with values associated (for example the colors of a picture) and rotating this grid so that the picture becomes rotated as well.

 

[.Scott]

There's a wiki article that describes the different types of interpolation:
https://en.wikipedia.org/wiki/Interpolation

But I will show you an example using simple linear interpolation:

Say you want f(2.3,4.6,6.9)
lololo = f(2,4,6)
lolohi = f(2,4,7)
lohilo = f(2,5,6)
lohihi = f(2,5,7)
hilolo = f(3,4,6)
hilohi = f(3,4,7)
hihilo = f(3,5,6)
hihihi = f(3,5,7)

So f(2.4,4.6,6.8) will be
0.7*(0.4*(0.1*lololo+0.9*lolohi)+0.6(0.1*lohilo+0.9*lohihi))
+0.3*(0.4*(0.1*hilolo+0.9*hilohi)+0.6(0.1*hihilo+0.9*hihihi))

So the closer you get to a corner, the closer your value will be to the value at that corner.

 

[.Scott]

Do you mean that I will lose information after doing the transformation by giving some default values?

No. The default value is what you had all along. If you were not rotating, what would you have returned is you needed f(-7,3,-3)? I assume you would not have just allowed the array bounds to be violated. Perhaps you would have asserted an exception. Or perhaps you would have provided a default value. Whatever you did before, do the same now.

But I do not think I will have definite values at integer positions in the new grid.

I was telling you how to put definite values at integer positions in the new grid (a g array).
But you don't really even need a new grid - since what I gave you will allow you to transform x',y',z' to x,y,x and then look it up in the f array with interpolation. So you really don't need to create the g array - unless you want to.

 

[.Scott]

Or even in simpler cases like having a 2D grid with values associated (for example the colors of a picture) and rotating this grid so that the picture becomes rotated as well.

Here is the best example I found. Although he doesn't use a matrix for rotation:
http://polymathprogrammer.com/2008/10/06/image-rotation-with-bilinear-interpolation/

Hope it's worth the effort to read through it.

 

[Stephen Tashi]

Mathematically I would like g to be a function of floating point variables.

The same general idea would apply. You can un-rotate (2.6,3.8,5.9) to obtain a point in the coordinate system used by f. Then you must interpolate the data at that point.

But I need to store it in an array, so it looks impossible.

You could have a function of floating point variables that looked-up values in an array. That would still involve interpolating.

Or even in simpler cases like having a 2D grid with values associated (for example the colors of a picture) and rotating this grid so that the picture becomes rotated as well.

That's a special situation and I don't know whether it applies to you data. Data for picture at (i,j) usually represents data for the area of a pixel. To find find data at (i', j') in the rotated image, a good way to do it to imagine the 4 (floating point) vertices of the pixel that has data for (i',j') unrotated so they sit on the orginal image. The unrotated pixel is a rectangle (or square) that sits at an angle to the pixels on the original image and it usually contains parts of several pixels on the original image. So, in principle, you can take the average of the data in the original pixels, weighting each value by the area that that the unrotated pixel covers. (I'm sure rotating images is a well studied process and we can look-up quick ways of doing it.)

On the other hand, your data at (i,j) might represent data only at that one location instead of for an area around it. What exactly is the data?

 

[Amentia]

So if I understand well so far, here is the situation. I define a rotation matrix R :

$(x',y',z') = R (x,y,z)$

Now I have my well defined function f(i,j,k) and I would like to create a rotated function g(i,j,k). So to keep the indices of g as integers, I can say:

$g(R(x,y,z)) = f(x,y,z) \Leftrightarrow g(x,y,z) = f(R^{-1}(x,y,z))$

From there, the problem is that f has now real indices. So I need to make some for loops such as, in a cube of length n:
$\text{for } k=1,n\\ \text{for } j=1,n\\ \text{for } i=1,n\\ (x,y,z) = R^{-1}(i,j,k)\\ \text{find the 8 integers indices closest to x,y,z}\\ g[i,j,k] = interpolation(f(int[x], int[y], int[z]), f(int[x], int[y]+1, int[z]), f(int[x], int[y], int[z]+1), etc.)$

And here I have to select a suitable interpolation function. I know there is a b-spline function in fortran but I don't know if it is easy to use there or if my problem is simple enough to write one from scratch.

Also I did not understand well this part:

No. The default value is what you had all along. If you were not rotating, what would you have returned is you needed f(-7,3,-3)? I assume you would not have just allowed the array bounds to be violated. Perhaps you would have asserted an exception. Or perhaps you would have provided a default value. Whatever you did before, do the same now.

I did not do anything before because I did not write the original code, but I need to perform this rotation to make it suitable for my goals.

Here is the best example I found. Although he doesn't use a matrix for rotation:
http://polymathprogrammer.com/2008/10/06/image-rotation-with-bilinear-interpolation/

Hope it's worth the effort to read through it.

Thank you for the link, I have not read it yet but it looks nice.

 

[Amentia]

On the other hand, your data at (i,j) might represent data only at that one location instead of for an area around it. What exactly is the data?

The data is some wave functions, similar to atomic orbitals for example, which undergo a lot of transformations before like Fast Fourier Transform, so that I could even think about doing the rotation before in the k-space, but I am not sure it makes things simpler. So maybe from the way the functions look, it is similar enough to the 2D images.

 

[.Scott]

From there, the problem is that f has now real indices. So I need to make some for loops such as, in a cube of length n:
$\text{for } k=1,n\\ \text{for } j=1,n\\ \text{for } i=1,n\\ (x,y,z) = R^{-1}(i,j,k)\\ \text{find the 8 integers indices closest to x,y,z}\\ g[i,j,k] = interpolation(f(int[x], int[y], int[z]), f(int[x], int[y]+1, int[z]), f(int[x], int[y], int[z]+1), etc.)$

And here I have to select a suitable interpolation function. I know there is a b-spline function in fortran but I don't know if it is easy to use there or if my problem is simple enough to write one from scratch.

Only the simplest interpolation algorithms work with just the closest 8 points (or closest 2 points in 1D). Most required 4 points in 1D, which could be as many as 64 in 3D.
I would implement with the simple linear interpolation that I described - then move on to a more sophisticated transform if it was needed.

Also I did not understand well this part: <default value>
I did not do anything before because I did not write the original code, but I need to perform this rotation to make it suitable for my goals.

You haven't told us what the purpose of all this is - so I don't know for certain whether there was ever a problem with violating the array bounds of the f array.
But now, you are computing indeces to the f array, so you definitely need to check those indeces.
So you tell me: What are you going to do is the indeces you calculate come out to -4.4, 5.5, -6.6?
If you can answer that question, you should be able to understand what I was saying.[/QUOTE]

 

[.Scott]

The data is some wave functions, similar to atomic orbitals for example, ...

OK. So these objects are surrounded by "white space". So if you try to exceed your array bounds, you want f() to return whatever corresponds to the white space.

 

[Amentia]

OK. So these objects are surrounded by "white space". So if you try to exceed your array bounds, you want f() to return whatever corresponds to the white space.

In fact they are not surrounded by white space but I think in the end, it is the same. If we take a colour scale where blue is some small value close to 0 and red is a high value close to 1, the first ball in the image will be red in the centre and then will start going orange, yellow, green and finally blue after its boundary. And then it is not suddenly white but it should be a uniform colour. So if I now understand, since the point (0,0,0) in my case will always have a low value, I could choose to assign this value to the function whenever the problem occurs.

I don't know if this will happen, but what if an important point exceed the array bounds? There is nothing to do about this?

 

[.Scott]

In fact they are not surrounded by white space but I think in the end, it is the same. If we take a colour scale where blue is some small value close to 0 and red is a high value close to 1, the first ball in the image will be red in the centre and then will start going orange, yellow, green and finally blue after its boundary. And then it is not suddenly white but it should be a uniform colour. So if I now understand, since the point (0,0,0) in my case will always have a low value, I could choose to assign this value to the function whenever the problem occurs.

I don't know if this will happen, but what if an important point exceed the array bounds? There is nothing to do about this?

It's going to happen - so what do you want to do about it? It's your program. What do you want to do.
Instead of using f(0,0,0) as your default, perhaps you should used f( max(0,min(xmax,x)), max(0,min(ymax,y)), max(0,min(zmax,z)) ).
Perhaps you want to stop the program and issue an error?
If you don't know what the requirements are, you're not going to get them from PhysicsForum - you would need to go to the stakeholders.
Who are your stakeholders? Who will use or judge the final result?

Also, when you described this as a simple "rotation", I let it slide. If you notice in my earlier description, I not only provided a rotation matrix (Ur) but a translation vector as well (Ud). If you don't have that translation vector, you will always be rotating about 0,0,0. And that will certainly give you negative numbers - and worse yet, you may not even scan through the actual region with the electron shell so you will have nothing to display (or whatever you are doing with the data).

 

[Amentia]

It's going to happen - so what do you want to do about it? It's your program. What do you want to do.
Instead of using f(0,0,0) as your default, perhaps you should used f( max(0,min(xmax,x)), max(0,min(ymax,y)), max(0,min(zmax,z)) ).
Perhaps you want to stop the program and issue an error?
If you don't know what the requirements are, you're not going to get them from PhysicsForum - you would need to go to the stakeholders.
Who are your stakeholders? Who will use or judge the final result?

I am the only one who will use it and if it does not work, nobody else will be concerned in the near future. But I really would like to make it work because it is a good tool for the analysis in my work. So there is no specific requirements, but it would be important to make this rotation as good as possible. If half of the relevant points are lost each time I use the code, maybe it is not worth it to continue. So that is why I am asking here. I am not sure we are on the same wavelength because I did not mean to ask about requirements that stakeholders want from me, only how to do best this rotation and avoid losing important points if it is possible. Most of my functions will look quite similar and have important points in the same zone when nothing unexpected happens. So I came here to know if what I want to do is not out of reach. Raising exceptions 90% of the time might not help me because the calculations performed after will be wrong.

Also, when you described this as a simple "rotation", I let it slide. If you notice in my earlier description, I not only provided a rotation matrix (Ur) but a translation vector as well (Ud). If you don't have that translation vector, you will always be rotating about 0,0,0. And that will certainly give you negative numbers - and worse yet, you may not even scan through the actual region with the electron shell so you will have nothing to display (or whatever you are doing with the data).

Yes, sure, you are right, I was thinking to use the centre of my cube to translate the origin. Usually the important points will be closer to the centre but some can be a bit far as well.

 

[.Scott]

I am the only one who will use it and if it does not work, nobody else will be concerned in the near future. But I really would like to make it work because it is a good tool for the analysis in my work. So there is no specific requirements, but it would be important to make this rotation as good as possible.

OK. So you are the only stakeholder - and you are looking for suggestions. My suggestion would be to use the max/min code I showed in the last post. If it works, it works. And when it doesn't work, it should be easy to tell what is wrong. But after you're "done", experiment with it and see for yourself.

If half of the relevant points are lost each time I use the code, maybe it is not worth it to continue.

You don't want to loose half the points. When I said "It's going to happen", I didn't mean every time. I meant that in a lot of cases, you're going to specify rotation/translation parameters that don't do what you expect - or they do do what you expect but they clip a lot of data. So, in that sense, "It's going to happen" and therefore, you need to write code to handle the situation.

I am not sure we are on the same wavelength because I did not mean to ask about requirements that stakeholders want from me, only how to do best this rotation and avoid losing important points if it is possible.

I was using the terms commonly used in software design. One of my main purposes was to introduce you to the terms. Hopefully you can appreciate that people often start software projects without a firm grasp of what the requirements are. Often, that can not be avoided.

Yes, sure, you are right, I was thinking to use the centre of my cube to translate the origin. Usually the important points will be closer to the centre but some can be a bit far as well.

Let's say that range of your cube is f(0,0,0) to f(20,20,20). Then the center point will be 10,10,10. So you could go (x',y',z') = (((x,y,z)-(10,10,10)) * U ) + (10,10,10). That way the center of the cube will not move. Of course, it will still be possible to stay within you g() cube while running out of the bounds of your f() cube.

 

[Amentia]

Thank you very much for your help. I will look into that tomorrow and I will come back here if I have new questions.

 

[Tom.G]

One way to avoid running out-of-bounds is to size the (perhaps conceptual) g cube such that its major diagonal always fits 'inside' the f cube regardless of rotation. That still leaves the translation problem un-solved though.