33 views (last 30 days)
Show older comments
Joan B on 28 Mar 2022
Commented: Joan B on 2 Apr 2022
Accepted Answer: Kevin Holly
Open in MATLAB Online
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!
0 Comments Show -2 older commentsHide -2 older comments
Show -2 older commentsHide -2 older comments
Sign in to comment.
Sign in to answer this question.
Accepted Answer
Kevin Holly on 28 Mar 2022
Edited: Kevin Holly on 28 Mar 2022
Open in MATLAB Online
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
5 Comments Show 3 older commentsHide 3 older comments
Show 3 older commentsHide 3 older comments
Joan B on 29 Mar 2022
Direct link to this comment
https://www.mathworks.com/matlabcentral/answers/1682834-how-can-i-create-a-3d-intensity-map-or-3d-scatter-plot-where-colour-represents-intensity#comment_2070704
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 on 29 Mar 2022
Direct link to this comment
https://www.mathworks.com/matlabcentral/answers/1682834-how-can-i-create-a-3d-intensity-map-or-3d-scatter-plot-where-colour-represents-intensity#comment_2071579
Edited: Kevin Holly on 29 Mar 2022
Open in MATLAB Online
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 on 30 Mar 2022
Direct link to this comment
https://www.mathworks.com/matlabcentral/answers/1682834-how-can-i-create-a-3d-intensity-map-or-3d-scatter-plot-where-colour-represents-intensity#comment_2075004
Open in MATLAB Online
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 on 31 Mar 2022
Direct link to this comment
https://www.mathworks.com/matlabcentral/answers/1682834-how-can-i-create-a-3d-intensity-map-or-3d-scatter-plot-where-colour-represents-intensity#comment_2076574
Open in MATLAB Online
% 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 on 2 Apr 2022
Direct link to this comment
https://www.mathworks.com/matlabcentral/answers/1682834-how-can-i-create-a-3d-intensity-map-or-3d-scatter-plot-where-colour-represents-intensity#comment_2079234
Amazing! Thank you so much!
Sign in to comment.
More Answers (0)
Sign in to answer this question.
See Also
Categories
MATLABGraphics2-D and 3-D PlotsDiscrete Data Plots
Find more on Discrete Data Plots in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- Deutsch
- English
- Français
- United Kingdom(English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)
Contact your local office
