tak0kadaの何でもノート

発声練習、生存確認用。

医学関連は 医学ノート

2値、2次元でのMarching Cubes実装

Polygonising a scalar field (Marching Cubes)http://users.polytech.unice.fr/~lingrand/MarchingCubes/algo.htmlを見つつ実装した。ただし2つめのページにあるように表面に穴が開く場合がある。

この例では4x4のマス目に値がある場合について計算している。実際にはimgのところを適当なデータに置き換える。1になっている領域が連結であるかのチェックも本来は必要。

#include <cstdint>
#include <array>
#include <vector>
#include <iostream>

int main(void)
{
    // type definition
    using Point = std::array<int, 2>;

    // stores 2^4 pattern of borders(if none: -1)
    const std::vector<std::vector<int>> vtable = {
        {-1, -1, -1, -1}, {1, 4, -1, -1}, {1, 2, -1, -1}, {2, 4, -1, -1},
        {2, 3, -1, -1}, {1, 2, 3, 4}, {1, 3, 0, 0}, {3, 4, -1, -1},
        {3, 4, -1, -1}, {1, 3, 0, 0}, {1, 4, 2, 3}, {2, 3, -1, -1},
        {2, 4, -1, -1}, {1, 2, -1, -1}, {1, 4, -1, -1}, {-1, -1, -1, -1}
    };


    int W=4, H=4;
    std::vector<std::vector<double>> img =
    {{0,0,0,0},
     {0,1,1,0},
     {0,1,0,0},
     {0,0,0,0}};

    std::vector<std::array<Point, 2>> edges;

    double p1, p2, p3, p4;
    for (int w = 0; w < W-1; ++w)
    {
        for (int h = 0; h < H-1; ++h)
        {
            p1 = img[h][w]; p2 = img[h][w+1];
            p3 = img[h+1][w+1]; p4 = img[h+1][w];
            std::int8_t index = 0;

            if (p1 == 1)
                index |= 1 << 0;
            else if (p2 == 1)
                index |= 1 << 1;
            else if (p3 == 1)
                index |= 1 << 2;
            else if (p4 == 1)
                index |= 1 << 3;

            for (int i = 0; vtable[index][i] != -1; i += 2)
            {
                std::array<Point, 2> edge{};
                if (vtable[index][i] == 1)
                    edge[0] = Point{{2*w+1, 2*h}};
                else if (vtable[index][i] == 2)
                    edge[0] = Point{{2*w+2, 2*h+1}};
                else if (vtable[index][i] == 3)
                    edge[0] = Point{{2*w+1, 2*h+2}};
                else if (vtable[index][i] == 4)
                    edge[0] = Point{{2*w, 2*h+1}};

                if (vtable[index][i+1] == 1)
                    edge[1] = Point{{2*w+1, 2*h}};
                else if (vtable[index][i+1] == 2)
                    edge[1] = Point{{2*w+2, 2*h+1}};
                else if (vtable[index][i+1] == 3)
                    edge[1] = Point{{2*w+1, 2*h+2}};
                else if (vtable[index][i+1] == 4)
                    edge[1] = Point{{2*w, 2*h+1}};

                edges.push_back(edge);
            }
        }
    }

    for (auto edge: edges)
    {
        std::cout << "(" << edge[0][0] << ", " << edge[0][1] << "), "
                  << "(" << edge[1][0] << ", " << edge[1][1] << ")\n";
    }

    return 0;
}

出力は以下のようになる

(2, 1), (1, 2)
(1, 2), (2, 3)
(1, 4), (2, 5)
(4, 1), (3, 2)
(3, 2), (2, 3)
(3, 4), (2, 5)
(5, 2), (4, 1)
(5, 2), (4, 3)