LS3/NS3球体生成算法及其实现

给定球体的中心及其半径。您的任务是有效地存储在像素的计算机屏幕上显示球体所需的所有整数点，并且该算法将通过其幼稚弧创建体素来生成幼稚球体。 NS3代表平方和的朴素球体，LS3代表平方和的格子球体。

例子：

Input : Center(0, 0, 0) Radius 1
Output : (-1, 0, 0), (0, 0, 1), (0, 0, -1), (0, 1, 0), 
         (0, -1, 0), (1, 0, 0)
         6(Total no. of points)

Input : Center(0, 0, 0) Radius 2
Output :(-2, 0, 0)(-1, 1, 1)(-1, -1, 1)(-1, 1, -1)
        (-1, -1, -1)(0, 0, 2)(0, 0, -2)(0, 2, 0)
        (0, -2, 0)(0, 1, 2)(0, 2, 1)(0, -1, 2)
        (0, 1, -2)(0, -1, -2)(0, -2, 1)(0, 2, -1)
        (0, -2, -1)(1, 2, 0)(1, 0, 2)(1, -2, 0)
        (1, 2, 0)(1, -2, 0)(1, 0, -2)(1, 0, 2)
        (1, 0, -2)(1, 1, 1)(1, -1, 1)(1, 1, -1)
        (1, -1, -1)(2, 0, 0)(2, 0, 1)(2, 0, 1)
        (2, 0, -1)(2, 0, -1)(2, 1, 0)(2, -1, 0)
        (2, 1, 0)(2, -1, 0) 
         38 (Total no. of points)

我强烈建议您，如果您在该领域没有任何先决条件，请转到以下链接。这将帮助您获得理解以下实现和算法所需的一些先决条件和概念。

算法 LS3 主要侧重于生成prima quadraginta的体素。 prima quadraginta 是将一个离散的球体分成 48 个部分，使得所有 48 个部分都具有某种对称关系，并且它们都可以通过使用其中的任何一个部分来生成时形成的。这一部分称为 prima quadraginta（见图 2）。类似地，quadraginta 八分圆的定义与 prima qudraginta 相同，只是它被划分为 8 个相等的部分，而不是 48 个。 quadaginta 八分圆包含 6 个 prima quadraginta，如下图所示。 Q 表示四边形（q-octant），C 表示四边形（c-octant）。一个朴素的弧是那些在四边形中具有相同 x 坐标的体素集（在上图中，红色体素正在形成一个天真的弧）。

借助48对称性质，可以从它的初四边形发展出一个完整的离散球体。使用 48 对称的限制是 q-octant 的一些体素将重复，因为它们与其他相邻的 prima quadraginta 共享一个边缘或一个体素。例如，在上图中，以灰色显示的体素属于两个或多个 q-octants 。

我们的算法会产生一些重复的体素，因此我们需要对输入体素做一些限制。
为了做到这一点，我在我的实现中包含了一些功能，它们是：

所有这些方法主要关注灰色体素（见上图），除了一个 putpixelall()，它将通过使用 48 对称性而没有任何限制地产生除灰色体素之外的所有体素。

putpixeledge1 ：这个函数将产生所有具有一个零坐标的体素
如果我们继续使用对称性找到另一个体素并拍摄关于该零坐标的图像
轴（如果我们相对于 x 轴拍摄图像，则此 (0, 2, -1) 体素将具有相同的图像）然后它来了
成为相同的体素。

putpixeledge2 ：这些具有 x 和 y 坐标的体素集是相同的（参见 prima quadraginta 图 2）。为了存储这些体素，我们必须寻找 q1 和 q2 八分圆之间的关系（参见上表）。您可以通过先交换 (x, y) 然后交换 (y, z) 来从 q1 获得 q2 中任何体素的一般坐标。

putpixelsingle ：这些像素的数量非常少。它们位于四边形八分圆（c 八分圆）的中间，并与位于同一 c 八分圆中的所有 q 八分圆共享。为了获得这些体素，我们可以看到上表并写出与其c-八分圆对应的所有体素。

putpixeldouble ：这些具有 y 和 z 坐标的体素集是相同的（参见 prima qudraginta 图 2）。为了存储这些体素，我们需要找到 q1 和 q6 八分圆之间的关系（见表）。您可以通过先交换 (y, z) 然后交换 (x, y) 来从 q1 获得 q6 中任何体素的一般坐标。

putpixelmiddle ：这些体素集位于正好有两个坐标的点上
value 为 NULL 可以帮助我们找到获取这些体素的坐标（你可以写坐标
这些体素手动）。

我们需要有效存储将由该算法形成的所有整数点（体素）。为了存储这些点，我们将使用 hashmap，hashmap 的键值将是点的 x 坐标，数据值将是 y 和 z 坐标对的链表。

#include 
#include 
  
// this header file contains get<0> and get<1> function
#include 
  
#include 
#include 
using namespace std;
  
// map to store the voxels
map > > mymap;
map > >::iterator itr;
  
// this function will store single voxel
void putpixelone(int m, int n, int p)
{
    mymap[m].push_back(make_pair(n, p));
}
  
// function to store some particular type of voxel
void putpixelmiddle(int a, int b, int c, int x, int y, int z)
{
    putpixelone(a + x, b + y, c + z);
    putpixelone(a + x, b + y, -c + z);
    putpixelone(a + x, c + y, b + z);
    putpixelone(a + x, -c + y, b + z);
    putpixelone(c + x, a + y, b + z);
    putpixelone(-c + x, a + y, b + z);
}
  
// This program will generate 24 voxels
void putpixeldouble(int a, int b, int c, int x, int y, int z)
{
    for (int j = 0; j < 3; j++) {
        putpixelone(a + x, b + y, c + z);
        swap(b, c);
        swap(a, b);
    }
    for (int j = 0; j < 3; j++) {
        putpixelone(-a + x, b + y, c + z);
        swap(b, c);
        swap(a, b);
    }
    for (int j = 0; j < 3; j++) {
        putpixelone(a + x, -b + y, c + z);
        swap(b, c);
        swap(a, b);
    }
    for (int j = 0; j < 3; j++) {
        putpixelone(-a + x, -b + y, c + z);
        swap(b, c);
        swap(a, b);
    }
    for (int j = 0; j < 3; j++) {
        putpixelone(a + x, b + y, -c + z);
        swap(b, c);
        swap(a, b);
    }
    for (int j = 0; j < 3; j++) {
        putpixelone(-a + x, b + y, -c + z);
        swap(b, c);
        swap(a, b);
    }
    for (int j = 0; j < 3; j++) {
        putpixelone(a + x, -b + y, -c + z);
        swap(b, c);
        swap(a, b);
    }
    for (int j = 0; j < 3; j++) {
        putpixelone(-a + x, -b + y, -c + z);
        swap(b, c);
        swap(a, b);
    }
}
  
// Explained above
void putpixelsingle(int a, int b, int c, int x,
                    int y, int z)
{
    putpixelone(a + x, b + y, c + z);
    putpixelone(-a + x, b + y, c + z);
    putpixelone(a + x, -b + y, c + z);
    putpixelone(-a + x, -b + y, c + z);
    putpixelone(a + x, b + y, -c + z);
    putpixelone(-a + x, b + y, -c + z);
    putpixelone(a + x, -b + y, -c + z);
    putpixelone(-a + x, -b + y, -c + z);
}
// This program will generate 24 voxels.
void putpixeledge2(int a, int b, int c, int x,
                   int y, int z)
{
  
    for (int j = 0; j < 3; j++) {
        putpixelone(a + x, b + y, c + z);
        swap(a, b);
        swap(b, c);
    }
    for (int j = 0; j < 3; j++) {
        putpixelone(-a + x, b + y, c + z);
        swap(a, b);
        swap(b, c);
    }
    for (int j = 0; j < 3; j++) {
        putpixelone(a + x, -b + y, c + z);
        swap(a, b);
        swap(b, c);
    }
    for (int j = 0; j < 3; j++) {
        putpixelone(-a + x, -b + y, c + z);
        swap(a, b);
        swap(b, c);
    }
    for (int j = 0; j < 3; j++) {
        putpixelone(a + x, b + y, -c + z);
        swap(a, b);
        swap(b, c);
    }
    for (int j = 0; j < 3; j++) {
        putpixelone(-a + x, b + y, -c + z);
        swap(a, b);
        swap(b, c);
    }
    for (int j = 0; j < 3; j++) {
        putpixelone(a + x, -b + y, -c + z);
        swap(a, b);
        swap(b, c);
    }
    for (int j = 0; j < 3; j++) {
        putpixelone(-a + x, -b + y, -c + z);
        swap(a, b);
        swap(b, c);
    }
}
// after excluding the repeated pixels from the set all 48
// pixels we will left with only 24 pixels so, we have
// defined all the voxels given below.
void putpixeledge1(int a, int b, int c, int x, int y, int z)
{
    putpixelone(a + x, b + y, c + z);
    putpixelone(a + x, c + y, b + z);
    putpixelone(a + x, -b + y, c + z);
    putpixelone(a + x, b + y, -c + z);
    putpixelone(a + x, -b + y, -c + z);
    putpixelone(a + x, -c + y, b + z);
    putpixelone(a + x, c + y, -b + z);
    putpixelone(a + x, -c + y, -b + z);
    putpixelone(b + x, c + y, a + z);
    putpixelone(b + x, a + y, c + z);
    putpixelone(b + x, -c + y, a + z);
    putpixelone(b + x, c + y, -a + z);
    putpixelone(b + x, -c + y, -a + z);
    putpixelone(b + x, a + y, -c + z);
    putpixelone(b + x, -a + y, c + z);
    putpixelone(b + x, -a + y, -c + z);
    putpixelone(c + x, a + y, b + z);
    putpixelone(c + x, -a + y, b + z);
    putpixelone(c + x, a + y, -b + z);
    putpixelone(c + x, -a + y, -b + z);
    putpixelone(c + x, b + y, a + z);
    putpixelone(c + x, -b + y, a + z);
    putpixelone(c + x, b + y, -a + z);
    putpixelone(c + x, -b + y, -a + z);
}
  
// voxel formation by 48 symmetry.
void putpixelall(int a, int b, int c, int x,
                 int y, int z)
{
    putpixelone(a + x, b + y, c + z);
    putpixelone(b + x, a + y, c + z);
    putpixelone(c + x, b + y, a + z);
    putpixelone(a + x, c + y, b + z);
    putpixelone(c + x, a + y, b + z);
    putpixelone(b + x, c + y, a + z);
    putpixelone(-a + x, -b + y, c + z);
    putpixelone(-b + x, -a + y, c + z);
    putpixelone(c + x, -b + y, -a + z);
    putpixelone(-a + x, c + y, -b + z);
    putpixelone(c + x, -a + y, -b + z);
    putpixelone(-b + x, c + y, -a + z);
    putpixelone(a + x, -b + y, -c + z);
    putpixelone(-b + x, a + y, -c + z);
    putpixelone(-c + x, -b + y, a + z);
    putpixelone(a + x, -c + y, -b + z);
    putpixelone(-c + x, a + y, -b + z);
    putpixelone(-b + x, -c + y, a + z);
    putpixelone(-a + x, b + y, -c + z);
    putpixelone(b + x, -a + y, -c + z);
    putpixelone(-c + x, b + y, -a + z);
    putpixelone(-a + x, -c + y, b + z);
    putpixelone(-c + x, -a + y, b + z);
    putpixelone(b + x, -c + y, -a + z);
    putpixelone(-a + x, b + y, c + z);
    putpixelone(b + x, -a + y, c + z);
    putpixelone(c + x, b + y, -a + z);
    putpixelone(-a + x, c + y, b + z);
    putpixelone(c + x, -a + y, b + z);
    putpixelone(b + x, c + y, -a + z);
    putpixelone(a + x, -b + y, c + z);
    putpixelone(-b + x, a + y, c + z);
    putpixelone(c + x, -b + y, a + z);
    putpixelone(a + x, c + y, -b + z);
    putpixelone(c + x, a + y, -b + z);
    putpixelone(-b + x, c + y, a + z);
    putpixelone(a + x, b + y, -c + z);
    putpixelone(b + x, a + y, -c + z);
    putpixelone(-c + x, b + y, a + z);
    putpixelone(a + x, -c + y, b + z);
    putpixelone(-c + x, a + y, b + z);
    putpixelone(b + x, -c + y, a + z);
    putpixelone(-a + x, -b + y, -c + z);
    putpixelone(-b + x, -a + y, -c + z);
    putpixelone(-c + x, -b + y, -a + z);
    putpixelone(-a + x, -c + y, -b + z);
    putpixelone(-c + x, -a + y, -b + z);
    putpixelone(-b + x, -c + y, -a + z);
}
  
// detailed explanation of this algorithm is
// given in above link.
void drawsphere(int x, int y, int z, int r)
{
    int i = 0, j = 0;
    int k0 = r;
    int k = k0;
    int s0 = 0;
    int s = s0;
    int v0 = r - 1;
    int v = v0;
    int l0 = 2 * v0;
    int l = l0;
  
    // this while loop will form naive arcs
    while (i <= k) {
  
        // this while will form voxels in naive arcs
        while (j <= k) {
            if (s > v) {
                k = k - 1;
                v = v + l;
                l = l - 2;
            }
            if ((j <= k) && ((s != v) || (j != k))) {
  
                // this if, else and else if condition
                // can easily build using figure 2, 3
                if (i == 0 && j != 0)
  
                    // First naive arc pixels and each
                    // pixel is shared with two q-octants
                    putpixeledge1(i, j, k, x, y, z);
  
                // voxels shared between q1 and q2
                else if (i == j && j != k && i != 0)
                    putpixeledge2(i, j, k, x, y, z);
  
                // center voxel of c-octants
                else if (i == j && j == k)
                    putpixelsingle(i, j, k, x, y, z);
  
                // voxels shared between q1 and q6
                else if (j == k && i < k && i < j)
                    putpixeldouble(i, j, k, x, y, z);
  
                // starting voxel of q-octant
                else if (i == 0 && j == 0)
                    putpixelmiddle(i, j, k, x, y, z);
  
                // voxels inside the q-octant
                else
                    putpixelall(i, j, k, x, y, z);
            }
            s = s + 2 * j + 1;
            j = j + 1;
        }
        s0 = s0 + 4 * i + 2;
        i = i + 1;
  
        while ((s0 > v0) && (i <= k0)) {
            k0 = k0 - 1;
            v0 = v0 + l0;
            l0 = l0 - 2;
        }
        j = i;
        k = k0;
        v = v0;
        l = l0;
        s = s0;
    }
}
// Print out all the voxels of discrete sphere
void showallpoints(map > >& mymap)
{
    int count = 0;
  
    for (itr = mymap.begin(); itr != mymap.end(); ++itr) {
        list > l = itr->second;
        list >::iterator it;
        for (it = l.begin(); it != l.end(); ++it) {
            cout << itr->first << ", " << get<0>(*it)
                 << ", " << get<1>(*it) << endl;
            count += 1;
        }
    }
    cout << endl;
    cout << count << endl;
}
// Driver program
int main()
{
    int cx, cy, cz;
    cin >> cx >> cy >> cz;
    int r;
    cin >> r;
    drawsphere(cx, cy, cz, r);
    showallpoints(mymap);
}

输入：Center(0, 0, 0) 半径 3
输出：

(-3, 0, 0)(-3, 1, 1)(-3, -1, 1)(-3, 1, -1)(-3, -1, -1)
(-2, 1, 2)(-2, 2, 1)(-2, -1, 2)(-2, -2, 1)(-2, 1, -2)
(-2, 2, -1)(-2, -1, -2)(-2, -2, -1)(-1, 1, 3)(-1, 3, 1)
(-1, -1, 3)(-1, -3, 1)(-1, 1, -3)(-1, 3, -1)(-1, -1, -3)
(-1, -3, -1)(-1, 2, 2)(-1, -2, 2)(-1, 2, -2)(-1, -2, -2)
(0, 0, 3)(0, 0, -3)(0, 3, 0)(0, -3, 0)(0, 1, 3)(0, 3, 1)
(0, -1, 3)(0, 1, -3)(0, -1, -3)(0, -3, 1)(0, 3, -1)
(0, -3, -1)(0, 2, 2)(0, 2, 2)(0, -2, 2)(0, 2, -2)
(0, -2, -2)(0, -2, 2)(0, 2, -2)(0, -2, -2)(1, 3, 0)
(1, 0, 3)(1, -3, 0)(1, 3, 0)(1, -3, 0)(1, 0, -3)(1, 0, 3)
(1, 0, -3)(1, 1, 3)(1, 3, 1)(1, -1, 3)(1, -3, 1)(1, 1, -3)
(1, 3, -1)(1, -1, -3)(1, -3, -1)(1, 2, 2)(1, -2, 2)(1, 2, -2)
(1, -2, -2)(2, 2, 0)(2, 0, 2)(2, -2, 0)(2, 2, 0)(2, -2, 0)
(2, 0, -2)(2, 0, 2)(2, 0, -2)(2, 0, 2)(2, 0, 2)(2, 0, -2)
(2, 0, -2)(2, 2, 0)(2, -2, 0)(2, 2, 0)(2, -2, 0)(2, 1, 2)
(2, 2, 1)(2, -1, 2)(2, -2, 1)(2, 1, -2)(2, 2, -1)(2, -1, -2)
(2, -2, -1)(3, 0, 0)(3, 0, 1)(3, 0, 1)(3, 0, -1)(3, 0, -1)
(3, 1, 0)(3, -1, 0)(3, 1, 0)(3, -1, 0)(3, 1, 1)(3, -1, 1)
(3, 1, -1)(3, -1, -1)
102 (Total no. of points)

参考：
http://www.sciencedirect.com/science/article/pii/S0304397515010178#tl0010

如果您希望与专家一起参加现场课程，请参阅DSA 现场工作专业课程和学生竞争性编程现场课程。