VTK Learning Twenty-three- PointInterpolator2D（一）

最新推荐文章于 2023-09-26 10:39:08 发布

wblong_cs

最新推荐文章于 2023-09-26 10:39:08 发布

阅读量483

点赞数

分类专栏： VTK 文章标签： VTK 插值高斯核函数 Interpolator

本文链接：https://blog.csdn.net/mrbaolong/article/details/105883499

版权

VTK 专栏收录该内容

72 篇文章 39 订阅

订阅专栏

VTK Learning Twenty-three- PointInterpolator2D

Description

use terrain dataset. Just for kicks we’ll treat the elevations as a “point cloud” and interpolate the elevation onto a plane.
使用DEM数据。在平面上内插高程值(忽略z值)。vtkPointInterpolator2D在XY平面上插值属性值或z值。插值结果默认保存到以Elevation为名称的属性数组中。

Code

python

#!/usr/bin/env python
import vtk
from vtk.util.misc import vtkGetDataRoot
VTK_DATA_ROOT = vtkGetDataRoot()
VTK_DATA_ROOT = "D:/codelearnning/vtk/VTK_BUILD/ExternalData/Testing/"
# Parameters for debugging
res = 200
math = vtk.vtkMath()

# create pipeline: use terrain dataset. Just for kicks we'll treat the elevations
# as a "point cloud" and interpolate the elevation onto a plane.
#
# Read the data: a height field results
demReader = vtk.vtkDEMReader()
demReader.SetFileName(VTK_DATA_ROOT + "/Data/SainteHelens.dem")
demReader.Update()
# 获取高程值的最大最小值
lo = demReader.GetOutput().GetScalarRange()[0]
hi = demReader.GetOutput().GetScalarRange()[1]

bds = demReader.GetOutput().GetBounds()
center = demReader.GetOutput().GetCenter()

# Create a plane of onto which to map the elevation data
plane = vtk.vtkPlaneSource()
plane.SetResolution(res,res)
plane.SetOrigin(bds[0],bds[2],bds[4])# minx,miny,minz
print(bds[0],bds[2],bds[4])
plane.SetPoint1(bds[1],bds[2],bds[4])# maxx,miny,minz
print(bds[1],bds[2],bds[4])
plane.SetPoint2(bds[0],bds[3],bds[4])# mixx,maxy,minz
print(bds[0],bds[3],bds[4])
plane.Update()


# Gaussian kernel-------------------------------------------------------
gaussianKernel = vtk.vtkGaussianKernel()
gaussianKernel.SetSharpness(2)
gaussianKernel.SetRadius(200)

interp = vtk.vtkPointInterpolator2D()
interp.SetInputConnection(plane.GetOutputPort())
interp.SetSourceConnection(demReader.GetOutputPort())
interp.SetKernel(gaussianKernel)
interp.SetNullPointsStrategyToClosestPoint()
interp.GetLocator().SetNumberOfPointsPerBucket(1)
interp.InterpolateZOff()

# Time execution
timer = vtk.vtkTimerLog()
timer.StartTimer()
interp.Update()
timer.StopTimer()
time = timer.GetElapsedTime()
print("Interpolate Terrain Points (Gaussian): {0}".format(time))

scalarRange = interp.GetOutput().GetPointData().GetArray("Elevation").GetRange()

lut = vtk.vtkLookupTable()
lut.SetHueRange(0.6, 0)
lut.SetSaturationRange(1.0, 0)
lut.SetValueRange(0.5, 1.0)

intMapper = vtk.vtkPolyDataMapper()
intMapper.SetInputConnection(interp.GetOutputPort())
intMapper.SetScalarModeToUsePointFieldData()
intMapper.SelectColorArray("Elevation")
intMapper.SetScalarRange(scalarRange)
intMapper.SetLookupTable(lut)

intActor = vtk.vtkActor()
intActor.SetMapper(intMapper)

# Create some contours
cf = vtk.vtkContourFilter()
cf.SetInputConnection(interp.GetOutputPort())
cf.GenerateValues(20,scalarRange)

cfMapper = vtk.vtkPolyDataMapper()
cfMapper.SetInputConnection(cf.GetOutputPort())

cfActor = vtk.vtkActor()
cfActor.SetMapper(cfMapper)

# Create the RenderWindow, Renderer and both Actors
#
ren0 = vtk.vtkRenderer()
renWin = vtk.vtkRenderWindow()
renWin.AddRenderer(ren0)
iren = vtk.vtkRenderWindowInteractor()
iren.SetRenderWindow(renWin)

# Add the actors to the renderer, set the background and size
#
ren0.AddActor(intActor)
ren0.AddActor(cfActor)
ren0.SetBackground(0.1, 0.2, 0.4)

renWin.SetSize(300, 300)

cam = ren0.GetActiveCamera()
cam.SetFocalPoint(center)

fp = cam.GetFocalPoint()
cam.SetPosition(fp[0]+.2,fp[1]+.1,fp[2]+1)
ren0.ResetCamera()

iren.Initialize()

# render the image
#
renWin.Render()
iren.Start()

cplusplus

#include <vtkSmartPointer.h>

#include <vtkPolyDataMapper.h>
#include <vtkProperty.h>
#include <vtkActor.h>
#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>

#include <vtkDEMReader.h>
#include <vtkImageData.h>
#include <vtkPlaneSource.h>
#include <vtkGaussianKernel.h>
#include <vtkPointInterpolator2D.h>
#include <vtkPointData.h>
#include <vtkLookupTable.h>
#include <vtkStaticPointLocator.h>
#include <vtkTimerLog.h>
#include <vtkContourFilter.h>

int main()
{
	
	//加载dem
	vtkNew<vtkDEMReader>reader;
	reader->SetFileName("D:/codelearnning/vtk/VTK_BUILD/ExternalData/" 
		"Testing/Data/SainteHelens.dem");
	reader->Update();
	//高程的范围
	double lo=reader->GetOutput()->GetScalarRange()[0];
	double hi = reader->GetOutput()->GetScalarRange()[1];
	
	double *bds=reader->GetOutput()->GetBounds();
	double *center = reader->GetOutput()->GetCenter();
	double*space=reader->GetOutput()->GetSpacing();
	std::cout << "Dem space :" << space[0] << ";" << space[1] << ";" << space[2] << std::endl;
	vtkNew<vtkPlaneSource>plane;
	plane->SetResolution(200, 200);
	plane->SetOrigin(bds[0], bds[2], bds[4]);
	plane->SetPoint1(bds[1], bds[2], bds[4]); 
	plane->SetPoint2(bds[0], bds[3], bds[4]);
	plane->Update();

    //高斯函数内核
	vtkNew<vtkGaussianKernel>gaussianKernel;
	//峰的窄宽，值越大对距离的敏感会减少
	gaussianKernel->SetSharpness(2);
	//搜索半径
	gaussianKernel->SetRadius(200);

	vtkNew<vtkPointInterpolator2D> interp;
	//待插值平面
	interp->SetInputConnection(plane->GetOutputPort());
	interp->SetSourceConnection(reader->GetOutputPort());
	interp->SetKernel(gaussianKernel);
	//没有零近点的时候取最近的点
	interp->SetNullPointsStrategyToClosestPoint();
	//建立插值索引定位器
	vtkStaticPointLocator *locator= vtkStaticPointLocator::SafeDownCast( interp->GetLocator());
	locator->SetNumberOfPointsPerBucket(1);
	//忽略z插值，属性插值
	interp->InterpolateZOff();

	vtkNew<vtkTimerLog>timer;
	timer->StartTimer();
	interp->Update();
	timer->StopTimer();
	double elapsedTime=timer->GetElapsedTime();
	std::cout << "Interpolate Terrain Points:" << elapsedTime << " seconds." << std::endl;
	//HSV
	double*scalarRange = interp->GetOutput()->GetPointData()->GetArray("Elevation")->GetRange();
	vtkNew<vtkLookupTable>lut;
	//色度，饱和度，亮度
	lut->SetHueRange(0.6, 0);//蓝色->红色
	lut->SetSaturationRange(.5,1);// 纯度（半蓝色）->近于白色
	lut->SetValueRange(0.5, 1.0);//亮度

	vtkNew<vtkPolyDataMapper>intMapper;
	intMapper->SetInputConnection(interp->GetOutputPort());
	intMapper->SetScalarModeToUsePointFieldData();
	intMapper->SelectColorArray("Elevation");
	intMapper->SetScalarRange(scalarRange);
	intMapper->SetLookupTable(lut);
	vtkNew<vtkActor>intActor;
	intActor->SetMapper(intMapper);

	//等高线
	vtkNew<vtkContourFilter>cf;
	cf->SetInputConnection(interp->GetOutputPort());
	cf->GenerateValues(20,scalarRange);
	vtkNew<vtkPolyDataMapper>cfMapper;
	cfMapper->SetInputConnection(cf->GetOutputPort());
    
	vtkNew<vtkActor>cfActor;
	cfActor->SetMapper(cfMapper);

	vtkSmartPointer<vtkRenderer> renderer  =
		vtkSmartPointer<vtkRenderer>::New(); 
	renderer->AddActor(intActor);
	renderer->AddActor(cfActor);
	vtkSmartPointer<vtkRenderWindow> renderWindow =
		vtkSmartPointer<vtkRenderWindow>::New();
	renderWindow->AddRenderer(renderer);
	renderWindow->SetSize(800, 600);

	vtkSmartPointer<vtkRenderWindowInteractor> interactor =
		vtkSmartPointer<vtkRenderWindowInteractor>::New();
	interactor->SetRenderWindow(renderWindow);
	renderer->SetBackground(0, 0, 0);
	renderWindow->Render();
	interactor->Start();
	return 0;

}

Result

在这里插入图片描述

Test

把Z作为属性来插值。vtkPoints作为插值的数据源，插值的结果数据以属性数据的形式存在于输出数据中。

//原始采样点数据
	vtkNew<vtkPoints>src_pts;
	src_pts->InsertNextPoint(0, 0, 99);
	src_pts->InsertNextPoint(100, 100, 100);
	src_pts->InsertNextPoint(0, 100, 98);
	src_pts->InsertNextPoint(100, 0, 97);
//待插值的点数据，z值不计
	vtkNew<vtkPoints>dest_pts;
	for (int i = 0; i < 100; i++) {
		dest_pts->InsertNextPoint(vtkMath::Random(-50, 150), vtkMath::Random(-50, 150), 0);
	}
	vtkNew<vtkPolyData>srcdata;
	srcdata->SetPoints(src_pts);

	vtkNew<vtkPolyData>destdata;
	destdata->SetPoints(dest_pts);

	//高斯函数内核
	vtkNew<vtkGaussianKernel>gaussianKernel;
	gaussianKernel->SetSharpness(2);
	gaussianKernel->SetRadius(100);

	vtkNew<vtkPointInterpolator2D> interp;
	//
	interp->SetInputData(destdata);
	interp->SetSourceData(srcdata);
	interp->SetKernel(gaussianKernel);
	//没有零近点的时候取最近的点
	interp->SetNullPointsStrategyToClosestPoint();
	//建立插值索引定位器
	vtkStaticPointLocator *locator = vtkStaticPointLocator::SafeDownCast(interp->GetLocator());
	locator->SetNumberOfPointsPerBucket(1);
	interp->InterpolateZOn();
	/**
	* Specify the name of the output array containing z values. This method is
	* only applicable when InterpolateZ is enabled. By default the output
	* array name is "Elevation".
	*/
	interp->SetZArrayName("Z");

	vtkNew<vtkTimerLog>timer;
	timer->StartTimer();
	interp->Update();
	timer->StopTimer();
	double elapsedTime = timer->GetElapsedTime();
	std::cout << "Interpolate Terrain Points:" << elapsedTime << " seconds." << std::endl;
	
	vtkDataArray* res = interp->GetOutput()->GetPointData()->GetArray("Z");
	double p[3];
    //插值结果
	for (int i = 0; i < res->GetNumberOfValues(); ++i) {
		dest_pts->GetPoint(i, p);
		//Z值
		p[2] = res->GetVariantValue(i).ToDouble();
		std::cout << i << ":" << p[0] << ";" << p[1] << ";" <<p[2]<< std::endl;
	}