How can I optimize this 4D contour plot?

I am trying to make a 4d contour plot using an array with ua, us, g, delta data plotted in {ua, us, g} 3D space with each point having a value of delta. The 4th value of delta will correspond to a color. The data file is attached as output.txt

There seem to be numerous similar questions to this one, most with answers that suggest using a meshgrid + isosurface + patch approach. I don't think this is applicable for me as I do not have a continuous surface over the 3D space. Rather, I have discrete data that may or may not be linearly spaced in each dimension. In this example, ua locations are exponentially spaced while us and g are linearly spaced.

Here is what I have so far:

name = sprintf('output.txt');
data = readtable(name);
data_arr = table2array(data);

ua = data_arr(:,1);
us = data_arr(:,2);
g = data_arr(:,3);
delta = data_arr(:,4);
    
figure
h = scatter3(ua, us, g, 12, delta,'MarkerFaceColor', 'Flat','MarkerFaceAlpha',.1,'MarkerEdgeAlpha',0.05);
xlabel('{\mu}_a (mm^{-1})');
ylabel('{\mu}_s (mm^{-1})');
zlabel('g');
grid on
colormap(flipud(turbo))
c=colorbar();
c.Label.String = 'delta';
caxis([0 0.6])

which results in the attached figure. I also pair this with the CaptureFigVid function (https://www.mathworks.com/matlabcentral/fileexchange/41093-create-video-of-rotating-3d-plot) to have a nice rotating video. The idea is that with enough sheets in the ua dimension, I can replicate a continuous 3D contour map. I'd like to optimize this method to create a figure with easily visible color contour in the 3D space. Perhaps someone has some experience with this or tips I can try? Perhaps even there is a better, more sophisticated method.

So far, I've tried to play with the values of

'MarkerFaceAlpha',.1,'MarkerEdgeAlpha',0.05

to make the volume transparent.

Perhaps also I could animate alternating through slices of constant ua, like the last example here: http://web.mit.edu/8.13/matlab/MatlabTraining_IAP_2012/AGV/DemoFiles/ScriptFiles/html/Part7_SlicesIsosurfaces.html. However, it looks like the slice() function is not compatible with data in this style.

Rather than slice, I can use

xslice = 0.230849; %% for instance.
scatter3(xslice,us(ua==xslice),g(ua==xslice),ones(size(us(ua==xslice)))*50,delta(ua==xslice))

and write a for loop to step through each value of ua.

Alternatively, I've added a loop to show the slices for each unique value of ua.

unique_ua=unique(ua);
for i = 1:length(unique_ua)
 
    ua_scatter = repmat(unique_ua(i),M,1);
    scatter3(ua_scatter,us(ua==unique_ua(i)),g(ua==unique_ua(i)),ones(size(us(ua==unique_ua(i))))*50,delta(ua==unique_ua(i)))
    xlim([0 0.675]);
    ylim([0.75 14]);
    zlim([0 0.875]);
    xlabel('{\mu}_a (mm^{-1})');
    ylabel('{\mu}_s (mm^{-1})');
    zlabel('g');
    grid on
    colormap(flipud(turbo))
    c=colorbar();
    c.Label.String = 'delta';
    caxis([0 0.6])

    pause(0.01);
end

Any other ideas or tips? Thanks.

Answer 1

[Ran Yang]

It looks like you're trying to visualize volumetric data. One way to do that is with slice, but you need to first restructure your data_arr from an Nx4 array of coordinates and values to a 3D volumetric array. V is also more space efficient since you're not repeating the same coordinate information n = 55 times.

data = readtable('output.txt');

% some overhead
x = data{:, 1}; y = data{:, 2}; z = data{:, 3}; c = data{:, 4};
ux = unique(x); uy = unique(y); uz = unique(z);

% assuming you sampled each data point (x,y,z) exactly once,
% brute force convert data_arr to volumetric array, v
v = zeros(length(ux), length(uy), length(uz));
for i = 1:size(v, 1)
    for j = 1:size(v, 2)
        for k = 1:size(v, 3)
            v(i, j, k) = c(x == ux(i) & y == uy(j) & z == uz(k));
        end
    end
end

% visualize
figure, slice(ux, uy, uz, v, median(ux), median(uy), median(uz));

From this point you can do whatever you need to prettify your data (such as adjusting colormap, rescaling v, etc.), and it should be compatible with other functions (such as isosurface and contourslice). If you want to visualize an evenly spaced grid, see interp3.

[iontrap]

Thanks Ran,

I use the slightly modified version of your code:

data = readtable('output.txt');
x = data{:, 1}; y = data{:, 2}; z = data{:, 3}; c = data{:, 4};
ux = unique(x); uy = unique(y); uz = unique(z);
v = zeros(length(ux), length(uy), length(uz));
for i = 1:size(v, 1)
    for j = 1:size(v, 2)
        for k = 1:size(v, 3)
            v(i, j, k) = c(x == ux(i) & y == uy(j) & z == uz(k));
        end
    end
end

scatter3(x(:), y(:), z(:), 30, v(:), 'filled');
colormap(flipud(turbo))
c=colorbar();
c.Label.String = 'delta';
caxis([0 0.6])

however it looks like the v data no longer corresponds to the previous plot of "delta" - see attached picture, Yang_pic. The blue area of v should begin around z = ~0.45, like the figure attached as my_pic. The low delta region corresponds well in the x and y dimensions, but not z.

[Ran Yang]

For some reason that's not apparent to me at this moment, my code swaps the x and z dimensions. To fix this issue, use:

v2 = permute(v, [3 2 1]); % swap x and z

...or in the original code, swap the output order of v:

% reorder the allocation of v as well, for robustness
v = zeros(length(uz), length(uy), length(ux));
for i = 1:size(v, 1)
    for j = 1:size(v, 2)
        for k = 1:size(v, 3)
            % switch i and k (but why...)
            v(k, j, i) = c(x == ux(i) & y == uy(j) & z == uz(k));
        end
    end
end

Either way, you should now be able to run

scatter3(x(:), y(:), z(:), 30, v2(:), 'filled')
colormap(flipud(turbo)), caxis([0 0.6])

...and reproduce your result.

Note that you need to be careful plotting scatter3(x, y, z). It only works if data is row-sorted (which it just happened to be in this case). If it isn't row sorted, you will likely need to sort it:

data_sorted = sortrows(data, [1 2]);

[iontrap]

Thanks Ran, this fixed the issue I was having. Another question I have: this method of brute forcing the c values into a volume v. How can I ensure that this method is robust in the cases of unequal lengths of ux, uy, and uz?

For instance, if

length(uz)==length(uy)!=length(ux)

can we use a variation of this:

v = zeros(length(uz), length(uy), length(ux));
for i = 1:size(v, 1)
    for j = 1:size(v, 2)
        for k = 1:size(v, 3)
            v(k, j, i) = c(x == ux(i) & y == uy(j) & z == uz(k));
        end
    end
end

where we change the iteration limit for each case?

[Ran Yang]

That code should already be robust to unequal lengths of ux,.... The size of v and the number of iterations in each for loop is determined by how many unique points you sampled per dimension (that's what length(ux), ... is).

As an aside, I think a much faster way of forcing a volumetric array (assuming that you have exactly 1 data point per x,y,z) would be something like:

data_sorted = sortrows(data, 1:3);
v = reshape(data_sorted{:, 4}, length(ux), length(uy), length(uz));

Don't quote me on that but you're free to try it.

[iontrap]

Hi Ran, when I try the original "brute force" code on the case where ux is a 65x1 array while uy and uz are 55x1, I get the error that

"Index exceeds the number of array elements (55)."

When I run the reshape method, I have the error,

"Number of elements must not change. Use [] as one of the size inputs to automatically calculate the appropriate size for that dimension."

Also attached is the data file I use for testing.

[Ran Yang]

My mistake: the indexing error is solved by switching size(v, 1) and size(v, 3) in the nested for loops. When you see that kind of error, it's often because the for loop indices are mismatched with the data you're trying to loop it over (in this case, I accidentally told it to loop over the number of elements in uz while using the data in x and vice versa).

There is no issue with reshape. You're telling me that ux is 65x1, but your attached data is sampled 80x55x55 times. Perhaps you are mixing up your variables?

[iontrap]

Yes, I should've used a more general example than 65x1 for ux. Thanks for all of your great help.

Also, in case you were wondering, the reshape method takes ~ 2 seconds to plot data on my machine, while the loop method takes about 60 seconds.

[iontrap]

Hm, now when I try the code:

data = readtable('output_blood_2.txt');
x = data{:, 1}; y = data{:, 2}; z = data{:, 3}; c = data{:, 4};
ux = unique(x); uy = unique(y); uz = unique(z);

data_sorted = sortrows(data, 1:3);
v = reshape(data_sorted{:, 4}, length(ux), length(uy), length(uz));

figure, slice(ux, uy, uz, v, median(ux), median(uy), median(uz));

I get the error:

"GridVectors must define a grid whose size is compatible with the Values array."

[Ran Yang]

This whole x/z axis swapping thing is getting annoying, isn't it?

v = reshape(data_sorted{:, 4}, length(uz), length(uy), length(ux)); % swap x/z
v(v > 15) = 15; % for better visualization of output_blood_2.txt only
figure, slice(uz, uy, ux, v, median(uz), median(uy), median(ux)); % swap x/z

For future reference, that GridVectors error just means that your ux,uy,uz "rescaling" vectors (the first 3 input arguments into slice) don't agree with the dimensions of v. You should always able to see your entire volume just by omitting those 3 scaling arguments and directly calling:

figure, slice(v, ...)

If even that doesn't work, or produces resulting patterns that are different than what you expect, then you'd have a much bigger problem (such as reshape not behaving as expected or your data are not 1:1).

Side note: If you have the time, I would recommend you try to figure out what's up with the whole x/z axis swap. I don't really get why I can't just assign it the proper way (why reshape has to be done using z,y,x instead of x,y,z, nested for loop has to be done z,y,x instead of x,y,z, etc.) Maybe that could even be posted as a separate Question on its own, haha.

[iontrap]

I will look into the x/z axis issue. My knowledge of Matlab is obviously a little lacking, but I'll try to figure it out. Anyways, I can now use isosurface and similar functions to produce very nice looking figures - thanks so much for all of your help.

I think this question and answer may also be useful for others as I saw many unanswered similar versions of this question.

Thanks,

Zachary