Isosurface function, irregular coordinates (Input grid is not a valid MESHGRID)

I have a data file from a physics model that shows the density of a certain material around a object. The data file is in spherical coordinates and is reshaped to its original dimensions (Theta,Phi,R and density). Sph2cart is used to transform from spherical coordinates to cartesian (x,y,z) coordinates. To visualize the different densities around the object I found the isosurface function in MATLAB (example: http://www.mathworks.com/matlabcentral/answers/110977-3d-density-plot-multiple-isosurfaces-on-the-same-plot).

This is the code that I currently have:

input = importdata( 'denNa20.dat' );
input(1,:)=[]; 

theta = reshape(input(:,3),200,20,40);
phi = reshape(pi/2 - input(:,2),200,20,40);
r = reshape(input(:,1),200,20,40);
density = reshape(input(:,4),200,20,40);

[x,y,z] = sph2cart(theta,phi,r);

% This has ofcourse the complete wrong dimensions but then it works
% [x,y,z] = meshgrid(1:1:20,1:1:200,1:1:40);

p = patch(isosurface(y,x,z,density,0.00001));
isonormals(x,y,z,density,p);
set(p,'FaceColor','red','EdgeColor','none','FaceAlpha',0.15);
daspect([1 1 1])
view(3)
grid on
camlight; lighting phong

When I run the code I receive the following errors:

Error using interp3 (line 146) Input grid is not a valid MESHGRID.

Error in isonormals (line 75) n(:,1)=interp3(x, y, z, nx, verts(:,1), verts(:,2), verts(:,3));

Error in data_import_plot (line 16) isonormals(x,y,z,density,p);

If I create my own meshgrid with very easy x,y,z coordinates (check %comment in code) it works. I don't know what I am doing wrong. I hope some one can help me.

Cheers, Jeroen

P.S. if you want you can download the data file from this link https://www.dropbox.com/s/msphgmg2oyi91cx/denNa20.dat?dl=0

Answer 1

[Jeroen Terwisscha van Scheltinga]

I found out that the problem was that isosurface does not accept irregular x,y and z coordinates. This is what I found to be the solution. I used scatteredInterpolant to interpolate the density and then created my own meshgrid.

input = importdata( 'denNa20.dat' );
input(1,:)=[]; 

[x,y,z] = sph2cart(input(:,3),pi/2 - input(:,2),input(:,1));
density = input(:,4);

F = scatteredInterpolant(x,y,z,density);
xRange = linspace(min(x),max(x),75);
yRange = linspace(min(y),max(y),75);
zRange = linspace(min(z),max(z),75);
[x,y,z] = meshgrid(xRange,yRange,zRange);
density = F(x,y,z);

axis([-5000,5000,-5000,5000,-5000,5000])

p = patch(isosurface(x,y,z,density,0.00001));
isonormals(x,y,z,density,p);
set(p,'FaceColor','red','EdgeColor','none','FaceAlpha',0.50);
daspect([1 1 1])
view(3)
grid on
camlight; lighting phong