8

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

故事开始：

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

领带：

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

附加信息：

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

我试图做的事情：

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

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

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

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

2-修改行进立方体算法，以使其使用使用“范围”isolevel值的作品来工作

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

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

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

提前致谢。

编辑：

通过说

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

我的意思是我想检索一组三角形。每个三角形应表示为 3Point3D秒。如果您不知道 aPoint3D是什么，它是一个struct包含 3 个双精度数（X、Y 和 Z）的 a，用于表示 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;

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

另一个编辑：

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

再编辑： TriTable.LookupTable这里EdgeTable.LookupTable有表格吗

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

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

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

4

2 回答 2

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

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

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

ps：如果您想要的只是显示，您可以通过获取与显示的每个点的光线连接的第一个点的值来使用光线跟踪（根据您的硬件和图片大小，每隔几秒钟期望一个图像。）对于这个你可以使用画线算法（检查维基百科，有一个样本，虽然是 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;
      }
    }
  } 
}
于 2016-12-08T13:24:43.450 回答
2

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

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

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

多亏了这一点，我们可以使用一些巧妙的优化 - 对于 8 个角中的每一个，您检查角是在网格内部还是外部，您使用这 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检查哪个 isolevel 在 valp1 和 valp2 值之间
  • 像在当前代码中一样进行插值，但使用选定的等值线（最小值或最大值）

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

于 2016-12-12T17:51:18.110 回答