我怎么能多边形化 bool[,,]答案

【问题标题】：How could I polygonise a bool[,,]我怎么能多边形化 bool[,,]
【发布时间】：2016-12-01 17:30:06
【问题描述】：

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

故事开始：

我得到了一系列灰度图像（模型的切片）。用户输入灰度值的“范围”来构建 3D 模型。因此，我创建了一个 3D bool 数组，以便更轻松地构建 3D 模型。这个数组代表一个像素框，表示每个像素是否应该被构建/生成/填充。

领带：

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

其他信息：

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

我想做什么：

1- 我实现了行进立方体算法，但它似乎并没有被创建来接受值的“范围”。

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

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

~~1- 修改 Marching cubes 算法以使其使用 bool[,,]~~

2- 修改 Marching cubes 算法以使其使用使用 isolevel 值的“范围”的作品来工作

3- 展示如何在 WPF 中实现另一种适合我的目的的算法（可能是 Ray-Casting 算法）。

~~4- 请求我尝试实现的算法的源代码，然后向我展示如何解决它。~~（首先将bool[,,] 多边形化）

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

提前致谢。

编辑：

说出来

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

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

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

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 类：

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 类：

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 做的，所以请指出任何错别字。这是有效的，这会产生一个模型的外壳。我想做的是将Polygonise 函数中的isolevel 变量替换为2 个整数：isolevelMin 和isolevelMax。不仅是 1 个int，而是一个范围。

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

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

【问题讨论】：

“使用 bool[,,]，我想检索一个 List，其中每个 Point3D[] 的长度为 3，表示 3D 空间中的一个三角形”：1. On一个便签本，用手写的准确写出你的意思，你的猜测肯定比我的好很多。 2. 对其进行细化，直到它看起来像伪代码。 3. 简单的部分：将其翻译成真正的代码。 4. 如果遇到问题，请寻求帮助。
“我实现了行进立方体算法，但它似乎并没有被创建来接受一个“范围”的值”：当一个问题说“似乎没有”时，我死了小的。该句子中的信息四舍五入为零。如果您在实现该算法时遇到具体问题，请发布您的代码并说出问题所在。但是您不是在描述具体的实现问题；你在描述一种情绪状态。也许这里有人可以开处方，我不知道。我不能。
here 是一篇可能有助于讨论在使用行进立方体算法时如何平滑网格的论文
@投票接近的人 - 请不要关闭这个。提问者是一个新人。他显然已经在他的研究中付出了努力，并写了一个格式良好的问题，并且仅仅因为没有足够的代码（我确信这是因为复杂性而不是无法编码）它获得了接近的选票。 Ed Plunkett 是对的，我们需要准确的信息，但我相信 OP 很乐意提供它。请至少在你投票之前让赏金时间消逝。
@None 人们在有准确的代码要处理之前不会发布答案。您说您已经看到了您不理解的伪代码，将其链接到您处理过的其他源材料（即使它不起作用）。设置一个 github minimal complete sample 项目，其中包含您当前状态的所有代码并将其链接到此处，以便我们可以准确地看到您所做的以及需要修复的内容。就目前而言，您的赏金代表将被浪费。

标签： c# wpf 3d volume-rendering marching-cubes

【解决方案1】：

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

鉴于此，您可以将整个网格分割成立方体并分别求解。您可能还注意到，尽管最终网格取决于 3D 标量场的精确值，但网格的拓扑结构仅取决于某个给定立体角是否在网格内部或外面。

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

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

每个Cube Variant 都是一个三角形列表，用于这种特定的角配置。 Marching Cubes 中的三角形总是跨越立方体边缘，因此我们使用边缘索引来描述它。这些索引是您在 Paul Bourke 代码中的那些“魔术”表中看到的确切数据:)。

但这还不是最终的网格。你在这里得到的三角形只是基于一个事实，如果某个角在网格内部或外部，但它离表面有多远？对结果有影响吗？再次出现您的网格值 - 当您选择 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 内部检查哪个等值水平位于 valp1 和 valp2 值之间
按照当前代码进行插值，但使用选定的等值线（最小值或最大值）

如果它有效，请告诉我 :) 或者如果您有任何问题。

【讨论】：

【解决方案2】：

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

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

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

射线是一个向量，在这种情况下，如果你从一侧的 x,y 平面到另一侧，它会是这样的 [伪代码]

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 等...）

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

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

ps：如果您想要的只是显示，您可以通过获取与显示的每个点的光线连接的第一个点的值来使用光线跟踪（根据您的硬件和图片大小，预计每隔几秒钟出现一次图像.) 为此，您可以使用画线算法（查看 wikipedia，有一个示例，尽管是 2d）最好的事情是您可以直接从图片中获取颜色（只需设置像素无效的阈值（光线穿过））

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

使用穿透光线：

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

【讨论】：

问题是我需要根据最小预设值使图像是空心的。如果您向我展示如何获取球体内部边缘上的点，而不仅仅是外边界，我会给 +1。感谢“画线算法”。
@None 你可以尝试让光线通过（不要在第一次命中时中断），然后检查之前的点是否透明（假），我会为此添加一些代码在答案中。
@None ，当您解析所有侧面（以及顶部和底部）时，代码将检测之后是否有空格，然后一种方式在另一种方式之后。要创建可打印的模型，您可能只需逐行扫描或制作一个以每个点为中心的立方体，但这可能会使某些东西变得难看。所以我想，如果你需要一个更平滑的模型，你可以使用 Catmull–Clark 细分曲面算法。要制作 3d 可打印文件，您可以将其导出到 obj hodge.net.au/sam/blog/wp-content/uploads/obj_format.txt 并使用切片器，或者只生成 g 代码。