MATLAB球面经纬网格着色出现未填充区域的解决问询

想要创建匹配MATLAB传统球面经纬线的球面直方图，现有插件生成的是等面积四边形，不符合需求。尝试通过划分方位角（Azimuth）和仰角（Elevation）范围，循环为每个网格单元分配随机颜色，但无论n取20还是50，均有约1/3区域未被着色，推测是精度问题导致surf()函数未匹配到网格单元。请问如何遍历给定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);

效果说明：

• n=50时：球面存在大量未着色空白区域，仅约2/3的经纬网格被随机颜色填充
• n=20时：同样存在明显未着色区域，空白占比约1/3

问题根源 #

1. 边界精度误差：atan2d(y,x)生成的方位角phi范围是[-180,180]，原代码用左闭右闭区间判断，当区间上限等于180时，浮点数精度问题会导致phi=180的点无法被匹配。
2. 网格数量不匹配：sphere(n)生成的是(n+1)×(n+1)的网格点，对应n×n个面，原逻辑未完全覆盖所有网格单元。
3. 冗余计算：循环内重复计算theta和phi，额外引入精度误差并降低效率。

修正后的代码 #

n = 20;
r = 1;

% 提前生成球面网格与对应角度，避免循环内重复计算
[x,y,z] = sphere(n);
theta = acosd(z);    % 仰角范围[0,180]
phi = atan2d(y, x);  % 方位角范围[-180,180]

% 绘制基础球面
surf(r.*x,r.*y,r.*z,'FaceColor', 'white', 'FaceAlpha',0);
hold on;
axis equal;
view(40,40);

% 遍历所有n×n个球面网格单元
for countAz = 1:n
    % 方位角区间采用左闭右开，避免边界精度问题
    current_minAz = -180 + (countAz-1)*(360/n);
    current_maxAz = -180 + countAz*(360/n);
    
    % 处理最后一个方位角区间，确保包含phi=180的点
    if countAz == n
        current_maxAz = 180 + eps;
    end
    
    for countE = 1:n
        % 仰角区间同样采用左闭右开
        current_minEl = 0 + (countE-1)*(180/n);
        current_maxEl = 0 + countE*(180/n);
        
        % 匹配当前区间内的网格点
        ind = (phi >= current_minAz & phi < current_maxAz) & ...
              (theta >= current_minEl & theta < current_maxEl);
              
        % 复制网格数据并屏蔽非目标区域
        x2 = x; y2 = y; z2 = z;
        x2(~ind) = NaN; y2(~ind) = NaN; z2(~ind) = NaN;
        
        % 绘制当前网格单元
        surf(r.*x2,r.*y2,r.*z2,'FaceColor',rand(1,3), 'FaceAlpha',1);
    end
end

关键修正点 #

• 区间改为左闭右开：将判断逻辑从<= current_max改为< current_max，最后一个方位角区间额外加eps，确保180°边界点被正确匹配，解决浮点数精度问题。
• 匹配网格数量：sphere(n)生成(n+1)个点对应n个网格单元，循环次数保持n次即可覆盖所有区域。
• 提前计算角度：将theta和phi的计算移到循环外，减少误差并提升运行效率。

内容的提问来源于stack exchange，提问作者yes