3

I am trying to create a spherical histogram and I have found an addon that creates a spherical histogram, however it uses equal area quadrilaterals and for my purposes I would like to preserve the lines found on the traditional sphere created with MATLAB so they can match latitude and longitude lines.

I found a way to color a given subset of the sphere by setting the min/max Azimuth and Elevation values of the region I want to shade, and to try and test out coloring each cell of the sphere I have tried dividing the azimuth and elevation total angles by a given n value and looping through to assign a random color to each cell. This works for about 2/3s of the cells, however there are random azimuth/elevation sections where there is no color assigned, indicating some type of problem with the surf() function presumably. I thought the count n for the sphere being only 20 might cause rounding errors that would be the main contributing factor for this, but there are proportionally similar sized gaps found for a 50x50 cell sphere as well. My best guess is some kind of rounding precision error that causes the given cell region to skip assigning a color due to the given bounds being too far from matching an actual cell, essentially having the given bounds being in between two lines. How can I iterate through every cell based on its azimuth and elevation range for a given sphere of size n?



n = 20;

%// Change your ranges here
minAzimuth = -180;
maxAzimuth = 180;
minElevation = 0;
maxElevation = 180;

%// Compute angles - assuming that you have already run the code for sphere
[x,y,z] = sphere(n);
theta = acosd(z);
phi = atan2d(y, x);

%%%%%// Begin highlighting logic
ind = (phi >= minAzimuth & phi <= maxAzimuth) & ...
      (theta >= minElevation & theta <= maxElevation); % // Find those indices
x2 = x; y2 = y; z2 = z; %// Make a copy of the sphere co-ordinates
x2(~ind) = NaN; y2(~ind) = NaN; z2(~ind) = NaN; %// Set those out of bounds to NaN

%%%%%// Draw our original sphere and then the region we want on top
r = 1;
surf(r.*x,r.*y,r.*z,'FaceColor', 'white', 'FaceAlpha',0); %// Base sphere
hold on;
%surf(r.*x2,r.*y2,r.*z2,'FaceColor','red', 'FaceAlpha',0.5); %// Highlighted portion
 %// Adjust viewing angle for better view



for countAz = 1:1:n
    current_minAzimuth = -180 + (countAz-1) * (360/n);
    current_maxAzimuth = -180 + (countAz) * (360/n);
    for countE = 1:1:n
        current_minElevation = 0 + (countE-1) * (180/n);
        current_maxElevation = 0 + (countE) * (180/n);
        
        theta = acosd(z);
        phi = atan2d(y, x);

        %%%%%// Begin highlighting logic
        ind = (phi >= current_minAzimuth & phi <= current_maxAzimuth) & ...
            (theta >= current_minElevation & theta <= current_maxElevation); % // Find those indices
        x2 = x; y2 = y; z2 = z; %// Make a copy of the sphere co-ordinates
        x2(~ind) = NaN; y2(~ind) = NaN; z2(~ind) = NaN;
        
        random_color = rand(1,3);

        surf(r.*x2,r.*y2,r.*z2,'FaceColor',random_color, 'FaceAlpha',1);

    end

end

axis equal;
view(40,40);


Here is an n=50 sphere: n=50 sphere

And here is an n=20 sphere

n=20 sphere

Wolfie
  • 27,562
  • 7
  • 28
  • 55
yes
  • 67
  • 4

1 Answers1

2

You're doing more conversions to and from spherical/Cartesian coordinates than you need to, and subsequently doing more floating point comparisons than required.

You are essentially just looping through all of the 2x2 index blocks in the x, y and z arrays.

It's simpler to just directly create the index arrays you're trying to use and then pull out those values of the existing surface coordinates.

for countAz = 1:1:n
    for countE = 1:1:n
        % Linear index for a 2x2 patch of the matrices
        idx = countAz + [0,1,n+(1:2)] + (countE-1)*(n+1); 
        % Pull out the coordinates, reshape for surf
        x2 = x(idx); x2 = reshape(x2,2,2);
        y2 = y(idx); y2 = reshape(y2,2,2);
        z2 = z(idx); z2 = reshape(z2,2,2);
                        
        random_color = rand(1,3);
        surf(r*x2,r*y2,r*z2,'FaceColor',random_color, 'FaceAlpha',1);
    end
end

surface

Wolfie
  • 27,562
  • 7
  • 28
  • 55
  • This does fix most of the blank surfaces, however there is one sliver of cells along a particular Azimuth(on the back side of the figure in your picture) that is still blank, and I don't think it's due to the index not reaching there. – yes Nov 22 '22 at 22:19
  • @yes Maybe because my suggested indexing misses the case where you have one pair of coords on the left edge of the matrices and one pair on the right edge, "wrapping around"? You could handle this as a special case after the loop... – Wolfie Nov 22 '22 at 22:42
  • hey @Wolfie! I believe your `idx` is incorrect here. Also, there is no need to `wrap around`, `x,y,z` are already repeated in boundary. I have not tried to fix `idx`, but I believe this is what you were trying to make `idx` do, yet its just almost right: https://stackoverflow.com/questions/74540426/matlab-sphere-region-coloring-based-on-azimuth-and-elevation-input-not-index#74540426 . I suspect the bug is that there are `n` patches, but `n+1` points. – Ander Biguri Nov 23 '22 at 11:21
  • 1
    Thanks @AnderBiguri, indeed fixing the column offset to be `n+1` instead of `n` sorts out the indexing, and means the weird offsets shown in the linked question don't happen, i.e. you can fix the az or el and it plots correctly. Edited – Wolfie Nov 23 '22 at 11:32