[C++知识库]根据面积均分球面 & 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 ,?,r)$，然后将$\theta$和$?$均分。均分$?$（经度）确实能让分出来的小块面积相等，但$\theta$（维度）分块后是不均匀的，如下图

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

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

方法

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

• 首先将纬度划分成$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?(\pi /2?\theta )$
  • 赤道附近的纬度带分$n$块，则我们想要上下边的长度尽量接近高度，即
    $$\begin{gathered} l_{上}=\frac{2\pi rsin\theta }{n}\approx r?(\pi /2?\theta )\to n=\frac{2\pi sin\theta }{\pi /2?\theta } \\ {l}_{下}=\frac{2\pi r}{n}\approx r?(\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°$

$x=5$的划分方式

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

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;
  }
}