我怎么能多边形化一个 bool[]

Question

如果有人关心，我正在使用 WPF....

故事开头：

我得到了一系列灰度图像（模型的切片）。用户输入 "range" 个灰度值来构建 3D 模型。因此，我创建了一个 3D bool 数组，以便更轻松地构建 3D 模型。这个数组代表一盒像素，表示每个像素是否应该built/generated/filled。

领带：

使用 bool[,,]，我想检索一个 List<Point3D[]>，其中每个 Point3D[] 的长度为 3，代表一个 3D 三角形 space。

附加信息：

生成的模型将被 3D 打印。在bool[,,]中，true表示物质存在，而false表示物质不存在。我需要避免立方体模型，其中每个像素都被一个立方体替换（根据我的目的，这将是无用的）。模型要尽可能平滑和准确。

我尝试做的事情：

1- 我实现了行进立方体算法，但它似乎并没有被创建来接受 "range" 个值。

2- 几天来我一直在努力尝试制作自己的算法，但部分失败了。（真的很复杂，如果想了解更多请告知）

一些我不知道如何实现的预期解决方案：

~~1- 修改行进立方体算法以使其使用 bool[,,]~~

2- 修改 Marching 立方体算法，使其使用 "range" 的 isolevel 值

3- 显示如何在 WPF 中实现适合我的目的的另一种算法（也许是光线投射算法）。

~~4- 请求我尝试实现的算法的来源，然后向我展示如何解决它。~~（它是在第一个 bool[,,] 中进行多边形化地方)

5- 一些其他神奇的解决方案。

提前致谢。

编辑：

说

Using a bool[,,], I want to retrieve a List<Point3D[]> where each Point3D[] has a length of 3 and represents a triangle in 3D space.

我的意思是我要检索一组三角形。每个三角形应表示为 3 Point3Ds。如果您不知道 Point3D 是什么，它是一个包含 3 个双精度数（X、Y 和 Z）的 struct，用于表示 3D 中的位置 space.

行进立方体算法的问题在于它有点模糊。我的意思是，你这样做明白什么？？

cubeindex = 0; 
if (grid.val[0] < isolevel) cubeindex |= 1; 
if (grid.val[1] < isolevel) cubeindex |= 2; 
if (grid.val[2] < isolevel) cubeindex |= 4; 
if (grid.val[3] < isolevel) cubeindex |= 8; 
if (grid.val[4] < isolevel) cubeindex |= 16; 
if (grid.val[5] < isolevel) cubeindex |= 32; 
if (grid.val[6] < isolevel) cubeindex |= 64; 
if (grid.val[7] < isolevel) cubeindex |= 128;

因此，我根本不知道从哪里开始修改它。

另一个编辑：

在对行进立方体算法进行一些试验后，我注意到对 bool[,,] 进行多边形化不会产生我期待的平滑 3D 模型，所以我的第一个预期解决方案和第四个预期解决方案解都出局了。经过大量研究，我知道了体绘制的Ray-Casting方法。看起来真的很酷，但我不确定它是否符合我的需要。虽然我知道它是如何工作的，但我真的不明白它是如何实现的。

再编辑： TriTable.LookupTable 和 EdgeTable.LookupTable 是存在的表 here

这是 MarchingCubes class:

public class MarchingCubes
{
    public static Point3D VertexInterp(double isolevel, Point3D p1, Point3D p2, double valp1, double valp2)
    {
        double mu;
        Point3D p = new Point3D();

        if (Math.Abs(isolevel - valp1) < 0.00001)
            return (p1);

        if (Math.Abs(isolevel - valp2) < 0.00001)
            return (p2);

        if (Math.Abs(valp1 - valp2) < 0.00001)
            return (p1);

        mu = (isolevel - valp1) / (valp2 - valp1);

        p.X = p1.X + mu * (p2.X - p1.X);
        p.Y = p1.Y + mu * (p2.Y - p1.Y);
        p.Z = p1.Z + mu * (p2.Z - p1.Z);

        return (p);
    }

    public static void Polygonise (ref List<Point3D[]> Triangles, int isolevel, GridCell gridcell)
    {
        int cubeindex = 0;
        if (grid.val[0] < isolevel) cubeindex |= 1;
        if (grid.val[1] < isolevel) cubeindex |= 2;
        if (grid.val[2] < isolevel) cubeindex |= 4;
        if (grid.val[3] < isolevel) cubeindex |= 8;
        if (grid.val[4] < isolevel) cubeindex |= 16;
        if (grid.val[5] < isolevel) cubeindex |= 32;
        if (grid.val[6] < isolevel) cubeindex |= 64;
        if (grid.val[7] < isolevel) cubeindex |= 128;

        // Cube is entirely in/out of the surface 
        if (EdgeTable.LookupTable[cubeindex] == 0)
            return;

        Point3D[] vertlist = new Point3D[12];

        // Find the vertices where the surface intersects the cube 
        if ((EdgeTable.LookupTable[cubeindex] & 1) > 0)
            vertlist[0] = VertexInterp(isolevel, grid.p[0], grid.p[1], grid.val[0], grid.val[1]);

        if ((EdgeTable.LookupTable[cubeindex] & 2) > 0)
            vertlist[1] = VertexInterp(isolevel, grid.p[1], grid.p[2], grid.val[1], grid.val[2]);

        if ((EdgeTable.LookupTable[cubeindex] & 4) > 0)
            vertlist[2] = VertexInterp(isolevel, grid.p[2], grid.p[3], grid.val[2], grid.val[3]);

        if ((EdgeTable.LookupTable[cubeindex] & 8) > 0)
            vertlist[3] = VertexInterp(isolevel, grid.p[3], grid.p[0], grid.val[3], grid.val[0]);

        if ((EdgeTable.LookupTable[cubeindex] & 16) > 0)
            vertlist[4] = VertexInterp(isolevel, grid.p[4], grid.p[5], grid.val[4], grid.val[5]);

        if ((EdgeTable.LookupTable[cubeindex] & 32) > 0)
            vertlist[5] = VertexInterp(isolevel, grid.p[5], grid.p[6], grid.val[5], grid.val[6]);

        if ((EdgeTable.LookupTable[cubeindex] & 64) > 0)
            vertlist[6] = VertexInterp(isolevel, grid.p[6], grid.p[7], grid.val[6], grid.val[7]);

        if ((EdgeTable.LookupTable[cubeindex] & 128) > 0)
            vertlist[7] = VertexInterp(isolevel, grid.p[7], grid.p[4], grid.val[7], grid.val[4]);

        if ((EdgeTable.LookupTable[cubeindex] & 256) > 0)
            vertlist[8] = VertexInterp(isolevel, grid.p[0], grid.p[4], grid.val[0], grid.val[4]);

        if ((EdgeTable.LookupTable[cubeindex] & 512) > 0)
            vertlist[9] = VertexInterp(isolevel, grid.p[1], grid.p[5], grid.val[1], grid.val[5]);

        if ((EdgeTable.LookupTable[cubeindex] & 1024) > 0)
            vertlist[10] = VertexInterp(isolevel, grid.p[2], grid.p[6], grid.val[2], grid.val[6]);

        if ((EdgeTable.LookupTable[cubeindex] & 2048) > 0)
            vertlist[11] = VertexInterp(isolevel, grid.p[3], grid.p[7], grid.val[3], grid.val[7]);

        // Create the triangle 
        for (int i = 0; TriTable.LookupTable[cubeindex, i] != -1; i += 3)
        {
            Point3D[] aTriangle = new Point3D[3];

            aTriangle[0] = vertlist[TriTable.LookupTable[cubeindex, i]];
            aTriangle[1] = vertlist[TriTable.LookupTable[cubeindex, i + 1]];
            aTriangle[2] = vertlist[TriTable.LookupTable[cubeindex, i + 2]];

            theTriangleList.Add(aTriangle);
        }
    }
}

这是 GridCell class:

public class GridCell
{
    public Point3D[] p = new Point3D[8];
    public Int32[] val = new Int32[8];
}

好的，这就是我正在做的：

List<int[][]> Slices = new List<int[][]> (); //Each slice is a 2D jagged array. This is a list of slices, each slice is a value between about -1000 to 2300

public static void Main ()
{
    List <Point3D[]> Triangles = new List<Point3D[]> ();
    int RowCount;//That is my row count
    int ColumnCount;//That is my column count
    //Here I fill the List with my values
    //Now I'll polygonise each GridCell
    for (int i = 0; i < Slices.Count - 1; i++)
    {
        int[][] Slice1 = Slices[i];
        int[][] Slice2 = Slices[i + 1];
        for (int j = 0; j < RowCount - 1; j++)
        {
            for (int k = 0; k < ColumnCount - 1; k++)
            {
                GridCell currentCell = GetCurrentCell (Slice1, Slice2, j, k);
                Polygonise (ref Triangles, int isoLevel, GridCell currentCell);
            }
        }
    }
    //I got the "Triangles" :D
}

//What this simply does is that it returns in 3D space from RI and CI
public static GridCell GetCurrentCell (int[][] CTSliceFront, int[][] CTSliceBack, int RI, int CI)//RI is RowIndex and CI is ColumnIndex
{
    //Those are preset indicating X,Y or Z coordinates of points/edges
    double X_Left_Front;
    double X_Right_Front;
    double X_Left_Back;
    double X_Right_Back;
    double Y_Top_Front;
    double Y_Botton_Front;
    double Y_Top_Back;
    double Y_Botton_Back;
    double Z_Front;
    double Z_Back;

    GridCell currentCell = new GridCell();
    currentCell.p[0] = new Point3D(X_Left_Back, Y_Botton_Back, Z_Back);
    currentCell.p[1] = new Point3D(X_Right_Back, Y_Botton_Back, Z_Back);
    currentCell.p[2] = new Point3D(X_Right_Front, Y_Botton_Front, Z_Front);
    currentCell.p[3] = new Point3D(X_Left_Front, Y_Botton_Front, Z_Front);
    currentCell.p[4] = new Point3D(X_Left_Back, Y_Top_Back, Z_Back);
    currentCell.p[5] = new Point3D(X_Right_Back, Y_Top_Back, Z_Back);
    currentCell.p[6] = new Point3D(X_Right_Front, Y_Top_Front, Z_Front);
    currentCell.p[7] = new Point3D(X_Left_Front, Y_Top_Front, Z_Front);

    currentCell.val[] = CTSliceBack[theRowIndex + 1][theColumnIndex];
    currentCell.val[] = CTSliceBack[theRowIndex + 1][theColumnIndex + 1];
    currentCell.val[] = CTSliceFront[theRowIndex + 1][theColumnIndex + 1];
    currentCell.val[] = CTSliceFront[theRowIndex][theColumnIndex];
    currentCell.val[] = CTSliceBack[theRowIndex][theColumnIndex + 1];
    currentCell.val[] = CTSliceFront[theRowIndex][theColumnIndex + 1];
    currentCell.val[] = CTSliceFront[theRowIndex][theColumnIndex];

    return currentCell;
}

我是从 Stack-Overflow 做的，所以请指出任何拼写错误。这很有效，这导致模型的 shell。我想要做的是用 2 个整数替换 Polygonise 函数中的 isolevel 变量：isolevelMin 和 isolevelMax。不只是 1 int，而是一个范围。

编辑： 伙计们，请帮忙。 50 个代表几乎是我声誉的 1/2。我希望我能改变我的期望。现在我只想知道如何实现我想要的东西，关于如何实现的任何片段，或者如果有人给出一些我可以 google 我用来解决我的问题的关键字，他会得到代表。我只需要一些光 - 这里很黑。

编辑很棒： 经过越来越多的研究，我发现体积光线行进算法实际上并不能生成表面。由于我需要打印生成的3D模型，所以我现在只有一个办法。我必须让行进立方体算法从最大和最小等值水平值生成一个空心模型。

Answer 1

免责声明：这不会是一个完整的答案，因为我现在没有时间给出一个，我不急于获得赏金。

我将开始构思一种进行光线追踪的方法，它应该会生成一个在您的形状外部可见的点列表。

基本上对于光线追踪，您从每个面的每个点（可能是对角线）追踪一条光线，这不是最好的方法，但应该会输出一些东西。

射线是一个矢量，在这种情况下，如果你从一侧的 x,y 平面到另一侧，它将是这样的 [pseudocode]

bool[,,] points = getpoints();
bool[,,] outsidepoints = new bool[width,height,depth];
//this is traceray from x,y,0 to x,y,depth
for(int x=0;x<width;x++){
  for(int y=0;y<height;y++){
    for(int z=0;z<depth;z++){
      if (points[x,y,z]==true){
        outsidepoints[x,y,z]=true
        z=depth;//we break the last for
      }
    }
  } 
}

对最后一个循环中的所有 6 种组合执行相同的操作（x=0..width、x=width..0、y=0..height 等...）

如果没有悬垂物（例如单个凸形），那应该可以很好地工作，如果这样做，它们可能会被切断。

一旦你有了外部点，你只需要一个算法将它们变成三角形（就像在每个平面上沿轴制作连接切片，然后缝合最后两个连接。）

ps：如果你想要的只是显示，你可以通过获取与显示器每个点的光线连接的第一个点的值来使用光线追踪（每几秒期望一个图像，具体取决于在你的硬件和图片尺寸上。）为此你可以使用画线算法（检查维基百科，有一个样本，虽然是二维的）最好的是你可以直接从你的图片中获取颜色（只需将阈值设置为其中像素为空（光线穿过））

另一项编辑：如果您需要任何精确度，请私信我，或在此处发表评论

使用直通光线：

bool lastwasempty=true;
for(int x=0;x<width;x++){
  for(int y=0;y<height;y++){
    for(int z=0;z<depth;z++){
      if (points[x,y,z]==true){
        if(lastwasempty){
           outsidepoints[x,y,z]=true
        }
        lastwasempty=false;
      }else{
        lastwasempty=true;
      }
    }
  } 
}

Answer 2

Marching Cube 算法的假设之一是每个立方体几何体都是一个单独的实体，不依赖于其他立方体（这不完全正确，但让我们暂时坚持这一点）。

鉴于此，您可以将整个网格拆分为立方体并分别求解每个立方体。您可能还会注意到，虽然最终的网格取决于 3D 标量场的确切值，但网格的拓扑结构仅取决于某些给定立方角是否位于网格的内部或外面.

例如： 假设您的 isolevel 设置为 0 并且您的 3D 字段有正数和负数。如果您现在将某些值从 0.1 更改为 1.0，那么您只会更改某些三角形的比例，而不是它们的数量 and/or 拓扑。另一方面，如果您将某些值从 0.1 更改为 -0.1，那么您的三角形可能会切换到另一种排列，并且它们的数量可能会改变。

多亏了这一点，我们可以使用一些巧妙的优化 - 对于 8 个角中的每一个，您检查角是 inside 还是 outside一个网格，您使用这 8 位来构建一个数字 - 它成为您预先计算的 table 可能立方体解决方案的索引，我们称它们为 Variants.

每个 Cube Variant 都是针对此特定角配置的三角形列表。 Marching Cubes 中的三角形总是跨越立方体的边，所以我们使用边索引来描述它。这些索引是您在 Paul Bourke 代码中的 'magic' table 中看到的确切数据 :)。

但这还不是最终的网格。你在这里得到的三角形只是基于一个事实，如果某个角在网格的内部或外部，但它离表面有多远？对结果有影响吗？再次出现您的网格值 - 当您选择了 Variant 时，获得了三角形列表和每个三角形的 3 条边，然后您获得了边末端的值并将它们插值以找到 0 值（我们将其设置为 isolevel，还记得吗？）。这些最终的插值顶点应该用于您的网格。

这就变成了一个很长的介绍，我希望你能理解我写的内容。

我的建议： 最简单的改变是角，而不是这样做：

if (grid.val[0] < isolevel) cubeindex |= 1;

试试这样改：

if (grid.val[0] > isolevelMin && grid.val[0] < isolevelMax) cubeindex |= 1;
// And the same for other lines...

最后的插值有点棘手，我没有办法测试它，但让我们试试：

修改您的 VertexInterp 以通过两个等值线 - 最小值和最大值
在 VertexInterp 中检查哪个 isolevel 在 valp1 和 valp2 值之间
按照当前代码进行插值，但使用选定的等值线（最小值或最大值）

让我知道它是否有效 :) 或者如果您有任何问题。

我怎么能多边形化一个 bool[]

How could I polygonise a bool[,,]

c#

3d

wpf

marching-cubes

volume-rendering