将海量数据握在手中——使用 ENVI LiDAR 和 IDL 从点云到 3D 打印模型
20235 为本文评分：
4.8
将海量数据握在手中——使用 ENVI LiDAR 和 IDL 从点云到 3D 打印模型
Jim Pendleton 2017年5月2日，星期二
对我们中的一些人来说，将一个物体握在手中，比在屏幕上甚至虚拟现实环境中观察它，能更直观地理解其在三维空间中的形态。
土地利用或军事规划者、建筑师和土木工程师有时会依赖其环境比例模型，以更好地理解地形特征及其决策的影响。
本篇博客将介绍一种技术，使用 ENVI LiDAR 以及一些自定义 IDL 代码，将公开可用的 LiDAR 数据转换为可 3D 打印的物体——一个比例模型。
本研究的主体是风景如画的 Mt. Princeton 的一部分，这座海拔 14,196 英尺的山峰位于科罗拉多州 Buena Vista 镇附近的阿肯色河谷西侧。以下是该地区的 Google Earth 视图，以便获得整体印象。
该山脉部分区域的 LiDAR 数据可以通过 美国地质调查局的地球探测器在线门户 公开获取。注册一个免费账户很简单。
登录后，使用地图显示导航到感兴趣的区域。
最初的目标是确定普林斯顿山山顶周围区域有哪些可用的数据集。在本例中，我手动绘制了一个红色边框。点击“Data Sets（数据集）”按钮将进入下一页，您可以从众多数据源中进行选择。
为了生成可打印的表面，我们将从与 LiDAR 采集相关的点云数据开始。查看此感兴趣区域内可用数据的清单。
对于每个数据集，足迹图标将在图像显示中高亮显示相关数据的地面投影。在本例中，我选择了数据集编号 8，它显示了相当大的高程变化。
点击数据集的“download（下载）”按钮，将弹出一个对话框让您选择输出数据的格式。
ENVI LiDAR 可以读取压缩和未压缩的数据文件格式，因此我选择了“LAZ Product（压缩 LiDAR）”选项以加快传输速度。
将生成的文件下载到本地计算机。运行 ENVI 5.4 并选择下载的“.laz”文件作为输入。由于点云数据在此格式中不附带地图投影元数据，系统会首先提示您为点云数据选择坐标系。
我知道此数据集是在 UTM 13N 区坐标的 WGS84 坐标系下采集的。实际上，就我们打印数据的目的而言，坐标系并不是特别重要。除非您拥有一台打印体积非常庞大的打印机，否则我们最终将重新缩放和重新定位单位，从米级缩放到几厘米级。
该区域的原始数据在地面上大约为 1500 x 400 米，海拔范围为 500 米。
应用投影后，点数据将显示在 ENVI LiDAR 用户界面中。最初，这将包括从树木、岩石、徒步旅行者和土拨鼠等物体返回的点。
点云不能直接打印，因为点不占据物理体积。在本例中，我们将首先创建一个代表表面的多边形网格。表面本身也不足以进行 3D 打印，因为表面没有物理深度。但我们将在后续步骤中增加第三个维度。
为了创建初始表面模型，我们将使用 ENVI LiDAR 中的 Process Data（处理数据）选项。
在本例中，我们将以不规则三角网（TIN）的形式创建地形表面模型。TIN 将为我们打印物理物体的顶部提供所需的顶点和连接性。
默认情况下，顶点和连接性将写入一组 ShapeFile 格式的文件。处理过程中的进度显示在 ENVI LiDAR 的日志窗口中。
关闭 ENVI LiDAR 用户界面，返回到主 ENVI 用户界面。选择刚刚创建的 ShapeFile 并显示它。
在 ENVI 中，此视图可以作为图像导出到 Google Earth，并可以覆盖在地表之上。这是一个相关的导览文件。
为了将此表面模型输入打印机，我们必须首先从 ShapeFile 中提取数据。由 ENVI LiDAR 生成的 ShapeFile 仅将三角剖分存储为每个三角形一个独立实体。因此，在多个三角形共享一个顶点的地方存在大量重复。在后续步骤中，定位边界顶点时，我们需要考虑这种共享连接性。
以下 IDL 例程接受 .shp 文件的路径作为输入，并生成顶点列表和连接性作为输出，同时移除重复顶点。欢迎改进其效率。
pro shptopoly, file, v, c os = idlffshape(file) os.getproperty, n_entities = nen c = ulonarr(nen*4) tol = .005^2 ent = os.getentity(0) vc = *ent.vertices c[0] = [3L, 0, 1, 2] xs = reform(vc[0, 0:2]) ys = reform(vc[1, 0:2]) zs = reform(vc[2, 0:2]) vcount = 3LL for x = 1l, nen - 1 do begin ent = os.getentity(x) vc = *ent.vertices c[x*4] = 3 for i = 0l, 2 do begin r2 = min((vc[0, i] - xs)^2 + (vc[1, i] - ys)^2 + (vc[2, i] - zs)^2, mr2) if (r2 gt tol) then begin c[x*4 + i + 1] = vcount xs = [xs, reform(vc[0, i])] ys = [ys, reform(vc[1, i])] zs = [zs, reform(vc[2, i])] vcount++ endif else begin c[x*4 + i + 1] = mr2 endelse endfor endfor v = transpose([[xs], [ys], [zs]]) end
我们可以使用此例程将新模型与 EarthExplorer 中的感兴趣区域进行比较。
IDL> shptopoly, pathto + 'terrainTIN.shp', v, c
IDL> xobjview, idlgrpolygon(v, poly = c, color = !color.gold, style = 2)
由于表面本身没有深度，因此无法打印。增加深度的一种方法是在边界周围创建一个“裙边”。首先，我们必须定位外部顶点。外部顶点是其相连三角形的角度之和小于 360 度的顶点。
以下例程返回外部顶点的索引。同样欢迎优化此例程。
function exteriorvertices, v, c compile_opt strictarr on_error, 2 indices = [] c0 = c[1:*:4] c1 = c[2:*:4] c2 = c[3:*:4] tol = 10./180*!dpi twopi = !dpi*2. for i = 0L, v.length/3 - 1 do begin hasvertex = where(c0 eq i or c1 eq i or c2 eq i, n) if (n eq 0) then continue allthetas = 0.d for j = 0L, n - 1 do begin v0 = double((v[*, c0[hasvertex[j]]])[0:1]) v1 = double((v[*, c1[hasvertex[j]]])[0:1]) v2 = double((v[*, c2[hasvertex[j]]])[0:1]) h = hasvertex[j] case 1 of c0[hasvertex[j]] eq i : begin aa = sqrt(total((v2 - v1)^2)) bb = sqrt(total((v0 - v1)^2)) cc = sqrt(total((v2 - v0)^2)) end c1[hasvertex[j]] eq i : begin aa = sqrt(total((v2 - v0)^2)) bb = sqrt(total((v1 - v0)^2)) cc = sqrt(total((v2 - v1)^2)) end c2[hasvertex[j]] eq i : begin aa = sqrt(total((v1 - v0)^2)) bb = sqrt(total((v2 - v1)^2)) cc = sqrt(total((v2 - v0)^2)) end else : endcase theta = acos((bb^2 + cc^2 - aa^2)/(2.*bb*cc)) allthetas += abs(theta) endfor if (abs(twopi - allthetas) gt tol) then indices = [indices, i] endfor return, indices end
将之前创建的顶点和连接性列表作为输入传递。索引向量是函数的返回值。
IDL> e = exteriorvertices(v, c)
（代码中未解决一些边界情况，例如被 TIN 过程计算为与 XY 平面垂直的边界三角形。您会如何处理这些？）
有了外部索引，我们现在可以通过沿 Z 维度向下投影一组新的顶点到任意基底深度（低于原始数据中最低 Z 值），为模型生成多边形网格裙边。然后，我们用三角形将上部和下部的外部多边形缝合以提供深度，并对基底进行三角剖分以填充它。
pro addskirt, v, c, newv, newc, exterior_indices = e compile_opt strictarr newv = !null newc = !null newv = v mins = min(newv, dimension = 2, max = maxs) for i = 0, 1 do begin newv[i, *] -= mins[i] endfor if (e eq !null) then begin triangulate, newv[0, *], newv[1, *], tr, b endif else begin b = e endelse xmid= total(newv[0, b])/b.length ymid = total(newv[1, b])/b.length xx = newv[0, b] - xmid yy = newv[1, b] - ymid theta = atan(yy, xx) s = sort(theta) b = b[s] plots, newv[0, [b, b[0]]], newv[1, [b, b[0]]] oldx = reform(v[0, *]) oldy = reform(v[1, *]) oldz = reform(v[2, *]) oldc = oldx.length newx = [oldx, oldx[[b, b[0]]]] newy = [oldy, oldy[[b, b[0]]]] newz = [oldz, replicate(mins[2] - (maxs[2] - mins[2])/100., b.length + 1)] b = [b, b[0]] newc = c for i = 0, b.length - 2 do begin newc = [newc, 3, oldc + i, oldc + i + 1, b[i]] newc = [newc, 3, oldc + i + 1, b[i + 1], b[i]] endfor t = idlgrtessellator() t.addpolygon, oldx[b], oldy[b] r = t.tessellate(basev, basep) for i = 0, basep.length/4 - 1 do begin newc = [newc, 3, oldc + basep[[i*4 + 1 + lindgen(3)]]] endfor newv = transpose([[newx], [newy], [newz]]) newc mod= newv.length/3 opoly = idlgrpolygon(newv, poly = newc, style = 1, color = !color.red) xobjview, opoly end
此例程接受原始顶点和连接性列表以及外部顶点索引作为输入，并创建一组新的顶点和连接性。
IDL> addskirt, v, c, newv, newc, exterior_indices = e
出于调试目的，此调用还在 IDL 中可视化结果。
现在我们有了大多数打印设备可以接受的顶点集和连接性列表。但我们必须创建一个文件，以打印机可以使用的格式包含这些坐标，或者通过一个接口将数据转换为特定打印机可接受的格式。
大多数打印机驱动程序都接受的一种非常简单的文件格式是 STL，它是“STereoLithography（立体光刻）”的缩写。其最简单的形式只是一个三角形列表。
以下代码使用我们更新后的顶点和连接性列表编写 ASCII 格式的 STL 文件。
Pro MakeSTL, vertices, polygons, outfile On_Error, 2 Catch, ErrorNumber If (ErrorNumber ne 0) then Begin Catch, /Cancel If (LUN ne !null) then Begin Free_LUN, LUN EndIf If (C ne !null) then Begin CD, C Endif Message, /Reissue_Last Endif L = List() L.Add, 'solid mntprinceton' V = vertices For I = 0, 2 do begin V[I, *] -= min(Vertices[I, *]) EndFor Poly = Polygons For I = 0L, Poly.Length/4 - 1 Do Begin V1 = Reform(V[*, Poly[I*4 + 1]]) V2 = Reform(V[*, Poly[I*4 + 2]]) V3 = Reform(V[*, Poly[I*4 + 3]]) uu = v2 - v1 vv = v3 - v1 nx = uu[1]*vv[2] - uu[2]*vv[1] ny = uu[2]*vv[0] - uu[0]*vv[2] nz = uu[0]*vv[1] - uu[1]*vv[0] R = Sqrt(Total([nx,ny,nz]^2)) L.Add, 'facet normal ' + ((String([nx,ny,nz]/R, Format = '(f8.5)')).Trim()).Join(' ') L.Add, ' outer loop' L.Add, ' vertex ' + $ ((String(V1)).Trim()).Join(' ') L.Add, ' vertex ' + $ ((String(V2)).Trim()).Join(' ') L.Add, ' vertex ' + $ ((String(V3)).Trim()).Join(' ') L.Add, ' endloop' L.Add, 'endfacet' EndFor L.Add, 'endsolid mntprinceton' OpenW, LUN, OutFile, /Get_LUN ForEach Line, L Do Begin PrintF, LUN, Line EndForEach Free_LUN, LUN End
使用新的顶点和连接性列表以及输出 .stl 文件的路径调用此例程。
IDL> makestl, newv, newc, "path-to-your-output-file.stl"
输出文件将是简单的 ASCII 格式，内容类似于以下内容。
打印 STL 文件的一个选项是通过“3D Builder”应用程序，该应用程序可作为 Windows 10 的免费软件获得。
选择 STL 文件作为输入。这会将模型导入打印机驱动程序，该驱动程序负责将原始数据居中并缩放到打印机空间。
现在，如果您有打印机，就可以打印您的数据了。3D Builder 应用程序还有一个方便（尽管相当昂贵）的选项，可以将您的模型导出到第三方在线服务，该服务将远程打印您的模型并邮寄给您。
这是最终产品。请注意，细小的徒步小径清晰可见。
如果您在 ENVI 或 IDL 中有一个更复杂的模型，并且希望获得一些帮助来生成超出本文所涵盖范围的实体渲染，请联系 NV5 Geospatial 的 Custom Solutions Group（自定义解决方案组）。我们很乐意提供帮助！