上一节课内容回顾 4.1 可视化数据类型 4.2 科学数据管理 4.3 科学数据格式 4.4 图形图像数据格式 4.5 VTK中的数据表示

网格数据: 笛卡儿网格 规整网格 矩形网格 结构网格 非结构网格 曲线网格 块结构网格 混合结构网格 散乱点

科学数据格式: HDF格式 CDF和NetCDF

VTK中的数据对象: vtkDataObject vtkPiecewiseFunction vtkDataSet vtkImageData vtkStructuredPoints vtkImageCanvasSource2D vtkPointSet vtkPolyData vtkStructuredGrid vtkUnstructuredGrid vtkRectilinearGrid vtkDataSetAttributes vtkPointData vtkCellData vtkFieldData vtkAttributeData

vtkCell vtkHexahedron vtkLine vtkPixel vtkPolyLine vtkPolyVertex vtkPolygon vtkQuad vtkTetra vtkTriangle vtkTriangleStrip vtkVertex vtkVoxel vtkWedge vtkPyramid

本节课内容 5.1 网格序列法 5.2 网格无关法 5.3 区域填充法 5.4 曲面标量场可视化

标量场数据是环境科学与工程中应用最多的类型。 如污染物浓度分布数据、温度分布数据、高程数据等。 标量场可视化按数据划分为： 二维数据  等值线法，分层设色法 三维数据  等值面法 断层数据  表面重构 体素数据  体可视化

二维标量场可看成是位于某二维面上的一维标量函数F(x,y). 常见的如数字高程模型DEM、TIN、数字地形模型等。基本上是以规整格网或非规整格网表示某二维面，二维面网格上每一点有二维空间坐标（x，y），在网格点或网格拓扑单元上附属有标量数据。 二维标量场的可视化主要任务是把二维面上的标量数据的分布特征表现出来。主要有等值线法和分层设色法。

二维标量场的数据表示： （1） 结构点法 VTK DataFile Version 3.0 DEM ASCII DATASET STRUCTURED_POINTS DIMENSIONS 512 512 1 ORIGIN 0 0 0 SPACING 100 100 1 POINT_DATA 262144 SCALARS scalars double LOOKUP_TABLE default 22620.2 22620.3 22620.1 22620.2 22632 22632 22631.9 22631.8 ...

（2） 多边形法 vtk DataFile Version 3.0 TIN ASCII DATASET POLYDATA POINTS 13046 float 1428.22 2526.86 22620.2 1452.64 2526.86 22620.3 1428.22 2502.44 22620.1 ... POLYGONS 24112 96448 3 1129 1128 1127 3 10155 1664 10153 ... POINT_DATA 13046 SCALARS scalars double LOOKUP_TABLE default 22620.2 22620.3 22620.1 22620.2 22632 22632 22631.9 22631.8 ...

二维标量场的数据生成： 离散点插值 构造三角网 delny = vtk.vtkDelaunay2D() delny.SetInput(profile) delny.SetSource(profile) delny.SetTolerance(0.001) #delny.SetAlpha(0.1) #delny.BoundingTriangulationOn() #delny.SetOffset(10.00) 演示

三角网： mapMesh = vtk.vtkPolyDataMapper() # Generate some random points math = vtk.vtkMath() points = vtk.vtkPoints() for i in range(0, 100): points.InsertPoint(i, math.Random(0, 1), math.Random(0, 1), 0.0) polys = vtk.vtkCellArray() polys.InsertNextCell(3) polys.InsertCellPoint(4) polys.InsertCellPoint(5) polys.InsertCellPoint(7) # Create a polydata with the points we just created. profile = vtk.vtkPolyData() profile.SetPoints(points) profile.SetLines(polys) # Perform a 2D Delaunay triangulation on them. delny = vtk.vtkDelaunay2D() delny.SetInput(profile) #delny.SetSource(profile) delny.SetTolerance(0.001) #delny.SetAlpha(0.1) #delny.BoundingTriangulationOn() #delny.SetOffset(10.00) mapMesh = vtk.vtkPolyDataMapper() mapMesh.SetInput(delny.GetOutput()) meshActor = vtk.vtkActor() meshActor.SetMapper(mapMesh) meshActor.GetProperty().SetColor(.1, .2, .4) profileMesh = vtk.vtkPolyDataMapper() profileMesh.SetInput(profile) profileActor = vtk.vtkActor() profileActor.SetMapper(profileMesh) profileActor.GetProperty().SetColor(1, 0, 0) extract = vtk.vtkExtractEdges() extract.SetInput(delny.GetOutput()) tubes = vtk.vtkTubeFilter() tubes.SetInput(extract.GetOutput()) tubes.SetRadius(0.0025) tubes.SetNumberOfSides(6) mapEdges = vtk.vtkPolyDataMapper() mapEdges.SetInput(tubes.GetOutput())

edgeActor = vtk.vtkActor() edgeActor.SetMapper(mapEdges) edgeActor.GetProperty().SetColor(peacock) edgeActor.GetProperty().SetSpecularColor(1, 1, 1) edgeActor.GetProperty().SetSpecular(0.3) edgeActor.GetProperty().SetSpecularPower(20) edgeActor.GetProperty().SetAmbient(0.2) edgeActor.GetProperty().SetDiffuse(0.8) ball = vtk.vtkSphereSource() ball.SetRadius(0.005) ball.SetThetaResolution(12) ball.SetPhiResolution(12) balls = vtk.vtkGlyph3D() balls.SetInput(delny.GetOutput()) balls.SetSource(ball.GetOutput()) mapBalls = vtk.vtkPolyDataMapper() mapBalls.SetInput(balls.GetOutput()) ballActor = vtk.vtkActor() ballActor.SetMapper(mapBalls) ballActor.GetProperty().SetColor(hot_pink) ballActor.GetProperty().SetSpecularColor(1, 1, 1) ballActor.GetProperty().SetSpecular(0.3) ballActor.GetProperty().SetSpecularPower(20) ballActor.GetProperty().SetAmbient(0.2) ballActor.GetProperty().SetDiffuse(0.8) # Create the rendering window, renderer, and interactive renderer ren = vtk.vtkRenderer() renWin = vtk.vtkRenderWindow() renWin.AddRenderer(ren) iren = vtk.vtkRenderWindowInteractor() iren.SetRenderWindow(renWin) # Add the actors to the renderer, set the background and size ren.AddActor(profileActor) ren.AddActor(ballActor) ren.AddActor(edgeActor) ren.SetBackground(1, 1, 1) renWin.SetSize(150, 150) ren.GetActiveCamera().Zoom(1.5) # Interact with the data. iren.Initialize() renWin.Render() iren.Start()

离散点生成TIN或DEM del =vtkDelaunay2D() del.SetInput( polyData) vtkPoints points points InsertPoint 0 1 4 0 points InsertPoint 1 3 4 0 points InsertPoint 2 7 4 0 ..... points InsertPoint 32 7.5 9 0 points InsertPoint 33 7.5 11 0 points InsertPoint 34 6.5 11 0 points InsertPoint 35 6.5 9 0 points InsertPoint 36 5 9 0 points InsertPoint 37 4 6 0 points InsertPoint 38 3 9 0 points InsertPoint 39 2 9 0 vtkCellArray polys polys InsertNextCell 12 polys InsertCellPoint 0 polys InsertCellPoint 11 polys InsertNextCell 28 polys InsertCellPoint 39 .... polys InsertCellPoint 12 vtkPolyData polyData polyData SetPoints points polyData SetPolys polys del =vtkDelaunay2D() del.SetInput( polyData) mapMesh =vtkPolyDataMapper() mapMesh.SetInput(del.GetOutput()) meshActor=vtkActor() meshActor.SetMapper( mapMesh ) 例子

等值线生成TIN或DEM 例子 del =vtkDelaunay2D() del.SetInput( polyData) del.SetSource( polyData) mapMesh =vtkPolyDataMapper() mapMesh.SetInput(del.GetOutput()) meshActor=vtkActor() meshActor.SetMapper( mapMesh ) 例子

vtkDelaunay2D等说明： virtual void SetInput (vtkPointSet *input) 设置离散点集数据 void SetSource (vtkPolyData *) 设置多边形数据，对生成的三角网进行约束 virtual void SetAlpha (double) 设置输出结果的范围（MESH结点为中心的圆的半径） virtual void SetTransform (vtkAbstractTransform *) 设置变换矩阵 要把生成的三角网转成DEM可以使用 vtkShepardMethod 等类

二维标量场可视化方法： 等值线法 等值线是由F(xi,yi) == Ft的所有点(xi,yi)组成。这些点 按一定顺序连接组成函数F(x,y)的值为Ft时的等值线。 等值线的提取分成： 网格序列法 网格无关法 分层设色法 线性插值 非线性插值

5.1 网格序列法 正规化网格等值线生成 单元剖分法 连续光滑等值线生成 + - - + + - + - + - - + + - +

5.2 网格无关法 通过给定等值线的起始点（或求出），利用该点附近的 局部几何性质，计算该等值线的下一点，然后利用新点， 重复计算下一点直至到达区域边界，或回到始点。 步进法－－跟踪法 适应法－－逐步精化方法 P0 P1 P0 P1 S E

5.3 区域填充法 线性插值下的分层设色法 非线性插值下的分层设色法

5.4 曲面标量场可视化 曲面分布数据的插值方法 球面上的等值线

5.5 VTK中的二维标量场可视化 等值线可视化 ×vtkContourValues 管理、设置等值线的值 ×vtkLookupTable 在标量值和颜色之间进行匹配 ×vtkContourFilter 从标量数据中生成等值线 contours = vtkContourFilter() contours.SetInput(dataModel.GetOutput()) contours.GenerateValues(10,zmin,zmax) contMapper = vtkPolyDataMapper() contMapper.SetInput(contours.GetOutput())

地形数据等值线可视化实例 import os try: VTK_DATA = os.environ['VTKDATA'] except KeyError: VTK_DATA = '.' from vtk import * ren = vtkRenderer() renWin = vtkRenderWindow() renWin.AddRenderer(ren) iren = vtkRenderWindowInteractor() iren.SetRenderWindow(renWin)

demModel = vtkDEMReader() demModel.SetFileName( VTK_DATA + "/albany-w" ) demModel.Update() print demModel ren.SetBackground( 0.4, 0.4,0.4) ze = demModel.GetElevationBounds() contours = vtkContourFilter() contours.SetInput(demModel.GetOutput()) contours.GenerateValues(10,ze[0],ze[1]) contMapper = vtkPolyDataMapper() contMapper.SetInput(contours.GetOutput()) contMapper.SetScalarRange(ze[0],ze[1]) contMapper.Update()

contActor = vtkActor() contActor.SetMapper(contMapper) ren.AddActor( contActor ) iren.Initialize() ren.GetActiveCamera().SetViewUp( 0, 0, 1) ren.GetActiveCamera().SetPosition( -99900, -21354, 131801) ren.GetActiveCamera().SetFocalPoint( 41461, 41461, 2815) ren.GetActiveCamera().ComputeViewPlaneNormal() ren.ResetCamera() ren.GetActiveCamera().Dolly( 1.2 ) ren.ResetCameraClippingRange() renWin.Render() iren.Start() 运行

调整颜色： lut = vtkLookupTable() lut.SetHueRange( 0.6, 0 ) lut.SetSaturationRange( 1.0, 0) lut.SetTableRange( ze[0],ze[1]) contMapper = vtkPolyDataMapper() contMapper.SetInput(contours.GetOutput()) contMapper.SetScalarRange(ze[0],ze[1]) contMapper.SetLookupTable(lut) contMapper.Update() 运行 选择待显示的等高线： contours.GenerateValues(10,ze[0],ze[1]) contours. SetValue (i, ContourValue)

分成色设法： 运行 曲面标量场可视化： vtkElevationFilter SetLowPoint( 0, 0, lo) SetHighPoint( 0, 0, hi) SetScalarRange( lo, hi) vtkBandedPolyDataContourFilter GenerateValues ( Num, lo, hi) SetScalarModeToIndex /Value GenerateContourEdgesOn /Off ClippingOn

分成色设法1: toPolylod = vtkCastToConcrete() demModel = vtkDEMReader() demModel.SetFileName( VTK_DATA + "/albany-w" ) shrinklod = vtkImageShrink3D() shrinklod.SetShrinkFactors( lod,lod, 1) shrinklod.SetInput(demModel.GetOutput()) shrinklod.AveragingOn() geomlod = vtkStructuredPointsGeometryFilter() geomlod.SetInput(shrinklod.GetOutput()) geomlod.ReleaseDataFlagOn() warplod = vtkWarpScalar() warplod.SetInput(geomlod.GetOutput()) warplod.SetNormal( 0, 0, 1) warplod.UseNormalOn() warplod.SetScaleFactor(Scale) warplod.ReleaseDataFlagOn() elevationlod = vtkElevationFilter() elevationlod.SetInput(warplod.GetOutput()) elevationlod.SetLowPoint( 0, 0, lo) elevationlod.SetHighPoint( 0, 0, hi) elevationlod.SetScalarRange( lo, hi) elevationlod.ReleaseDataFlagOn() toPolylod = vtkCastToConcrete() toPolylod.SetInput(elevationlod.GetOutput()) normalslod = vtkPolyDataNormals() normalslod.SetInput(toPolylod.GetPolyDataOutput()) normalslod.SetMaxRecursionDepth(1000) normalslod.SetFeatureAngle(60) normalslod.ConsistencyOff() normalslod.SplittingOff() normalslod.ReleaseDataFlagOn() demMapperlod = vtkPolyDataMapper() demMapperlod.SetInput(normalslod.GetOutput()) demMapperlod.SetScalarRange(lo, hi) demMapperlod.SetLookupTable(lut) demMapperlod.ImmediateModeRenderingOn()

分成色设法2: demModel = vtk.vtkDEMReader() demModel.SetFileName("./SainteHelens.dem") demModel.Update() shrinkFactor = 4 shrink = vtk.vtkImageShrink3D() shrink.SetShrinkFactors(shrinkFactor, shrinkFactor, 1) shrink.SetInput(demModel.GetOutput()) shrink.AveragingOn() # Convert the image into polygons. geom = vtk.vtkImageDataGeometryFilter() geom.SetInput(shrink.GetOutput()) # Warp the polygons based on elevation. warp = vtk.vtkWarpScalar() warp.SetInput(geom.GetOutput()) warp.SetNormal(0, 0, 1) warp.UseNormalOn() warp.SetScaleFactor(Scale) # Create the contour bands. bcf = vtk.vtkBandedPolyDataContourFilter() bcf.SetInput(warp.GetPolyDataOutput()) bcf.GenerateValues(15, demModel.GetOutput().GetScalarRange()) bcf.ClippingOn() bcf.SetScalarModeToIndex() bcf.GenerateContourEdgesOn()

# Compute normals to give a better look. normals = vtk.vtkPolyDataNormals() normals.SetInput(bcf.GetOutput()) normals.SetFeatureAngle(60) normals.ConsistencyOff() normals.SplittingOff() demMapper = vtk.vtkPolyDataMapper() demMapper.SetInput(normals.GetOutput()) demMapper.SetScalarRange(0, 10) demMapper.SetLookupTable(lutLand) demMapper.SetScalarModeToUseCellData() demActor = vtk.vtkLODActor() demActor.SetMapper(demMapper) ## Create contour edges edgeMapper = vtk.vtkPolyDataMapper() edgeMapper.SetInput(bcf.GetContourEdgesOutput()) edgeMapper.SetResolveCoincidentTopologyToPolygonOffset() edgeActor = vtk.vtkActor() edgeActor.SetMapper(edgeMapper) edgeActor.GetProperty().SetColor(0, 0, 0)

## Test clipping # Create the contour bands. bcf2 = vtk.vtkBandedPolyDataContourFilter() bcf2.SetInput(warp.GetPolyDataOutput()) bcf2.ClippingOn() bcf2.GenerateValues(10, 0, 1100) bcf2.SetScalarModeToValue() # Compute normals to give a better look. normals2 = vtk.vtkPolyDataNormals() normals2.SetInput(bcf2.GetOutput()) normals2.SetFeatureAngle(60) normals2.ConsistencyOff() normals2.SplittingOff() lut = vtk.vtkLookupTable() lut.SetNumberOfColors(10) demMapper2 = vtk.vtkPolyDataMapper() demMapper2.SetInput(normals2.GetOutput()) demMapper2.SetScalarRange(demModel.GetOutput().GetScalarRange()) demMapper2.SetLookupTable(lut) demMapper2.SetScalarModeToUseCellData() demActor2 = vtk.vtkLODActor() demActor2.SetMapper(demMapper2) demActor2.AddPosition(0, 15000, 0)

# Create the RenderWindow, Renderer and both Actors ren = vtk.vtkRenderer() renWin = vtk.vtkRenderWindow() renWin.AddRenderer(ren) iren = vtk.vtkRenderWindowInteractor() iren.SetRenderWindow(renWin) # Add the actors to the renderer, set the background and size ren.AddActor(demActor) ren.AddActor(demActor2) ren.AddActor(edgeActor) ren.SetBackground(.4, .4, .4) renWin.SetSize(375, 200) cam = vtk.vtkCamera() cam.SetPosition(-17438.8, 2410.62, 25470.8) cam.SetFocalPoint(3985.35, 11930.6, 5922.14) cam.SetViewUp(0, 0, 1) ren.SetActiveCamera(cam) ren.ResetCameraClippingRange() iren.Initialize() renWin.Render() iren.Start()

5.6 纹理可视化 可以把2维标量场数据（如污染物浓度分布）转成图像，以纹理的 形式贴在某个区域上。

纹理可视化基本用法 #读入图片 pnmReader = vtkPNMReader() pnmReader.SetFileName(VTK_DATA+"/earth.ppm") #纹理对象 aText = vtkTexture() aText.SetInput(pnmReader.GetOutput()) aText.InterpolateOn() #将纹理贴到演员对象上 aActor.SetTexture(atext) 演示

下一节课内容： 第六讲 断层间表面重构 拓扑重构 三角片表面几何重构 连续表面重构 体的重构