快速读取并展示新文件格式数据——使用IDL进行古生物学研究
17374 为本文评分:
3.3
快速读取并展示新文件格式数据——使用IDL进行古生物学研究
Jim Pendleton 2016年3月17日，星期四
随着时间的推移，你可能会遇到想在IDL中读取的数据集，但这些文件的格式在IDL语言中缺乏内置的实用工具来支持。
在为定制解决方案组工作时，我经常使用以下技术为IDL构建自定义数据读取器。
最近我读到的一篇文章描述了一个新创建的古生物学数据存储库，其中包含了已灭绝物种骨骼的3D模型。MorphoSource.org是杜克大学的一个项目，部分资金由美国国家科学基金会资助。根据其网站说明，对所提供媒体的商业使用是被严格禁止的。
我的目标是看看读取其中一些数据并在交互式IDL显示窗口中展示它有多容易。事实证明，这非常简单！
在数据库中搜索后，我发现了许多有趣的数据类型，包括3D模型。上图所示即为一个例子，版权归属于佛罗里达自然历史博物馆，这是Calippus elachistus（一种早期类马哺乳动物，可能在中新世时期出没于佛罗里达地区）的一颗牙齿。
我向MorphoSource.org申请了一个用户账户，他们的机器人给我回复了许可，允许我抓取一些数据，当然，仅限于教学和非商业用途。
进一步深入查看后，我找到了一些ASCII格式的"Avizo PLY"文件。我以前从未处理过这种类型的文件，而且IDL也没有内置的读取器，但用标准文本编辑器打开文件可以看到，其格式非常直接。
有了这些信息，就很容易在IDL中构建一个既能读取数据又能展示3D对象的实用工具。
虽然我们不会尝试解析文件中的所有数据，但可以快速定位可视化所需的两大块主要数据：顶点坐标列表以及这些顶点的连通性。
文件头以以下内容开始：
ply format ascii 1.0 obj_info Format: Avizo PLY element vertex 249741 property float32 x property float32 y property float32 z element face 495182 property list uint8 int32 vertex_indices
这告诉我，在文件的某个地方有249741个带有浮点(X,Y,Z)坐标的顶点，以及495182个"面"或连通性列表。一个连通性列表告诉我们顶点列表中的索引是如何相互连接的。
（在MorphoSource数据集的其它文件中，你可能还会找到顶点着色的数据，这增加了更多趣味性！）
文件头中还有一些额外的行，我们将忽略它们。文件头以"end_header"行结束。
紧接在文件头之后的行是一系列三个浮点元素，例如：
end_header 34.6223 25.3743 11.5595 34.5367 25.4449 11.5501 34.6183 25.4799 11.5998
天哪，这些看起来像是顶点坐标，不是吗？现在我们凭直觉可以推测，如果顺着文件往下看249741行（顶点数量），我们会发现不同格式的数据。果然，我们找到了一个过渡点。
23.9367 23.4851 43.5914 24.1271 23.4578 43.5783 24.0427 23.5348 43.5879 3 0 4 1 0 3 1 2 0 0 3 3 8 7 0 3 8 3 0 0
每一行的连通性列表以一个多边形面包含的顶点数量开始。在这个例子中，它们都是三角形，因此值是"3"。顶点数量之后是构成每个三角形面的三个从零开始的顶点索引。每行的最后一个数字可能表示多边形是开放还是闭合的。它们的值都是0，所以我们直接忽略第五个值。
并非巧合，这些顶点列表和连通性列表的格式正好是IDL的IDLgrPolygon对象所使用的格式。
将文件下载到我的桌面后，我将把文件内容读入一个名为"val"的字符串数组变量。
file = 'UF-53577_M5794-5222_Calippus_elachistus.ply'
n = file_lines(file)
val = strarr(n)
openr, lun, file, /get_lun
readf, lun, val
free_lun, lun
定位以文本"element vertex"开头的行。我们将在空格处分割该行并提取第三个元素。那就是顶点数量。
vertexline = where(val.startswith('element vertex'))
vertexcount = long((val[vertexline[0]].split(' '))[2])
对以"element face"开头的多边形连通性数量重复此操作。
polygonline = where(val.startswith('element face'))
polygoncount = long((val[polygonline[0]].split(' '))[2])
创建用于存储顶点和多边形的IDL数组。
vertices = fltarr(3, vertexcount)
poly = lonarr(4, polygoncount)
从字符串数组中移除文件头，这样我们就定位到第一行顶点数据了。
last = (where(val.startswith('end_header')))[0]
val = val[last + 1: *]
循环遍历包含顶点数据的行，并将值提取到顶点数组中。
for i = 0L, vertexcount - 1 do $
vertices[0, i] = float(((val[i]).trim()).split(' '))
从字符串数组的开头移除顶点数据，并对连通性数据重复此过程。
val = val[vertexcount:*]
for i = 0L, polygoncount - 1 do $
poly[0, i] = (long(((val[i]).trim()).split(' ')))[0:3]
对于显示，我们只需要两行代码。用我们的数据填充一个IDLgrPolygon对象。
opolygon = idlgrpolygon(vertices, poly = poly, style = 2, color = !color.gold)
最后，在IDL的一个交互式工具XOBJVIEW.中显示数据。
xobjview, opolygon
这匹马有蛀牙吗？还是牙冠上的洞是为了取样而钻的？
下面是一个示例例程，它整合了上面展示的代码，并在你的临时目录中生成一个MP4动画。你需要申请一个MorphoSource账户才能下载你自己的PLY格式文件作为输入。欢迎你泛化此代码以获得一个更健壮的PLY格式解决方案。
pro read_ply, dir, vertices
file = filepath('UF-53577_M5794-5222_Calippus_elachistus.ply')
n = file_lines(file)
val = strarr(n)
openr, lun, file, /get_lun
readf, lun, val
free_lun, lun
vertexline = where(val.startswith('element vertex'))
vertexcount = long((val[vertexline[0]].split(' '))[2])
polygonline = where(val.startswith('element face'))
polygoncount = long((val[polygonline[0]].split(' '))[2])
vertices = fltarr(3, vertexcount)
poly = lonarr(4, polygoncount)
last = (where(val.startswith('end_header')))[0]
val = val[last + 1: *]
for i = 0L, vertexcount - 1 do $
vertices[0, i] = float(((val[i]).trim()).split(' '))
val = val[vertexcount:*]
for i = 0L, polygoncount - 1 do $
poly[0, i] = (long(((val[i]).trim()).split(' ')))[0:3]
;poly = reform(poly, poly.length, /overwrite)
om = idlgrmodel()
opolygon = idlgrpolygon(vertices, style = 2, poly = poly, color = !color.gold)
om->add, opolygon
mx = mean(vertices[0, *])
my = mean(vertices[1, *])
mz = mean(vertices[2, *])
om->translate, -mx, -my, -mz
xobjview, om, tlb = tlb, xsize = 256, ysize = 256
frames = list()
temppng = filepath(/tmp, 'frame.png')
for i = 0, 360, 3 do begin ; rotate
om->rotate, [1, 1, 0], 3
xobjview, refresh = tlb
xobjview_write_image, temppng, 'png'
frames.add, read_png(temppng)
file_delete, temppng
endfor
for i = 1, 100, 3 do begin ; zoom in
om->scale, 1.03, 1.03, 1.03
xobjview, refresh = tlb
xobjview_write_image, temppng, 'png'
frames.add, read_png(temppng)
file_delete, temppng
endfor
for i = 1, 100, 3 do begin ; zoom out
om->scale, 1/1.03, 1/1.03, 1/1.03
xobjview, refresh = tlb
xobjview_write_image, temppng, 'png'
frames.add, read_png(temppng)
file_delete, temppng
endfor
tempvid = filepath(/tmp, 'horsetooth.mp4')
ovid = IDLffVideoWrite(tempvid)
fs = (frames[0]).dim
vidstream = ovid.AddVideoStream(fs[1], fs[2], 24)
foreach frame, frames do $
t = ovid.put(vidstream, frame)
obj_destroy, ovid
spawn, tempvid
end