Slicing 3D model data without using griddata

I'm trying to plot slices through quite large data sets (about 40K points at the moment, but probably at least 4 times that size in the future) which I receive in a specific output from my solver.

The data is output for each single timestep is a matrix (which I'm reading in from a txt file) like so:

NodeID NodePosX NodePosY NodePosZ VarA VarB VarC .....
1        0.2      0.1     0.1       1    1   2 
2        0.3      0.1     0.3       4    0   2

where the can be many variables. The data is sorted by NodeID and not the x,y or z position.

What I want to be able to do is use subplot to display a few slices through the volume for any given variable for each timestep so that I have a full set of images to make a video with eventually.

Currently my approach is painfully slow:

• use meshgrid to define a even grid in all three directions based on the max and min values of the node positions
• use griddata to generate the volume data for the meshed grid for each variable I'm interested in
• use slice on a few planes in a subplot for each variable

This currently works....but it's terribly slow. Currently I have data for 250 timesteps and I tried running my script on it and it takes about 3-4 minutes per timestep, per variable. So I want to plot the results for 3 variables (say components of velocity) I can generate only 6 timesteps worth of data in an hour.

My bottleneck is having to go from the vector node data to a volume data using griddata, but I can't think of a better way of doing this for the format the data I receive is in. Is there a way to plot the data with scatter3 and then interpolate slices through that?

Any ideas on how to do this more efficiently would be appreciated. I can't think of a quicker way at the moment.

Answer 1

[Mike Garrity]

You're inserting your data points into a 3D grid and then pulling a 2D slice out of that grid. That could make sense if you want to get a lot (most) of the slices from that 3D volume.

But if you're only getting a few slices, it would probably make more sense to take a subset of your data points which are in (near) your slice, and insert those into a 2D grid. Does that make sense?

[Mike Garrity]

Here's an example showing the two approaches:

%%Create example data
npoints = 100000;
x = randi(100,[npoints 1]);
y = randi(100,[npoints 1]);
z = randi(100,[npoints 1]);
v = randn(npoints,1);
data = [x,y,z,v];

%%Insert into 3D grid
[xq,yq,zq] = meshgrid(1:100,1:100,1:100);
vq = griddata(data(:,1),data(:,2),data(:,3),data(:,4),xq,yq,zq);

%%Slice 3D -> 2D and display
figure
xs = xq(:,:,23);
ys = yq(:,:,23);
vs = vq(:,:,23);
surf(xs,ys,vs)

%%Insert directly into 2D grid
figure
mask = find(z==23);
[x2,y2] = meshgrid(1:100,1:100);
v2 = griddata(data(mask,1),data(mask,2),data(mask,4),x2,y2);
surf(x2,y2,v2)

There are some differences in areas where there are no data points at Z==23. If this is a problem, you could grid to a 3D grid that only had a few slices in Z and then pull out the middle slice. You really only need 3 slices for the common interpolation schemes in griddata.

[John]

Thanks for the suggestion - I've just been checking whether or not it's suitable. I've just modified your code to highlight the problem of the 2D grid method. Once you plot the interpolated grid values there's quite a difference in the result in terms of both trend and value, which is just too large a variation when looking at such model results.

%%Create example data
npoints = 100000;
x = randi(100,[npoints 1]);
y = randi(100,[npoints 1]);
z = randi(100,[npoints 1]);
v = randn(npoints,1);
data = [x,y,z,v];
maxlimit = 1.1 * max(v);
minlimit = 0.9 * min(v);

%%Insert into 3D grid
[xq,yq,zq] = meshgrid(1:100,1:100,1:100);
vq = griddata(data(:,1),data(:,2),data(:,3),data(:,4),xq,yq,zq);

%%Define Slices
%Slice on X and Y axis and one at base of cell

xslice = [0]; 
yslice = 0; 
zslice = [23]; 

figure
h1 = slice(xq,yq,zq,vq,xslice,yslice,zslice,'linear');
set(h1,'FaceColor','interp',...
  'EdgeColor','none',...
  'DiffuseStrength',.8)
colorbar EastOutside
view(0,090)
% daspect([1,1,1])
axis tight
box on
colormap (jet(64))
set(gcf,'Renderer','zbuffer')
caxis([minlimit maxlimit])

%%Slice 3D -> 2D and display
figure
xs = xq(:,:,23);
ys = yq(:,:,23);
vs = vq(:,:,23);
surf(xs,ys,vs)
    lighting phong
    shading interp
    colorbar EastOutside
view(0,90)
colormap (jet(64))
set(gcf,'Renderer','zbuffer')
caxis([minlimit maxlimit])

%%Insert directly into 2D grid
figure
mask = find(z==23);
[x2,y2] = meshgrid(1:100,1:100);
v2 = griddata(data(mask,1),data(mask,2),data(mask,4),x2,y2);
surf(x2,y2,v2)
    lighting phong
    shading interp
    colorbar EastOutside
view(0,90)
colormap (jet(64))
set(gcf,'Renderer','zbuffer')
caxis([minlimit maxlimit])

I'm only changing the visualisation here but have run lots of random example data sets and the differences is always quite noticeable, which will affect the correctness of the results I'm plotting.

Any thoughts on how to improve the match between the griddata method and the masked method? maybe increase the mask to z=22,23,24 and do some interpolation?

[Mike Garrity]

>> maybe increase the mask to z=22,23,24 and do some interpolation?

Yes, that's my guess. In areas where now datapoints fall in the slice, the 3D version is able to interpolate from the neighboring slices. But I don't think that you'll need a lot of slices to get the same effect.

I think this picture suggests that:

hold on
scatter(data(mask,1),data(mask,2), ...
        [],data(mask,4),'filled', ...
        'MarkerEdgeColor','black')
set(gca,'SortMethod','childorder')