4D spherical shell heat map

15 次查看(过去 30 天)
I have written a code (attached) to plot particles striking a sphere of radius r - I can produce the image in example_1. I divide the shell into segments and calculate the number of particles falling within each.
I would like to represent the data as a smooth 3d heatmap over the spherical shell. I can produce the image shown in example_2 using scatter3 and a color bar, but I would like a smooth surface. I have tried to follow the example of my previous question using an isosurface -
however ux != uy != uz and I have the "Number of elements must not change" error.
I have also tried meshgrid(), but the size of photons results in an array that excedes maximum size preference.
Is there an easier way to smooth this data?
photons = readmatrix('input.txt');
r = 70; % in unit mm
prop = zeros(length(photons),3);
for i = 1:length(photons)
c(i) = ((photons(i,1)*photons(i,1)) + (photons(i,2)*photons(i,2)) + (photons(i,3)*photons(i,3)) - (r*r));
b(i) = 2*((photons(i,1)*photons(i,7)) + (photons(i,2)*photons(i,8)) + (photons(i,3)*photons(i,9)));
a(i) = ((photons(i,7)*photons(i,7)) + (photons(i,8)*photons(i,8)) + (photons(i,9)*photons(i,9)));
t(i) = (-b(i) + sqrt(b(i).^2 - 4*a(i)*c(i)))/(2*a(i));
prop(i,1) = photons(i,1) + t(i)*photons(i,7);
prop(i,2) = photons(i,2) + t(i)*photons(i,8);
prop(i,3) = photons(i,3) + t(i)*photons(i,9);
prop_sph(i,3) = sqrt(prop(i,1).^2 + prop(i,2).^2 + prop(i,3).^2); % calculate spherical radius (should be r)
prop_sph(i,1) = acos((prop(i,3))/(prop_sph(i,3))); % calculate spherical theta (elevation)
if prop(i,1) > 0 && prop(i,2) >= 0
prop_sph(i,2) = atan(prop(i,2)/prop(i,1)); % calculate spherical phi (azimuth)
end
if prop(i,1) > 0 && prop(i,2) < 0
prop_sph(i,2) = atan(prop(i,2)/prop(i,1)) + 2*pi;
end
if prop(i,1) < 0 && prop(i,2) >=0
prop_sph(i,2) = atan(prop(i,2)/prop(i,1)) + pi;
end
if prop(i,1) < 0 && prop(i,2) >=0
prop_sph(i,2) = atan(prop(i,2)/prop(i,1)) + pi;
end
if prop(i,1) < 0 && prop(i,2) < 0
prop_sph(i,2) = (atan(prop(i,2)/prop(i,1)) - pi) + 2*pi ;
end
if prop(i,1) == 0 && prop(i,2) > 0
prop_sph(i,2) = pi / 2;
end
if prop(i,1) == 0 && prop(i,2) < 0
prop_sph(i,2) = (-pi / 2) + 2*pi;
end
if prop(i,1) == 0 && prop(i,2) == 0
prop_sph(i,2) = NaN;
end
end
% plot the photon distribution in Cartesian coordinates.
% scatter3(prop(:,1),prop(:,2),prop(:,3), 10,'k','filled','MarkerFaceAlpha',.9);
% sort photons into bins of size d_th*d_phi in spherical coordinates
max_th = pi;
max_phi = 2*pi;
d_th = max_th / 100;
d_phi = max_phi / 100;
l_th = round((max_th / d_th) * (max_phi / d_phi));
theta = (0:d_th:max_th)';
phi = (0:d_phi:max_phi)';
mat = [zeros(l_th,1)'; zeros(l_th,1)'; zeros(l_th,1)';]';
n = 1;
for j = 1:length(phi)
for i = 1:length(theta)
mat(n,1) = theta(i);
mat(n,2) = phi(j);
n = n+1;
end
end
for w = 1:length(mat)
for k = 1:length(prop_sph)
if (((prop_sph(k,1) >= mat(w,1)) && (prop_sph(k,1) < mat(w,1) + d_th)) && ((prop_sph(k,2) >= mat(w,2)) && (prop_sph(k,2) < mat(w,2) + d_phi)))
mat(w,3) = mat(w,3) + 1;
end
end
end
% ------- convert back to cartesian coordinates and include color data
for i = 1:length(mat)
prop_map_fin(i,1) = prop_sph(i,3)*sin(mat(i,1) + d_th/2)*cos(mat(i,2) + d_phi/2);
prop_map_fin(i,2) = prop_sph(i,3)*sin(mat(i,1) + d_th/2)*sin(mat(i,2) + d_phi/2);
prop_map_fin(i,3) = prop_sph(i,3)*cos(mat(i,1) + d_th/2);
prop_map_fin(i,4) = mat(i,3);
% X = prop_sph(:,3).*sin(prop_sph(:,1)).*cos(prop_sph(:,2));
% Y = prop_sph(:,3).*sin(prop_sph(:,1)).*sin(prop_sph(:,2));
% Z = prop_sph(:,3).*cos(prop_sph(:,1));
end
% --- visualize --------------------------------------------------------
%surf(prop_map_fin(:,1),prop_map_fin(:,2),prop_map_fin(:,3),prop_map_fin(:,4))
scatter3(prop_map_fin(:,1),prop_map_fin(:,2),prop_map_fin(:,3), 10,prop_map_fin(:,4), 'filled');
cb = colorbar; % create and label the colorbar
cb.Label.String = '';
% xslice = [5 9.9]; % define the cross sections to view
% yslice = 3;
% zslice = ([-3 0]);
%slice(prop_map_fin(:,1), prop_map_fin(:,2), prop_map_fin(:,3), prop_map_fin(:,4), xslice, yslice, zslice)
%
T = array2table(prop_map_fin);
x = T{:, 1}; y = T{:, 2}; z = T{:, 3}; c2 = T{:, 4};
ux = unique(x); uy = unique(y); uz = unique(z);
data_sorted = sortrows(T, 1:4);
v = reshape(data_sorted{:, 4}, length(ux),length(uy),length(uz));
Error using reshape
Number of elements must not change. Use [] as one of the size inputs to automatically calculate the appropriate size for that dimension.
isosurface(uz, uy, ux, v, 0);
%[X2,Y2,Z2] = meshgrid(prop_map_fin(:,1), prop_map_fin(:,2), prop_map_fin(:,3));
function [x y] = GetCircle(r, h, k, a, b)
t = linspace(a, b, 40);
x = r*cos(t) + h;
y = r*sin(t) + k;
end

采纳的回答

Milan Bansal
Milan Bansal 2024-4-30
Hi iontrap,
I understand that you want to divide the sphere into a number of segments and calculate the number of "photons" falling in each of those segments. Then, you wish to color the segments using a color map according to the number of "photons" in them and get a smooth sphere as output.
Please follow the steps given below to achieve this:
1.) From your code, it seems you have calculated the x, y, and z coordinates of "photons" (stored in the "prop" matrix) from the "input.txt" file, and then you are converting them to spherical coordinates in the for loop. This step can be performed using the cart2sph function, as shown in the code snippet below.
photons = readmatrix('input.txt');
r = 70; % in unit mm
prop = zeros(length(photons),3);
for i = 1:length(photons)
c(i) = ((photons(i,1)*photons(i,1)) + (photons(i,2)*photons(i,2)) + (photons(i,3)*photons(i,3)) - (r*r));
b(i) = 2*((photons(i,1)*photons(i,7)) + (photons(i,2)*photons(i,8)) + (photons(i,3)*photons(i,9)));
a(i) = ((photons(i,7)*photons(i,7)) + (photons(i,8)*photons(i,8)) + (photons(i,9)*photons(i,9)));
t(i) = (-b(i) + sqrt(b(i).^2 - 4*a(i)*c(i)))/(2*a(i));
prop(i,1) = photons(i,1) + t(i)*photons(i,7);
prop(i,2) = photons(i,2) + t(i)*photons(i,8);
prop(i,3) = photons(i,3) + t(i)*photons(i,9);
end
x = prop(:,1); y = prop(:,2); z = prop(:,3);
% Convert Cartesian coordinates to spherical coordinates
[azimuth, elevation, ~] = cart2sph(x, y, z);
2.) Divide the sphere into a number of surface segments with equal intervals. This can be done by dividing the elevation and azimuth as shown in the code snippet below:
% Define the number of segments in azimuth and elevation
nAzimuthSegments = 72; % 5-degree intervals
nElevationSegments = 36; % 5-degree intervals
% Find the range of azimuth and elevation
azimuthRange = linspace(-pi, pi, nAzimuthSegments+1);
elevationRange = linspace(-pi/2, pi/2, nElevationSegments+1);
3.) Calculate the number of "photons" in each of those segments.
% Initialize the count matrix
counts = zeros(nElevationSegments, nAzimuthSegments);
% Count photons in each segment
for i = 1:nElevationSegments
for j = 1: nAzimuthSegments
counts(i, j) = sum(azimuth >= azimuthRange(j) & azimuth < azimuthRange(j+1) & ...
elevation >= elevationRange(i) & elevation < elevationRange(i+1));
end
end
4.) Create a meshgrid to plot the surface for each segments. Convert this meshgrid to cartesian coordinates for plotting. Use sph2cart for this conversion.
% Generate a meshgrid for plotting that matches the dimensions of 'counts'
[theta, phi] = meshgrid(linspace(-pi, pi, nAzimuthSegments), linspace(-pi/2, pi/2, nElevationSegments));
% Convert to Cartesian coordinates for plotting
[xPlot, yPlot, zPlot] = sph2cart(theta, phi, r);
5.) Get the RGB triplets to color the segments according to the number of "photons" stored in "counts" variable for each segment.
% Normalize 'counts' to the range [0, 1] for colormap mapping
countsNormalized = counts / max(counts(:));
N = 256;
colormapMatrix = parula(N);
ind = round(countsNormalized * (N-1)) + 1;
% Initialize the RGB matrix
rgbMatrix = zeros(size(counts, 1), size(counts, 2), 3);
% Fill the RGB matrix using the colormap
for i = 1:size(counts, 1)
for j = 1:size(counts, 2)
rgbMatrix(i, j, :) = colormapMatrix(ind(i, j), :);
end
end
6.) Plot surface plot for each segment with the respective color. This can be done using the surf function as shown in the code snippet below.
figure;
surf(xPlot, yPlot, zPlot, rgbMatrix, 'EdgeColor', 'none');
colormap('parula'); % Use a heat map color scheme
colorbar("TickLabels", linspace(0,250,11 ));
axis equal;
title('Photon Distribution on Sphere');
Please refer to the following documentation links to learn more about cart2sph, sph2cart and surf functions.
Hope this helps!
  2 个评论
iontrap
iontrap 2024-4-30
Thanks so much @Milan Bansal for cleaning up my approach.
This is simpler than my attempt - I did not know that meshgrid could be used to provide a spherical grid. Hence why I was trying to convert back to cartesian prior to the call of meshgrid.
Thanks
iontrap
iontrap 2024-5-1
There is an interesting tradeoff between resolution and each segment's position-dependent area. Near the poles, when using a grid of ~ 300 x 300, small segments result in lower counts and a speckle-like pattern.
Some smoothing would be nice, or maybe I can normalize for surface area variation. Thanks.

请先登录,再进行评论。

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Contour Plots 的更多信息

产品


版本

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by