High Performance Interpolation of 3D grid data

Hello everybody,

currently working on a project in which magnetic field (grid) data is rotated and shifted in 3 dimensions and afterwards evaluated at specific points in the COS.

So, what I have, is x-, y-, z- Coordinates with the corresponding values for Bx, By, Bz (magnetic induction vector components) data of a round shaped magnet.

First, I arranged my data in a meshgrid and then I performed some matrix operations ("transformations") to rotate and shift the magnetic field data (which works well and really fast).

But nooooow, I need to evaluate the Bx, By, Bz component values at specific points in my global COS. As example: I grid my data with a meshsize of 0.5 mm. The problem is that after rotation about for example 30 ° about z-axis there is no grid point at [0.5 0.5 0.5] anymore and I can't evaluate the B field data there.

I can bring my transformed field data back to a grid using griddata and interpolate the B values on the 0.5 mm meshsize to evaluate at this point... But this is inconceivably slow......

Does anyone have a better way to interpolate my transformed magnetic field data back to the original meshgrid OR to do the mentioned evaluation on any other way?

Many thanks in advance! I spent my whole last 2 days trying things like scatteredInterpolant (which is also not really faster) and griddedInterpolant, which seems to be fast, but can't handle my rotated field data....

[darova]

If i understood you correctly you want to extrapolate your data

Is this correct?

[Oliver Vogt]

Not really... I have griddata in x-, y- and z- direction: [-7.5 : 7.5; -7.5 : 7.5; 0 : 6] and the corresponding Bx, By, Bz component values.

If I grid this data on a meshsize of 0.5 mm, I can evaluate the B components for example at this point: [0.5; 0.5; 0.5]

But if I rotate my data in the COS with matrix transformations then there's no gridpoint at [0.5; 0.5; 0.5] anymore.

I need to interpolate my field data on the desired point e.g. [0.5; 0.5; 0.5].

And I'm looking for a faster way than doing this with griddata!

thanks!

[darova]

Is it possible to make a simple drawing or sketch? I still don't get what is going on

[Oliver Vogt]

Yes. I'll try :) I will just do it in 2D, but its actually the same issue for my 3D data...

[Oliver Vogt]

You might also have a look on my comment down there.. I posted my files and code aswell.

And if you want, we could also have a 10 mins Skype or Zoom Conversation, which may lead to success faster... Also easier for you maybe... I can show you my project and you can interact easier...Afterwards (if we found a solution), I would post it here. But this is just an offer, we can also do it the usual way here in this forum :).

Again, many thanks in advance!! Have a nice day!

[darova]

I saw your codes. They are long, i don't want to spend a whole day to trying understand it

1 question: do you do you calculations on one matrix a few times? Or one matrix (pointcloud) - one set of calculations?

My point is: can you griddata once and then use interp2 ?

[Oliver Vogt]

Hello darova,

I actually arrange my data with griddata first. Then I do matrix operations (like rotation) in maybe 100 or more steps (to get a fine resolution of a rotation over time).

Then I need to evaluate the Bx, By, Bz data on specific points in the x-,y-z-COS (for example [0.5 0.5 0.5]) each step!

[Oliver Vogt]

I tried (in my 3D grid) interp3... It works when there is no rotation applied. The problem is that after rotation (when the grid is not in the original state anymore) I receive this error:

Error using interp3 (line 150)

Input grid is not a valid MESHGRID.

[Oliver Vogt]

btw the "rotation" is applied with a matrix transformation:

% Transformation matrix
R_z = [...
    cosd(phi) -sind(phi) 0    0; ...
    sind(phi) cosd(phi)  0    0; ...
    0         0          1    0; ...
    0         0          0    1];

% Coordinate and vector transformation
XYZ1 = [...
    GridMagnetData.X(:), ...
    GridMagnetData.Y(:), ...
    GridMagnetData.Z(:), ...
    ones(numel(GridMagnetData.X),1)];

rotXYZ1 = XYZ1 * R_z';

BxByBz1 = [...
    GridMagnetData.Bx(:), ...
    GridMagnetData.By(:), ...
    GridMagnetData.Bz(:), ...
    ones(numel(GridMagnetData.Bx),1)];

rotBxByBz1 = BxByBz1 * R_z';

% Reshape in meshgrid
RotatedMagnet.X = reshape(rotXYZ1(:,1), size(GridMagnetData.X));
RotatedMagnet.Y = reshape(rotXYZ1(:,2), size(GridMagnetData.Y));
RotatedMagnet.Z = reshape(rotXYZ1(:,3), size(GridMagnetData.Z));

RotatedMagnet.Bx = reshape(rotBxByBz1(:,1), ...
    size(GridMagnetData.Bx));
RotatedMagnet.By = reshape(rotBxByBz1(:,2), ...
    size(GridMagnetData.By));
RotatedMagnet.Bz = reshape(rotBxByBz1(:,3), ...
    size(GridMagnetData.Bz));

[darova]

What if transformate point location?

You are transformating data (mesh) all time. But what if transformate positions of points. You can then use your original mesh

[Oliver Vogt]

Sounds good at the first moment...! Then I would need a really fine grid at the beginning right? And can you give me some practical tips how to do that, especially in contrast to the way I did the transformation in my comment above with the Transformation matrix R_z?

Thanks! I think this could be the right way!

Answer 1

[Oliver Vogt]

Actually, what I need, becomes clear when you look at my file...

Just open up example.m and run it.

For example type "1" [mm] if it asks for meshsize und then choose the rawMagneticFieldData.fld file...

then its computing...

You can see what is done in general when you look at the plots.

Actually, you only have to look at the file "globalGridMagnetData_working.m", here I need a working and FAST alternative to the griddata commands as they become superslow if I compute more steps than in the example....(I need like > 100 steps and in the example it only computes five)

Answer 2

[darova]

Here is an example. It's just a surface in 3D. I rotated X and Y and found new Z

See script and try to calculate value of point with your method (griddata)

clc,clear
n = 20;
[X,Y] = meshgrid(linspace(0,5,n));
Z = 10-X.^2-10*Y.^2;        % some random data
a = 30;                     % angle of rotation
R = @(a)[cosd(a) sind(a);-sind(a) cosd(a)]; % rotation matrix
P = [4 1];                  % original point
P1 = R(-a)*P';              % rotate point (backward)
PP = [P(:) P1(:)]';
PP(:,3) = interp2(X,Y,Z,PP(:,1),PP(:,2));   % get Z data for points
V = R(a)*[X(:) Y(:)]';      % rotate data (forward)
X1 = reshape(V(1,:),[n n]); % restore size
Y1 = reshape(V(2,:),[n n]); % resotre size
Z1 = Z;                     % Z is the same (only X and Y changed)

cla
plot3(PP(:,1),PP(:,2),PP(:,3),'or');
h1 = surface(X,Y,Z,'facecolor','b');
h2 = surface(X1,Y1,Z1,'facecolor','r');

set([h1 h2],'edgecolor','none')
alpha(0.3)

axis vis3d 
view(2)

I don't understand why do you interpolate magnetic values here? I thought you calculating them using griddata?

BxByBz1 = [...
    GridMagnetData.Bx(:), ...
    GridMagnetData.By(:), ...
    GridMagnetData.Bz(:), ...
    ones(numel(GridMagnetData.Bx),1)];
rotBxByBz1 = BxByBz1 * R_z';

• Then I would need a really fine grid at the beginning right?

Not really, usual grid

[Oliver Vogt]

Uhm... that was also not 100 % what I'm looking for... I don't get why you do the backward rotation first???

And I also dont have a function for my values in each point but I have explicit values Bx, By, Bz for every x-,y-z- Coordinate combination...

To your point:

I don't understand why do you interpolate magnetic values here? I thought you calculating them using griddata?

BxByBz1 = [...
    GridMagnetData.Bx(:), ...
    GridMagnetData.By(:), ...
    GridMagnetData.Bz(:), ...
    ones(numel(GridMagnetData.Bx),1)];
rotBxByBz1 = BxByBz1 * R_z';

I actually just do the matrix operation here?

I think I'll ask this question again in the forum and instantly use my sketch.

[darova]

• I actually just do the matrix operation here?

You can't apply geometrical transformation with magnetic field. You transformate your XYZ matrix and then calculate magnetic field in new position using griddata or interp2(whatever)

[Oliver Vogt]

I think you didn't understand my intention at all. I want to AVOID the griddata function as it is super slow!!!

That's why I reopened the question in another thread.

[darova]

Can i close this question?

[Oliver Vogt]

yes