球面三角网格绘制算法（附OpenGL代码）

本文给出的代码都是在原点处、半径为1的球面绘制代码。

对于以任意点为球心、任意长度为半径的情况，可以通过使用glTranslated实现球的平移、glScaled实现球的缩放。

编筐法

如图所示，球的绘制是一层一层完成的，如同编一个竹筐。

将球分为N层，所以共N+1条纬线，除了最顶上和最底下的两个纬线以外（因为已经退化为一个点），每条纬线均分M等份，相邻纬线交错均分，该处的均分指的是角度均分。

如果只是要绘制如图所示的线框球，将每一层上的两个点与下一层的对应的点用三角形连线起来即可；如果要绘制完整球面，则还需要反向将下一层的每两个点与本层连接。

尽管M的取值与N无关，但还是建议取M=2N，这样球体在纬度和经度方向上均匀程度相同。

下面是N=3和N=20的效果：

N=5	N=20

该方法的优点是简单，缺点是图元的分布不均匀，两极的图元密度比赤道的密度大很多。

下面是完整代码：

/**
 * @brief 扎篮子法
 * @param resolution 分辨率
 */
void
DrawSphere_1(GLuint resolution)
{
  glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);

  GLdouble dTheta = kPi / resolution;

  struct XY
  {
    GLdouble x, y;
  };
  auto a = new XY[2 * resolution + 1];
  auto b = new XY[2 * resolution + 1];

  GLdouble thetaZ = dTheta;  // 纬转角
  GLdouble thetaXY = 0;      // 经转角
  GLdouble az = cos(thetaZ); // a纬面高
  GLdouble bz;               // b纬面高
  GLdouble rr = sin(thetaZ); // 纬半径

  // 绘制顶盖
  glBegin(GL_TRIANGLE_FAN);
  glVertex3d(0, 0, 1);
  for (auto i = 0U; i < 2 * resolution; ++i) {
    a[i].x = cos(thetaXY) * rr;
    a[i].y = sin(thetaXY) * rr;
    thetaXY += dTheta;
    glVertex3d(a[i].x, a[i].y, az);
  }
  a[2 * resolution] = a[0];
  glVertex3d(a[2 * resolution].x, a[2 * resolution].y, az); // 封闭纬面
  glEnd();

  // 绘制中间层
  for (auto i = 1U; i < resolution - 1; ++i) {
    thetaZ += dTheta;
    bz = cos(thetaZ) * 1;
    rr = sin(thetaZ) * 1;
    thetaXY = dTheta / 2 * i;

    glBegin(GL_TRIANGLES);
    for (auto j = 0U; j < 2 * resolution; ++j) {
      b[j].x = cos(thetaXY) * rr;
      b[j].y = sin(thetaXY) * rr;
      thetaXY += dTheta;

      glVertex3d(a[j].x, a[j].y, az);
      glVertex3d(a[j + 1].x, a[j + 1].y, az);
      glVertex3d(b[j].x, b[j].y, bz);
    }
    b[2 * resolution] = b[0];
    glEnd();

    auto tmp = a;
    a = b;
    b = tmp;
    az = bz;
  }

  // 绘制底盖
  glBegin(GL_TRIANGLE_FAN);
  glVertex3d(0, 0, -1);
  for (auto i = 0U; i < 2 * resolution; ++i) {
    glVertex3d(a[i].x, a[i].y, az);
  }
  glVertex3d(a[2 * resolution].x, a[2 * resolution].y, az); // 封闭纬面
  glEnd();

  delete[] a;
  delete[] b;
}

面细分法

如上图所示，通过将原来的三角面划分为更细小的面即可实现对球面的逼近。

该算法的核心是细分操作，一次细分操作将一个三角面划分为四个更小的三角面。给定一个顶点都在球面上的三个点a、b、c，依次求出角aob、boc、coa的平分线与球面的交点d、e、f，则新的四个三角形为：afd、bde、cef、def。

细分的起点可以是四面体、八面体、十二面体等等。从理论上说，任何一个封闭的、由三角形构成的凸多面体都可以收敛为一个球，但是建议采用八面体作为细分的起点，因为八面体的顶点坐标非常简单，且由于其对称性，只需计算第一象限的点，然后做镜像变换即可得到整个球。示例图片就是第一象限的细分过程。

细分的终点可以是多种多样的，例如，可以当三角面的面积小于某个给定值时就不再细分，也可以当三角形的边长不长于某个给定值时就不再细分，还可以直接指定细分的次数。

该方法的一个好处就是可以得到较均匀的球面近似结果，例如，如果约束条件为三角形的边长不大于X，那么最终结果绘制的线端之间的长度之差不会超过X/2，使用面积也同理。

一个较高精细度下的效果，这种方法绘制出的球体极性不明显，看着有点晕：

下面的函数midArcPoint计算角平分线与球面的交点（也即球面上弧的中点）：

struct XYZ
{
  GLdouble x, y, z;
};

XYZ
midArcPoint(const XYZ& a, const XYZ& b)
{
  XYZ c{ a.x + b.x, a.y + b.y, a.z + b.z };
  GLdouble mod = sqrt(c.x * c.x + c.y * c.y + c.z * c.z);
  c.x /= mod, c.y /= mod, c.z /= mod;
  return c;
}

下面是一种绘制代码：

#include <queue>

GLdouble
distSquare(const XYZ& a, const XYZ& b)
{
  GLdouble dx = a.x - b.x, dy = a.y - b.y, dz = a.z - b.z;
  return dx * dx + dy * dy + dz * dz;
}

/**
 * @brief 面细分法
 * @param resolution 分辨率
 */
void
DrawSphere_2(GLdouble resolution)
{
  glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);

  struct Triangle
  {
    XYZ a, b, c;
  };

  std::queue<Triangle> triangles;
  triangles.push({ 1, 0, 0, 0, 1, 0, 0, 0, 1 });

  resolution *= resolution;
  while (!triangles.empty()) {
    auto& t = triangles.front();

    // 当三角形各边长度都不大于resolution时就不再进一步细分
    if (distSquare(t.a, t.b) > resolution ||
        distSquare(t.b, t.c) > resolution ||
        distSquare(t.c, t.a) > resolution) {
      auto d = midArcPoint(t.a, t.b), e = midArcPoint(t.b, t.c),
           f = midArcPoint(t.c, t.a);
      triangles.push({ t.a, f, d });
      triangles.push({ t.b, d, e });
      triangles.push({ t.c, e, f });
      triangles.push({ d, e, f });
    } else {
      glBegin(GL_TRIANGLES);

      // 第一象限
      glVertex3dv((GLdouble*)&t.a);
      glVertex3dv((GLdouble*)&t.b);
      glVertex3dv((GLdouble*)&t.c);

      // 第二象限
      t.a.x = -t.a.x, t.b.x = -t.b.x, t.c.x = -t.c.x;
      glVertex3dv((GLdouble*)&t.a);
      glVertex3dv((GLdouble*)&t.b);
      glVertex3dv((GLdouble*)&t.c);

      // 第三象限
      t.a.y = -t.a.y, t.b.y = -t.b.y, t.c.y = -t.c.y;
      glVertex3dv((GLdouble*)&t.a);
      glVertex3dv((GLdouble*)&t.b);
      glVertex3dv((GLdouble*)&t.c);

      // 第四象限
      t.a.x = -t.a.x, t.b.x = -t.b.x, t.c.x = -t.c.x;
      glVertex3dv((GLdouble*)&t.a);
      glVertex3dv((GLdouble*)&t.b);
      glVertex3dv((GLdouble*)&t.c);

      // 第五象限
      t.a.z = -t.a.z, t.b.z = -t.b.z, t.c.z = -t.c.z;
      glVertex3dv((GLdouble*)&t.a);
      glVertex3dv((GLdouble*)&t.b);
      glVertex3dv((GLdouble*)&t.c);

      // 第六象限
      t.a.x = -t.a.x, t.b.x = -t.b.x, t.c.x = -t.c.x;
      glVertex3dv((GLdouble*)&t.a);
      glVertex3dv((GLdouble*)&t.b);
      glVertex3dv((GLdouble*)&t.c);

      // 第七象限
      t.a.y = -t.a.y, t.b.y = -t.b.y, t.c.y = -t.c.y;
      glVertex3dv((GLdouble*)&t.a);
      glVertex3dv((GLdouble*)&t.b);
      glVertex3dv((GLdouble*)&t.c);

      // 第八象限
      t.a.x = -t.a.x, t.b.x = -t.b.x, t.c.x = -t.c.x;
      glVertex3dv((GLdouble*)&t.a);
      glVertex3dv((GLdouble*)&t.b);
      glVertex3dv((GLdouble*)&t.c);
      glEnd();
    }
    triangles.pop();
  }
}