Mathematica:3D线框图

9

Mathematica是否支持线框图的隐藏线消除?如果不支持,这里有没有人知道如何实现?让我们从这个问题开始:

Plot3D[Sin[x+y^2], {x, -3, 3}, {y, -2, 2}, Boxed -> False]

output

创建线框图的方法如下:
Plot3D[Sin[x+y^2], {x, -3, 3}, {y, -2, 2}, Boxed -> False, PlotStyle -> None]

output

我们可以做的一件事是将所有表面都涂成白色,以达到效果。然而,这样做是不可取的。原因是,如果我们将这个隐藏线框架模型导出为pdf,我们将拥有Mathematica用于渲染图像的所有那些白色多边形。我希望能够获得一个pdf和/或eps格式的隐藏线框架。

更新:

我已经发布了解决这个问题的方案。问题在于代码运行非常缓慢。在当前状态下,它无法为这个问题中的图像生成线框图。随意尝试我的代码。我在帖子末尾添加了一个链接。你也可以在这里找到代码link


2
您可以像您建议的那样使用白色表面(以及照明->“中性”),然后栅格化到所需的分辨率。如果您想要矢量图形表示,则我怀疑您的问题的答案是否定的。 - Andrew Moylan
@Mr.Wizard,我正在研究一个想法:通过使用ExportStringImportStringGraphics3D更改为Graphics。在两条线段相交的位置拆分所有线段。删除位于多边形内部的所有线段。最后,删除所有剩余的多边形。我目前正在尝试弄清楚如何在交点处拆分线段。一旦我制定出问题,我可能会发布迷你问题。这听起来怎么样? - jmlopez
这听起来很“有趣”。 :-) 但它也听起来很复杂。你是建议完全使用表面,还是只使用网格线?在这个上下文中,“多边形”是什么意思? - Mr.Wizard
@Mr.Wizard,不确定您是否已经查看了我的新问题“Mathematica:Joining line segments”。在那里,我展示了如何先转换为图形。如果您在该Graphics对象上使用InputForm,您将看到PolygonJoinedCurve对象。我们的想法是在它们相交的地方断开线段,从而生成更多的线段。现在我们所要做的就是检查一个线段是否在Polygon对象内,我们已经拥有了所有这些信息,因此我们可以找到一种算法来确定它是否确实存在。然后我们删除这个多边形和位于多边形内的线段。 - jmlopez
@Mr.Wizard,我已经写了一个解决方案。一定要检查代码。也许你可以改进它。否则我将不得不花更多时间找出如何在Mathematica中编写C扩展。 - jmlopez
显示剩余4条评论
2个回答

5

我这里提供一个解决方案。首先,我会展示如何使用生成线框的函数,然后我会详细解释组成算法的其他函数。


wireFrame

wireFrame[g_] := Module[{figInfo, opt, pts},
   {figInfo, opt} = G3ToG2Info[g];
   pts = getHiddenLines[figInfo];
   Graphics[Map[setPoints[#] &, getFrame[figInfo, pts]], opt]
]

这个函数的输入是一个Graphics3D对象,最好没有坐标轴。
fig = ListPlot3D[
   {{0, -1, 0}, {0, 1, 0}, {-1, 0, 1}, {1, 0, 1}, {-1, 1, 1}},
   Mesh -> {10, 10},
   Boxed -> False,
   Axes -> False,
   ViewPoint -> {2, -2, 1},
   ViewVertical -> {0, 0, 1},
   MeshStyle -> Directive[RGBColor[0, 0.5, 0, 0.5]],
   BoundaryStyle -> Directive[RGBColor[1, 0.5, 0, 0.5]]
]

surface

现在我们使用函数wireFrame

wireFrame[fig]

线框图

如您所见,wireFrame获得了大部分线条和颜色。有一条绿线没有包含在线框图中。这很可能是由于我的阈值设置。

在我继续解释函数G3ToG2InfogetHiddenLinesgetFramesetPoints的详细内容之前,我将向您展示为什么带有隐藏线条的线框图可能会有用。

光栅线框图

上面显示的图像是使用三维图形中的光栅技术结合此处生成的线框图生成的pdf文件的屏幕截图。这可以有多种优势。无需保留三角形信息以显示彩色表面。相反,我们显示表面的光栅图像。所有的线都非常平滑,除了光栅图未被线覆盖的边界之外。我们还减少了文件大小。在这种情况下,使用光栅图和线框图的组合后,pdf文件大小从1.9mb减少到78kb。它在pdf查看器中显示的时间更短,图像质量也很好。

Mathematica在将3D图像导出为pdf文件方面做得非常好。当我们导入pdf文件时,我们获得由线段和三角形组成的图形对象。在某些情况下,这些对象重叠,因此我们有隐藏的线条。要制作一个没有表面的线框模型,我们首先需要删除这个重叠部分,然后删除多边形。我将从描述如何从Graphics3D图像中获取信息开始。


G3ToG2Info

getPoints[obj_] := Switch[Head[obj], 
   Polygon, obj[[1]], 
   JoinedCurve, obj[[2]][[1]], 
   RGBColor, {Table[obj[[i]], {i, 1, 3}]}
  ];
setPoints[obj_] := Switch[Length@obj, 
   3, Polygon[obj], 
   2, Line[obj], 
   1, RGBColor[obj[[1]]]
  ];
G3ToG2Info[g_] := Module[{obj, opt},
   obj = ImportString[ExportString[g, "PDF", Background -> None], "PDF"][[1]];
   opt = Options[obj];
   obj = Flatten[First[obj /. Style[expr_, opts___] :> {opts, expr}], 2];
   obj = Cases[obj, _Polygon | _JoinedCurve | _RGBColor, Infinity];
   obj = Map[getPoints[#] &, obj];
   {obj, opt}
  ]

这段代码适用于Mathematica 8,如果版本为7,则需要在函数getPoints中将JoinedCurve替换为Line。函数getPoints假定您提供的是原始Graphics对象。它会查看接收到的对象类型,然后从中提取所需的信息。如果是多边形,则获取3个点的列表;对于直线,则获得2个点的列表;如果是颜色,则获取包含3个点的单个列表的列表。这样做是为了保持与列表的一致性。
函数setPoints执行与getPoints相反的操作。您输入一个点列表,它将确定是否返回多边形、直线或颜色。
要获取三角形、线和颜色的列表,我们使用G3ToG2Info。此函数将使用ExportStringImportStringGraphics3D版本中获取Graphics对象。此信息存储在obj中。我们需要执行一些清理工作,首先获取obj的选项。这部分是必要的,因为它可能包含图像的PlotRange。然后,我们按照obtaining graphics primitives and directives中描述的方式获取所有PolygonJoinedCurveRGBColor对象。最后,我们对所有这些对象应用函数getPoints,以获取三角形、线和颜色的列表。此部分涵盖了行{figInfo, opt} = G3ToG2Info[g]

getHiddenLines

我们希望能够知道哪些线段的一部分不会被显示。为此,我们需要知道两个线段之间的相交点。我使用的算法可以在这里找到。
lineInt[L_, M_, EPS_: 10^-6] := Module[
  {x21, y21, x43, y43, x13, y13, numL, numM, den},
  {x21, y21} = L[[2]] - L[[1]];
  {x43, y43} = M[[2]] - M[[1]];
  {x13, y13} = L[[1]] - M[[1]];
  den = y43*x21 - x43*y21;
  If[den*den < EPS, Return[-Infinity]];
  numL = (x43*y13 - y43*x13)/den;
  numM = (x21*y13 - y21*x13)/den;
  If[numM < 0 || numM > 1, Return[-Infinity], Return[numL]];
 ]

lineInt 假设线段 LM 不重合。如果这两条直线平行,或者包含线段 L 的直线没有和线段 M 相交,则返回 -Infinity。如果包含线段 L 的直线与线段 M 相交,则返回一个标量值。假设这个标量值为 u,则相交点为L[[1]] + u (L[[2]]-L[[1]])。需要注意的是,u 可以是任意实数。你可以使用此 manipulate 函数来测试 lineInt 的工作原理。

Manipulate[
   Grid[{{
      Graphics[{
        Line[{p1, p2}, VertexColors -> {Red, Red}],
        Line[{p3, p4}]
       },
       PlotRange -> 3, Axes -> True],
      lineInt[{p1, p2}, {p3, p4}]
     }}],
   {{p1, {-1, 1}}, Locator, Appearance -> "L1"},
   {{p2, {2, 1}}, Locator, Appearance -> "L2"},
   {{p3, {1, -1}}, Locator, Appearance -> "M1"},
   {{p4, {1, 2}}, Locator, Appearance -> "M2"}
]

示例

现在我们知道从L[[1]]到线段M需要走多远,我们可以找出线段的哪一部分位于三角形内。


注:此处的代码片段中包含HTML标记,请勿删除。
lineInTri[L_, T_] := Module[{res},
  If[Length@DeleteDuplicates[Flatten[{T, L}, 1], SquaredEuclideanDistance[#1, #2] < 10^-6 &] == 3, Return[{}]];
  res = Sort[Map[lineInt[L, #] &, {{T[[1]], T[[2]]}, {T[[2]], T[[3]]},  {T[[3]], T[[1]]} }]];
  If[res[[3]] == Infinity || res == {-Infinity, -Infinity, -Infinity}, Return[{}]];
  res = DeleteDuplicates[Cases[res, _Real | _Integer | _Rational], Chop[#1 - #2] == 0 &];
  If[Length@res == 1, Return[{}]];
  If[(Chop[res[[1]]] == 0 && res[[2]] > 1) || (Chop[res[[2]] - 1] == 0 && res[[1]] < 0), Return[{0, 1}]];
  If[(Chop[res[[2]]] == 0 && res[[1]] < 0) || (Chop[res[[1]] - 1] == 0 && res[[2]] > 1), Return[{}]];
  res = {Max[res[[1]], 0], Min[res[[2]], 1]};
  If[res[[1]] > 1 || res[[1]] < 0 || res[[2]] > 1 || res[[2]] < 0, Return[{}], Return[res]];
 ]

该函数返回需要删除的行L的部分。例如,如果它返回{.5, 1},这意味着您将删除该行的50%,从半段开始到段的结束点。如果L = {A, B},并且该函数返回{u, v},则表示线段{A + (B-A)u,A + (B-A)v}是包含在三角形T中的线段。
实现lineInTri时,需要注意线L不是T的边之一,如果是这种情况,则该线不在三角形内。这就是舍入误差可能会产生负面影响的地方。当Mathematica导出图像时,有时一条线位于三角形的边缘,但这些坐标略有不同。我们需要决定线与边缘的距离多近才算在边缘上,否则该函数将看到该线几乎完全位于三角形内。这就是函数中第一行的原因。要查看线是否位于三角形的边缘上,我们可以列出三角形和线的所有点,并删除所有重复项。在这种情况下,您需要指定什么是重复项。最后,如果我们得到一个由3个点组成的列表,则表示该线位于边缘上。接下来的部分有点复杂。我们检查线L与三角形T的每条边是否相交,并将结果存储在列表中。然后,我们对列表进行排序,并找出线段在三角形内的哪个部分(如果有)。通过尝试使用一些测试,例如检查线端点是否是三角形的顶点,线是否完全位于三角形内部、部分内部或完全外部,可以理解它的意义。
Manipulate[
  Grid[{{
    Graphics[{
      RGBColor[0, .5, 0, .5], Polygon[{p3, p4, p5}],
      Line[{p1, p2}, VertexColors -> {Red, Red}]
     },
     PlotRange -> 3, Axes -> True],
    lineInTri[{p1, p2}, {p3, p4, p5}]
   }}],
 {{p1, {-1, -2}}, Locator, Appearance -> "L1"},
 {{p2, {0, 0}}, Locator, Appearance -> "L2"},
 {{p3, {-2, -2}}, Locator, Appearance -> "T1"},
 {{p4, {2, -2}}, Locator, Appearance -> "T2"},
 {{p5, {-1, 1}}, Locator, Appearance -> "T3"}
]

三角形测试

lineInTri 将用于查看线段的哪一部分不会被绘制。这条线段很可能会被许多三角形覆盖。因此,我们需要保留所有每条线段不会被绘制的部分的列表。这些列表没有顺序。我们只知道这些列表是一维的线段。每个线段都由 [0,1] 区间内的数字组成。我不知道有没有针对一维线段的并集函数,所以这里是我的实现。

union[obj_] := Module[{p, tmp, dummy, newp, EPS = 10^-3},
  p = Sort[obj];
  tmp = p[[1]];
  If[tmp[[1]] < EPS, tmp[[1]] = 0];
  {dummy, newp} = Reap[
    Do[
     If[(p[[i, 1]] - tmp[[2]]) > EPS && (tmp[[2]] - tmp[[1]]) > EPS, 
       Sow[tmp]; tmp = p[[i]], 
       tmp[[2]] = Max[p[[i, 2]], tmp[[2]]]
      ];
     , {i, 2, Length@p}
    ];
    If[1 - tmp[[2]] < EPS, tmp[[2]] = 1];
    If[(tmp[[2]] - tmp[[1]]) > EPS, Sow[tmp]];
   ];
  If[Length@newp == 0, {}, newp[[1]]]
 ]

该函数本来可以更简短,但是这里包含了一些if语句来检查一个数字是否接近于零或一。如果一个数字与零相差小于EPS,则将该数字变为零,对于数字一同样适用。我在这里还要涵盖的另一个方面是,如果需要显示的线段比较小,则很可能需要将其删除。例如,如果我们有{{0,.5}, {.500000000001}},这意味着我们需要绘制{{.5, .500000000001}}。但是这个线段非常小,在大的线段中甚至不会被注意到,我们不知道这两个数字是否相同。在实现union时,所有这些事情都需要考虑到。

现在我们准备好看看需要从线段中删除什么了。下一个步骤需要使用G3ToG2Info生成的对象列表、该列表中的一个对象和一个索引。

getSections[L_, obj_, start_ ] := Module[{dummy, p, seg},
  {dummy, p} = Reap[
    Do[
     If[Length@obj[[i]] == 3,
      seg =  lineInTri[L, obj[[i]]];
      If[Length@seg != 0, Sow[seg]];
     ]
     , {i, start, Length@obj}
    ]
   ];
  If[Length@p == 0, Return[{}], Return[union[First@p]]];
 ]

getSections 返回一个包含需要从 L 中删除的部分的列表。我们知道 obj 是三角形、线条和颜色的列表,我们知道列表中索引更高的对象将绘制在比索引更低的对象上面。因此,我们需要索引 start。这是我们将从 obj 中开始查找三角形的索引。一旦我们找到一个三角形,我们将使用函数 lineInTri 获取位于三角形内的线段的部分。最后,我们将得到一个部分列表,可以通过使用 union 组合。

最后,我们来到 getHiddenLines。所有这些需要做的就是查看由 G3ToG2Info 返回的列表中的每个对象,并应用函数 getSectionsgetHiddenLines 将返回一个列表的列表。每个元素都是需要删除的部分的列表。

getHiddenLines[obj_] := Module[{pts},
  pts = Table[{}, {Length@obj}];
  Do[
   If[Length@obj[[j]] == 2,
      pts[[j]] = getSections[obj[[j]], obj, j + 1]
    ];
    , {j, Length@obj}
   ];
   Return[pts];
  ]

getFrame

如果你已经理解了前面的概念,我相信你知道接下来要做什么。如果我们有三角形、线段和颜色的列表,以及需要删除的线段部分,我们只需要绘制可见的颜色和线段部分。首先,我们创建一个complement函数,这将告诉我们到底应该绘制什么。

complement[obj_] := Module[{dummy, p},
  {dummy, p} = Reap[
    If[obj[[1, 1]] != 0, Sow[{0, obj[[1, 1]]}]];
    Do[
     Sow[{obj[[i - 1, 2]], obj[[i, 1]]}]
     , {i, 2, Length@obj}
    ];
    If[obj[[-1, 2]] != 1, Sow[{obj[[-1, 2]], 1}]];
   ];
  If[Length@p == 0, {}, Flatten@ First@p]
 ]

现在,getFrame函数。
getFrame[obj_, pts_] := Module[{dummy, lines, L, u, d},
  {dummy, lines} = Reap[
    Do[
     L = obj[[i]];
     If[Length@L == 2,
      If[Length@pts[[i]] == 0, Sow[L]; Continue[]];
      u = complement[pts[[i]]];
      If[Length@u > 0, 
       Do[
        d = L[[2]] - L[[1]];
        Sow[{L[[1]] + u[[j - 1]] d, L[[1]] + u[[j]] d}]
        , {j, 2, Length@u, 2 }]
      ];
    ];
    If[Length@L == 1, Sow[L]];
    , {i, Length@obj}]
  ];
 First@lines
]

结束语

我对算法的结果有些满意,但我不喜欢执行速度。我用循环编写了这个算法,就像在C/C++/Java中一样。我尽力使用ReapSow创建增长列表,而不是使用函数Append。尽管如此,我仍然必须使用循环。需要注意的是,发布在这里的框架图片生成需要63秒。我尝试为问题中的图片制作线框图,但该3D对象包含约32000个对象。计算每行需要显示的部分大约需要13秒钟。如果我们假设有32000条线,并且需要13秒钟来完成所有计算,那么将需要大约116小时的计算时间。

我相信如果我们在所有例程上使用函数Compile,并找到不使用Do循环的方法,这个时间可以缩短。Stack Overflow能帮忙吗?

为了方便起见,我已经将代码上传到了网上。你可以在这里找到它:here。如果你能将这个代码的修改版本应用于问题中的图形,并显示出线框图,我会把你的解决方案标记为这篇文章的答案。

最好的祝愿, J Manuel Lopez


1
+1. 或许 codereview.stackexchange.com 是一个关于代码优化问题的正确场所。 - Alexey Popkov
谢谢Alexey,我会查看它。 - jmlopez
@jmlopez 可能最好让 Mathematica 带着 CompilationTarget->C 的选项为你完成它。 - Alexey Popkov
@Alexey,那将是不错的,但我只是在尝试编译一个函数时遇到了困难。你会如何编译函数“lineInt”? - jmlopez
谢谢Alexey,这应该让我开始看看我能在M8上跑多快。 - jmlopez
显示剩余2条评论

1

这不是正确的,但有些有趣:

Plot3D[Sin[x + y^2], {x, -3, 3}, {y, -2, 2}, Boxed->False, PlotStyle->{EdgeForm[None], FaceForm[Red, None]}, Mesh->False]

使用FaceForm为None时,多边形不会被渲染。我不确定是否有办法对网格线进行操作。


网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接