Having Trouble Plotting 4D data

I'm a new user to Matlab so I'm hoping someone can help.

I have the data output from an atmospheric modeling program I am developing that I import into Matlab. The data is in the format x,y,z,value. The values represent pollution concentration in the atmosphere and the x,y,z are kilometer coordinates in the 3D region. I have tried the following code, but my plot ends up empty (but I do not get any error or diagnostic messages from the isosurface function). I got the mat & reshape ideas from another answer as a way to get my data into 3D arrays for the isosurface function:

EDU>> mat = [x(:) y(:) z(:) q_new(:)];
EDU>> mat = sortrows(mat,[3,1,2]);
EDU>> x = reshape(mat(:,1), [61 61 20]);
EDU>> y = reshape(mat(:,2), [61 61 20]);
EDU>> z = reshape(mat(:,3), [61 61 20]);
EDU>> q_new = reshape(mat(:,4), [61 61 20]);
EDU>> isovalue = 0.1671

isovalue =

      0.1671

EDU>> isosurface(x,y,z,q_new,isovalue)

My values are quite small but I don't think that should make a difference, so I'm wondering if my data is too scattered or sparse and do I need to pre-process it in some other way to make it "plot-able"?

Any advice would be greatly appreciated. Thanks!

[Andrew Newell]

Could you please format the code? One line per command, double space before the command, remove EDU>>.

[Jessica]

Sorry about the formatting. I just cut and pasted and didn't realize it was joined up.

EDU>> mat = [x(:) y(:) z(:) q_new(:)];

EDU>> mat = sortrows(mat,[3,1,2]);

EDU>> x = reshape(mat(:,1), [61 61 20]);

EDU>> y = reshape(mat(:,2), [61 61 20]);

EDU>> z = reshape(mat(:,3), [61 61 20]);

EDU>> q_new = reshape(mat(:,4), [61 61 20]);

EDU>> isovalue = 0.1671

isovalue =

0.1671

EDU>> isosurface(x,y,z,q_new,isovalue)

Answer 1

[Andrew Newell]

I don't see any obvious problem with the code, but of course I don't have the data. Maybe 0.1671 is outside the data range? You could try these commands to get some idea of how the data are distributed:

max(q_new(:))
min(q_new(:))
hist(q_new(:))

The isovalue represents a value of q_new, and isosurface tries to plot a surface separating values below isovalue and above - like contour, but for surfaces. If the points are scattered, it may be hard to define where that surface is. If your data are scattered (x,y, and z don't form a regular grid) TriScatteredInterp might help by moving the points to a regular grid. You might even want to apply smooth3 to the output.

[Matt Tearle]

Not sure if I can add much to Andrew's explanation, but think of a weather map. Isobars are pressure contours. i.e. you measure pressure at a bunch of points at a given altitude, then join up the points with the same reading for pressure. Now imagine doing that at a bunch of altitudes. The isobars should morph slightly as you go from one height to another, so if you stacked all the weather maps on top of each other, the isobars would form a kind of cross-section/wireframe of a surface. (That's why I like to look at contours in slices.) In your case the surfaces would represent all the locations in space where the pollutant concentration was the given value (0.1671 or whatever).

As you can imagine, it's extremely difficult to make these surfaces by joining up all the same values if the data is sparse.

[Jessica]

x, y, and z are a regular grid. The pollution concentration values I want to plot are varied and sparsely populate the grid at times (many points are 0 at times). I've tried TriScatteredInterp and now get something on the plot. Now I'm trying to determine if it is what I expect based on the data input. Thanks for the definition of isovalue. I couldn't find one anywhere in the documentation or the help.

Answer 2

[Matt Tearle]

I concur with Andrew about the code. Another diagnostic you could use is to visualize a series of slices:

figure
for k = 1:9
  subplot(3,3,k)
  [~,h] = contour(x(:,:,2*k),y(:,:,2*k),q_new(:,:,2*k));
  set(h,'ShowText','on')
  title(['z = ',num2str(z(1,1,2*k))])
end

Answer 3

[Jessica]

0.1671 is the maximum value for q_new so it should be in range. I'll try the slices next.

[Andrew Newell]

If it is the maximum value, you'll probably get a single point in your plot. Try a smaller value.

[Jessica]

I tried an isovalue between the min and the max values for the data and still no visible plot.

[Andrew Newell]

Without seeing your data, it's hard to know what is going wrong. Try generating your inputs using the commands below. What do you see when you run your code on this data?

x = (1:61)/61; y = x; z = (1:20)/20;

[x,y,z] = meshgrid(x,y,z);

q_new = exp(-x-y-z)*0.1815;

[Jessica]

With an isovalue of the maximum for q_new (0.1671), the plot appears to be empty. With an isovalue of (0.05) there is something plotted.

[Andrew Newell]

O.k., so the problem has to do with your data. Have you tried several different isovalues?

[Matt Tearle]

Have you tried looking at the slices? It might give you an idea of isovalues that might work. Maybe your data is too sparse for the values you're trying?

[Jessica]

I have tried several different isovalues but no luck. My data is definitely sparse. I haven't been able to find any guidance on choosing an isovalue or a good description of what it is meant to represent. I think that is why I am having trouble selecting the correct value.

[Jessica]

I have also looked at a slice and it did show values in the plot

[Andrew Newell]

I think we need to see the actual data. We don't have a good way of doing this on Matlab Answers. You could email the data to me (click on my name to see the address) and probably @Matt too (but I can't speak for him).

[Matt Tearle]

Fine by me. The reason I suggested slices is because you're more likely to see the contours in each slice, and you can eyeball how the slices would join from one to another. From that you can get an idea of what the isosurfaces should look like... or if they're not going to look like anything!

[Andrew Newell]

I totally agree.

[Jessica]

Thanks for all your help. One last question before I send the data. I used the function TriScatteredInterp and I am now seeing a plot. Does that mean that my data was just too scattered to plot previously? Also, I was wondering what the isovalue variable is supposed to represent?

[Andrew Newell]

Yes, I think it would help. For clarity, I have added to my answer above.

Answer 4

[Matt Tearle]

Here's another possibility for looking at the sparsity of the data. You said many of the values were 0, so maybe try this:

idx = mat(:,4)>0;
figure
text(mat(idx,1),mat(idx,2),mat(idx,3),num2str(mat(idx,4)))

This will be hideously busy if there's a lot of nonzero data, but if not, it might give you a good idea of how spread the values are.

Similarly, you could also try plotting points where the pollutant concentration falls within a certain range. Something like

figure
idx = (mat(:,4)>0.1) & (mat(:,4)<0.12);
plot3(mat(idx,1),mat(idx,2),mat(idx,3),'o')

[Jessica]

this is great! while my data does not seem to display very well as an isosurface, when I plot it as points with the colors varying by value, it shows what I am trying to convey.

Thanks for all your help!

Answer 5

[Saeed68gm]

Hello Jessica,

I am also new to MATLAB and I have a dataset that is similar to yours((x,y,z,value) for each point). The dataset is from a 3d scanner and it I want to show the data as a surface. for example this data is from a scanned gun, and I want to plot the gun surface( i figured using patch might be useful).

but i don't know how to use isosurface in this case. can you put a more sophisticated version of your code for me? I think you may have the solution to my problem...:)

[Walter Roberson]

Please start a new question for this, and in that question specify whether your scan is over a regular grid or is scattered.

[Saeed68gm]

Ok, I already posted a new topic called Problem plotting 4D data.

Hope someone can help:)