根据面积均分球面 & c++示意代码

参考文献：Malkin Z. A new method to subdivide a spherical surface into equal-area cells[J]. arXiv preprint arXiv:1612.03467, 2016.

一种根据面积均分球面的方法是根据经纬度来分，即将直角坐标$(x,y,z)$转换成极坐标$(\theta ,\varphi ,r)$，然后将$\theta$和$\varphi$均分。均分$\varphi$（经度）确实能让分出来的小块面积相等，但$\theta$（维度）分块后是不均匀的，如下图

网上找到一篇博客介绍了两种均分方法：正多面体（柏拉图体regular polyhedron）、Hierarchical Equal Area isoLatitude Pixelization（这个方法甚至专门有个c++库）

但是这两种方法感觉在原理和实现上都很复杂，于是经过一番搜寻找到了论文16年的一篇论文《A new method to subdivide a spherical surface into equal-area cells》

方法

这篇论文的方法跟开头介绍的错误分法都是从经纬度（$\theta$和$\varphi$）下手，但它并不是均分。主要步骤分两步：

• 首先将纬度划分成$x$块
  • 相邻块之间的宽度不完全一样，但是大致相同，而且整个划分关于赤道对称
• 将每一个纬度带$i(i=1,2,...,x)$，均分成$y_i$块。
  • 不同纬度带划分的块数不一样。纬度越高块数越少，赤道附近的块接近正方形
  • 同样关于赤道对称
  • 要求划分完成后每一块的面积相等

划分的方法论文中并没有给，只给了将球面均分成46、130、406块的划分方式


由于我手边工作的需要，这些划分数量太多了，我就自己手推了$x=4$和$x=5$的划分方式，而且我并没有要求“赤道附近的块接近正方形”

$x=4$的划分方式

• 因为关于赤道对称，上半球有两个纬度带，最高纬度带对应角度为$\theta$，第二个为$\pi /2-\theta$
• 最高纬度带的面积可以用球冠面积公式得到
  $$S_1=2\pi rh=2\pi r(1-cos\theta )r=2\pi r^2(1-cos\theta )$$
• 根据半球面积可以计算得到第二个纬度带的面积
  $$S_2=S_{半球}-S_1=2\pi r^2-2\pi r^2(1-cos\theta )=2\pi r^{2}cos\theta$$
• 假设第一个纬度带分$m$块，第二个纬度带分$n$块，则要求
  $$\begin{gathered} \frac{2\pi r^2(1-cos\theta )}{m}=\frac{2\pi r^{2}cos\theta }{n} \\ \frac{(1-cos\theta )}{cos\theta }=\frac{m}{n} \end{gathered}$$
• （可选）如果想要让赤道附近的块尽量接近正方形，我们可以大致估算一下$n$的范围
  • 赤道附近的每个块的高度（纬度带的宽度）为$h=r\ast (\pi /2-\theta )$
  • 赤道附近的纬度带分$n$块，则我们想要上下边的长度尽量接近高度，即
    $$\begin{gathered} l_{上}=\frac{2\pi rsin\theta }{n}\approx r\ast (\pi /2-\theta )\to n=\frac{2\pi sin\theta }{\pi /2-\theta } \\ {l}_{下}=\frac{2\pi r}{n}\approx r\ast (\pi /2-\theta )\to n=\frac{2\pi }{\pi /2-\theta } \end{gathered}$$
  • 将这两条线画出来可以看到$n$越大，越接近正方形
• （可选）如果想要纬度带宽度大致相同，即要求$\theta \approx \pi /4$，带入之前的公式得
  $$\frac{m}{n}=\frac{(1-cos\theta )}{cos\theta }\approx \frac{(1-cos\pi /4)}{cos\pi /4}\approx 0.414$$
• 接下来只要按照想要的划分方式给$m$和$n$赋整数，就能求解$\theta$，例如
  • $m=2,n=3,\theta =53.130°$
  • $m=3,n=4,\theta =55.150°$
  • $m=3,n=5,\theta =51.318°$
  • $m=4,n=5,\theta =56.217°$
  • $m=3,n=6,\theta =48.190°$
  • $m=3,n=7,\theta =45.573°$

Total Cell	Cells	Latitudinal	Longitudinal
10	2 3 3 2	0-53.13° 53.13°-90° 90°-126.87° 126.87°-180°	180° 120° 120° 180°
14	3 4 4 3	0-55.15° 55.15°-90° 90°-124.85° 124.85°-180°	120° 90° 90° 120°
18	4 5 5 4	0-56.22° 56.22°-90° 90°-123.78° 123.78°-180°	90° 72° 72° 90°
20	3 7 7 3	0-45.573° / 0 - 0.7954 45.573°-90° / 0.7954 - 1.5708 90°-134.427° / 1.5708 - 2.3462 134.427°-180° / 2.3462 - 3.1416	120° / 2.0944 51.43° / 0.8976 51.43° / 0.8976 120° / 2.0944

$x=5$的划分方式

• 虽然$x=5$为奇数，但是因为关于赤道对称，这里算上半球有2.5个纬度带，即第三个纬度带只算一半。三个纬度带从高到底对应的角度分别为$\theta ,\varphi ,\pi /2-\theta -\varphi$
• 最高纬度带的面积可以用球冠面积公式得到
  $$S_1=2\pi rh=2\pi r(1-cos\theta )r=2\pi r^2(1-cos\theta )$$
• 根据球冠面积同样可以计算得到第二个纬度带的面积
  $$S_2=2\pi r^2(1-cos(\theta +\varphi ))-2\pi r^2(1-cos\theta )=2\pi r^2(cos\theta -cos(\theta +\varphi ))$$
• 根据半球面积可以计算得到最后半个纬度带的面积
  $$S_{2.5}=S_{半球}-(S_1+S_2)=2\pi r^2-2\pi r^2(1-cos(\theta +\varphi ))=2\pi r^{2}cos(\theta +\varphi )$$
• 假设第一个纬度带分$m$块，第二个纬度带分$n$块，最后半个纬度带分$o$块，则要求
  $$\begin{gathered} \frac{2\pi r^2(1-cos\theta )}{m}=\frac{2\pi r^2(cos\theta -cos(\theta +\varphi ))}{n}=2\ast \frac{2\pi r^{2}cos(\theta +\varphi )}{o} \\ \frac{1-cos\theta }{m}=\frac{cos\theta -cos(\theta +\varphi )}{n}=\frac{2cos(\theta +\varphi )}{o} \end{gathered}$$
• 接下来只要按照想要的划分方式给$m$、$n$和$o$赋整数，就能求解$\theta$和$\varphi$，例如
  • $m=3,n=4,o=5,\theta =46.826°,\varphi =27.916°$
  • $m=3,n=4,o=6,\theta =45.573°,\varphi =26.969°$
  • $m=4,n=5,o=6,\theta =48.190°,\varphi =27.333°$

Total Cell	Cells	Latitudinal	Longitudinal
15	1 3 7 3 1	0.0000 - 0.5223 0.5223 - 1.0853 1.0853 - 2.0563 2.0563 - 2.6193 2.6193 - 3.1416	360° / 6.2832 120° / 2.0944 51.43° / 0.8976 120° / 2.0944 360° / 6.2832
19	3 4 5 4 3	0-46.826° 46.826°-74.743° 74.743°-105.257° 105.257°- 133.174° 133.174°-180°	120° 90° 72° 90° 120°
19	1 5 7 5 1	0.0000 - 0.4630 0.4630 - 1.1931 1.1931 - 1.9485 1.9485 - 2.6786 2.6786 - 3.1416	360° / 6.2832 72° / 1.2567 51.43° / 0.8976 72° / 1.2567 360° / 6.2832
24	4 5 6 5 4	0-48.190° 48.190°-75.522° 75.522°-104.478° 104.478°- 131.810° 131.810°-180°	90° 72° 60° 72° 90°
25	1 7 9 7 1	0.0000 - 0.4027 0.4027 - 1.2025 1.2025 - 1.9391 1.9391 - 2.7389 2.7389 - 3.1416	360° / 6.2832 51.43° / 0.8976 40° / 0.6981 51.43° / 0.8976 360° / 6.2832
27	1 7 11 7 1	0.0000 - 0.3873 0.3873 - 1.1512 1.1512 - 1.9904 1.9904 - 2.7543 2.7543 - 3.1416	360° / 6.2832 51.43° / 0.8976 32.72° / 0.5712 51.43° / 0.8976 360° / 6.2832

c++代码

这里以“统计球形点云落在每个块中的数量”为例，划分方式为$x=4,m=3,n=7$

/**
 * 统计球形点云落在每个块中的数量
 * @param[in] points n*3的矩阵，存储三维点云
 * @param[out] nums  统计结果，存储顺序北极->南极 & 经度0°->360°
 */
void ClassifyCentroidVector(const Eigen::MatrixXf &points, vector<int> &nums) {
  direct.resize(20);
  Eigen::Vector4f block_theta(0., 0.7954, 1.5708, 2.3462);
  Eigen::Vector4f block_phi(2.0944, 0.8976, 0.8976, 2.0944);
  Eigen::Vector4i offset(0, 3, 10, 17);
  for (int k = 0; k < points.rows(); k++) {
    //笛卡尔坐标系转换为球坐标系
    const float x = points(k, 0), y = points(k, 1), z = points(k, 2);
    float r = sqrt(x * x + y * y + z * z);
    float theta = acos(z / r);
    float fai = atan2(y, x) + M_PI; // 取值范围0-2π
    
    int coord_theta = -1;
    for(int i = 0; i < block_theta.size(); i++)
      if(theta > block_theta(i)){
        coord_theta = i;
        break;
      }
    assert(coord_theta != -1);
    int coord_phi = int(fai/block_phi(coord_theta));
    direct[offset(coord_theta)+coord_phi] += 1;
  }
}