How can I create a 3D intensity map, or 3D scatter plot where colour represents intensity?

I am hoping to create a 3D plot for a data set (a 3D matrix) where each voxel has four pieces of data: an x position, y position, z position, and intensity. I would like to create a scatter plot (or other type of intensity map) where each point is appropriately located, and its colour relates its value. The result might look something like the example of 'Visualize function of 3 variables' on this page: https://www.mathworks.com/help/matlab/visualize/visualizing-four-dimensional-data.html. However, I have thus far been unable to replicate that result with my own data, which is organized in a 3D matrix.

For example, if I have this matrix M:

A = [1 2 3; 4 5 6; 7 8 9];
A(:,:,2) = [10 11 12; 13 14 15; 16 17 18];
B = cat(3,A,[3 2 1; 0 9 8; 5 3 7]);
M = B./max(B);

how do I create a 3D plot where the voxel with value 0.3333 is located at point 3, 2, 3 on the 3D plot and is blue to show that it has a relatively low value, while the points representing voxels with value 1.000 (the maximum) are red. I would also like to include a sidebar which shows the colour scale.

In case it is helpful, here is the code from the above-mentioned example (i.e., not my code):

cla
load accidents hwydata                             % load data

long = -hwydata(:,2);                              % longitude data
lat = hwydata(:,3);                                % latitude data
rural = 100 - hwydata(:,17);                       % percent rural data
fatalities = hwydata(:,11);                        % fatalities data

scatter3(long,lat,rural,40,fatalities,'filled')    % draw the scatter plot
ax = gca;
ax.XDir = 'reverse';
view(-31,14)
xlabel('W. Longitude')
ylabel('N. Latitude')
zlabel('% Rural Population')

cb = colorbar;                                     % create and label the colorbar
cb.Label.String = 'Fatalities per 100M vehicle-miles';

Thank you in advance for any help you can provide!

Answer 1

[Kevin Holly]

The dimensions of the matrix, M, is 3x3x3. I am going to assume the third dimension is your x, y, and z coordinates, where 1 is x, 2 is y, and 3 is z. With that established, what values do you want expressed in the 4th dimension (the colors).

Define variables

A = [1 2 3; 4 5 6; 7 8 9];
A(:,:,2) = [10 11 12; 13 14 15; 16 17 18];
B = cat(3,A,[3 2 1; 0 9 8; 5 3 7]);
M = B./max(B);

Create scatter plot and color bar

scatter3(M(:,:,1),M(:,:,2),M(:,:,3),40,'filled')
xlabel('x')
ylabel('y')
zlabel('z')
colorbar

[Joan B]

Thank you! This is exactly what I'm looking for. One further question, how do I adjust this if my matrix is not a cube (it's 61 x 61 x 101)?

[Kevin Holly]

Joan,

You need to input three different different vector arrays (one for each axis). The three vector arrays need to be the same size. I'm not sure how to get that from your 3D matrix as I don't know your data.

xarray = rand(1,61);
yarray = rand(1,61);
zarray = rand(1,61);
carray = rand(1,61); % An array with intensity values for colorbar
scatter3(xarray,yarray,zarray,40,carray,'filled')
colorbar
caxis([0.3,1])

[Joan B]

Hi Kevin,

Here is my code. Sorry that I did not include it initially. I thought it would be easier to include a shorter dummy version but I guess that backfired! I had section breaks in my code but I can't figure out how to include those here, sorry. 'tank' is what I want to plot.

% Equivalence of square and circular fields: R = 0.558 a
    % where R is the radius of the circle
    % and the side of the square has length a

R = 5; % radius of beam at surface of tank (cm)
a = R./0.558 % equivalent square side length (cm)

%%%%%%%%%%%%%%%%%%%%%%%%%

% get surface PDD
PDD_surface_mia_et_al = [0.3821, 0.4100, 0.4320, 0.4010]; % measured surface PDDs
field_size= [4, 6, 8, 10]; % square field side length
PDD_surface = interp1(field_size,PDD_surface_mia_et_al,a); % interpolate PDD at surface for field side length a

% This is the depth and PDD from surface up to 30 cm depth
% use steps of 0.5 cm since this is the voxel size
xmid = [0, 1.5:0.5:30]; % distance from surface of tank (cm)
fmid = [PDD_surface, 1, 0.988, 0.9692, 0.9505, 0.9295, 0.9084, 0.8874, 0.8664, 0.8454, 0.8244, 0.8041, 0.7838, 0.7641, 0.7443, 0.7251, 0.7058, 0.6875, 0.6693, 0.6515, 0.6338, 0.6170, 0.6002, 0.5847, 0.5692, 0.5537, 0.5382, 0.5240 ,0.5097, 0.4957, 0.4817, 0.4690, 0.4562, 0.4437, 0.4312, 0.4197, 0.4082, 0.3972, 0.3862, 0.3757, 0.3652, 0.3552, 0.3452, 0.3360, 0.3268, 0.3183, 0.3098, 0.3013, 0.2928, 0.2850, 0.2773, 0.2698, 0.2623, 0.2553, 0.2483, 0.2418, 0.2353, 0.2288, 0.2223];

% get PDD between 0 cm and 1.5 cm
pshallow = polyfit(xmid(1,1:2),fmid(1,1:2),1); % get fit equation
xshallow = 0:0.5:1.5; % depths in cm up to 1.5 cm
fshallow = polyval(pshallow,xshallow); % results of fit

% get extrapolated PDD for 30 to 50 cm depth
pdeep = polyfit(xmid(1,30:length(xmid)),fmid(1,30:length(fmid)),3); % get fit equation
xdeep = 30:0.5:50; % depth up to 50 cm
fdeep = polyval(pdeep,xdeep); % results of fit

% combine all of the PDDs
depth = 0:0.5:50; % depth in cm
PDD = [fshallow, fmid(1,3:length(fmid)), fdeep(1,2:length(fdeep))]; % PDD from 0 cm to 50 cm depth

%%%%%%%%%%%%%%%%%%%%%%%%%

% Create Gaussian distributions for lateral dose profiles
x = [-15:0.5:15]; % distances along x or y direction (cm), centred at 0 cm, 0.5cm voxel side length
lat_profile = normpdf(x,0,1); % Gaussian dist with mean = 0, sigma = 1 cm

% make the centre of the distribution match the dose on the CAX
% create profiles for 1.5 cm, 5 cm, 10 cm, 20 cm
x_lat_p_015 = lat_profile./max(lat_profile).*PDD(find(depth == 2)); % profile at 1.5 cm deep
x_lat_p_05 = lat_profile./max(lat_profile).*PDD(find(depth == 5)); % 5 cm deep
x_lat_p_10 = lat_profile./max(lat_profile).*PDD(find(depth == 10)); % 10 cm deep
x_lat_p_20 = lat_profile./max(lat_profile).*PDD(find(depth == 20)); % 20 cm deep

y_lat_p_015 = lat_profile./max(lat_profile).*PDD(find(depth == 2)); % profile at 1.5 cm deep
y_lat_p_05 = lat_profile./max(lat_profile).*PDD(find(depth == 5)); % 5 cm deep
y_lat_p_10 = lat_profile./max(lat_profile).*PDD(find(depth == 10)); % 10 cm deep
y_lat_p_20 = lat_profile./max(lat_profile).*PDD(find(depth == 20)); % 20 cm deep

%%%%%%%%%%%%%%%%%%%%%%%%%

% Put the lateral and CAX PDDs into one matrix
% first do the CAX
tank = zeros((30*2 + 1), (30*2 + 1), (50*2 + 1)); %pre-allocate 30cm x 30cm x 50cm matrix with 0.5 cm voxels
for k = 1:length(PDD)
    tank(k,31,51) = PDD(k);
end


% then do lateral

% 1.5 cm deep
tank = zeros((30*2 + 1), (30*2 + 1), (50*2 + 1));
for k = 1:length(x_lat_p_015)
    tank(31,k,find(depth == 1.5)) = x_lat_p_015(k);
end

for k = 1:length(x_lat_p_015)
    tank(k,31,find(depth == 1.5)) = x_lat_p_015(k);
end

% 5 cm deep
tank = zeros((30*2 + 1), (30*2 + 1), (50*2 + 1));
for k = 1:length(x_lat_p_015)
    tank(31,k,find(depth == 5)) = x_lat_p_05(k);
end

for k = 1:length(x_lat_p_015)
    tank(k,31,find(depth == 5)) = x_lat_p_05(k);
end

% 10 cm deep
tank = zeros((30*2 + 1), (30*2 + 1), (50*2 + 1));
for k = 1:length(x_lat_p_015)
    tank(31,k,find(depth == 10)) = x_lat_p_10(k);
end

for k = 1:length(x_lat_p_015)
    tank(k,31,find(depth == 10)) = x_lat_p_10(k);
end

% 20 cm deep
tank = zeros((30*2 + 1), (30*2 + 1), (50*2 + 1));
for k = 1:length(x_lat_p_015)
    tank(31,k,find(depth == 20)) = x_lat_p_20(k);
end

for k = 1:length(x_lat_p_015)
    tank(k,31,find(depth == 20)) = x_lat_p_20(k);
end

[Kevin Holly]

% Equivalence of square and circular fields: R = 0.558 a
% where R is the radius of the circle
% and the side of the square has length a
R = 5; % radius of beam at surface of tank (cm)
a = R./0.558; % equivalent square side length (cm)
%%%%%%%%%%%%%%%%%%%%%%%%%
% get surface PDD
PDD_surface_mia_et_al = [0.3821, 0.4100, 0.4320, 0.4010]; % measured surface PDDs
field_size= [4, 6, 8, 10]; % square field side length
PDD_surface = interp1(field_size,PDD_surface_mia_et_al,a); % interpolate PDD at surface for field side length a
% This is the depth and PDD from surface up to 30 cm depth
% use steps of 0.5 cm since this is the voxel size
xmid = [0, 1.5:0.5:30]; % distance from surface of tank (cm)
fmid = [PDD_surface, 1, 0.988, 0.9692, 0.9505, 0.9295, 0.9084, 0.8874, 0.8664, 0.8454, 0.8244, 0.8041, 0.7838, 0.7641, 0.7443, 0.7251, 0.7058, 0.6875, 0.6693, 0.6515, 0.6338, 0.6170, 0.6002, 0.5847, 0.5692, 0.5537, 0.5382, 0.5240 ,0.5097, 0.4957, 0.4817, 0.4690, 0.4562, 0.4437, 0.4312, 0.4197, 0.4082, 0.3972, 0.3862, 0.3757, 0.3652, 0.3552, 0.3452, 0.3360, 0.3268, 0.3183, 0.3098, 0.3013, 0.2928, 0.2850, 0.2773, 0.2698, 0.2623, 0.2553, 0.2483, 0.2418, 0.2353, 0.2288, 0.2223];
% get PDD between 0 cm and 1.5 cm
pshallow = polyfit(xmid(1,1:2),fmid(1,1:2),1); % get fit equation
xshallow = 0:0.5:1.5; % depths in cm up to 1.5 cm
fshallow = polyval(pshallow,xshallow); % results of fit
% get extrapolated PDD for 30 to 50 cm depth
pdeep = polyfit(xmid(1,30:length(xmid)),fmid(1,30:length(fmid)),3); % get fit equation
xdeep = 30:0.5:50; % depth up to 50 cm
fdeep = polyval(pdeep,xdeep); % results of fit
% combine all of the PDDs
depth = 0:0.5:50; % depth in cm
PDD = [fshallow, fmid(1,3:length(fmid)), fdeep(1,2:length(fdeep))]; % PDD from 0 cm to 50 cm depth
%%%%%%%%%%%%%%%%%%%%%%%%%
% Create Gaussian distributions for lateral dose profiles
x = [-15:0.5:15]; % distances along x or y direction (cm), centred at 0 cm, 0.5cm voxel side length
lat_profile = normpdf(x,0,1); % Gaussian dist with mean = 0, sigma = 1 cm
% make the centre of the distribution match the dose on the CAX
% create profiles for 1.5 cm, 5 cm, 10 cm, 20 cm
x_lat_p_015 = lat_profile./max(lat_profile).*PDD(find(depth == 2)); % profile at 1.5 cm deep
x_lat_p_05 = lat_profile./max(lat_profile).*PDD(find(depth == 5)); % 5 cm deep
x_lat_p_10 = lat_profile./max(lat_profile).*PDD(find(depth == 10)); % 10 cm deep
x_lat_p_20 = lat_profile./max(lat_profile).*PDD(find(depth == 20)); % 20 cm deep
y_lat_p_015 = lat_profile./max(lat_profile).*PDD(find(depth == 2)); % profile at 1.5 cm deep
y_lat_p_05 = lat_profile./max(lat_profile).*PDD(find(depth == 5)); % 5 cm deep
y_lat_p_10 = lat_profile./max(lat_profile).*PDD(find(depth == 10)); % 10 cm deep
y_lat_p_20 = lat_profile./max(lat_profile).*PDD(find(depth == 20)); % 20 cm deep
%%%%%%%%%%%%%%%%%%%%%%%%%

% Put the lateral and CAX PDDs into one matrix
% first do the CAX
tank = zeros((30*2 + 1), (30*2 + 1), (50*2 + 1)); %pre-allocate 30cm x 30cm x 50cm matrix with 0.5 cm voxels
for k = 1:length(PDD)
    tank(k,31,51) = PDD(k);
end
% then do lateral
% 1.5 cm deep
tank = zeros((30*2 + 1), (30*2 + 1), (50*2 + 1));
for k = 1:length(x_lat_p_015)
    tank(31,k,find(depth == 1.5)) = x_lat_p_015(k);
end
for k = 1:length(x_lat_p_015)
    tank(k,31,find(depth == 1.5)) = x_lat_p_015(k);
end
% 5 cm deep
tank = zeros((30*2 + 1), (30*2 + 1), (50*2 + 1));
for k = 1:length(x_lat_p_015)
    tank(31,k,find(depth == 5)) = x_lat_p_05(k);
end
for k = 1:length(x_lat_p_015)
    tank(k,31,find(depth == 5)) = x_lat_p_05(k);
end
% 10 cm deep
tank = zeros((30*2 + 1), (30*2 + 1), (50*2 + 1));
for k = 1:length(x_lat_p_015)
    tank(31,k,find(depth == 10)) = x_lat_p_10(k);
end
for k = 1:length(x_lat_p_015)
    tank(k,31,find(depth == 10)) = x_lat_p_10(k);
end
% 20 cm deep
tank = zeros((30*2 + 1), (30*2 + 1), (50*2 + 1));
for k = 1:length(x_lat_p_015)
    tank(31,k,find(depth == 20)) = x_lat_p_20(k);
end
for k = 1:length(x_lat_p_015)
    tank(k,31,find(depth == 20)) = x_lat_p_20(k);
end

Create vector coordinates for 30 cm x 30 cm x 50 cm tank

x_coord = 1:0.5:30+1;
y_coord = 1:0.5:30+1;
z_coord = 1:0.5:50+1;

Create vectors (x, y, and z) that represent all points within 3D grid

[X,Y,Z] = meshgrid(x_coord,y_coord,z_coord);
x=reshape(X,1,[]);
y=reshape(Y,1,[]);
z=reshape(Z,1,[]);

Reshape 3D tank matrix into a single array

tank_vector = reshape(tank,1,[]);

Create scatterplot

figure; scatter3(x,y,z,[],tank_vector,'filled')

% Remove Scatter points with value of zero
M = [x;y;z;tank_vector];
M(:,M(4,:)==0)=[];
figure; scatter3(M(1,:),M(2,:),M(3,:),[],M(4,:),'filled')
colorbar

[Joan B]

Amazing! Thank you so much!