4D spherical shell heat map

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

Answer 1

[Milan Bansal]

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.

• https://www.mathworks.com/help/matlab/ref/cart2sph.html
• https://www.mathworks.com/help/matlab/ref/sph2cart.html
• https://www.mathworks.com/help/matlab/ref/surf.html

Hope this helps!

[iontrap]

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]

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.