本文的主要内容：简单介绍VTK中有向包围盒的相关功能。 主要涉及vtkOBBTree类。
哪些人适合阅读本文：有一定VTK基础的人。
一、vtkOBBTree
vtkOBBTree是用于产生有向包围盒的类。
VTK官网描述：
vtkOBBTree是一个用于生成定向边界框（OBB）树的对象。定向边界框是指不一定沿坐标轴对齐的边界框。OBB树是这种盒子的层次树结构，其中更深层次的OBB限制了更小的空间区域。
为了构建OBB，使用了一个递归的、自上而下的过程。首先，通过找到定义数据集的单元格（及其点）的均值和协方差矩阵来构建根OBB。提取协方差矩阵的特征向量，给出一组三个正交向量，定义最紧拟合的OBB。要创建两个子OBB，需要找到一个（近似）将数字单元格一分为二的分割平面。然后将这些分配给子OBB。然后，此过程继续进行，直到MaxLevel ivar限制递归，或者找不到分割平面。
1.1 几种树结构的对比
对比的几种树有：vtkKdTreePointLocator，vtkOBBTree，vtkOctreePointLocator，vtkModifiedBSPTree。
核心代码如下：
#include <vtkActor.h>
#include <vtkCamera.h>
#include <vtkInteractorStyleTrackballCamera.h>
#include <vtkNew.h>
#include <vtkPointSource.h>
#include <vtkPolyData.h>
#include <vtkPolyDataMapper.h>
#include <vtkProperty.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkRenderer.h>
#include <vtkSmartPointer.h>
#include <vtkXMLPolyDataReader.h>
#include <vtkKdTreePointLocator.h>
#include <vtkModifiedBSPTree.h>
#include <vtkOBBTree.h>
#include <vtkOctreePointLocator.h>
#include <vtkPLYReader.h>
#include <vector>
namespace {
class KeyPressInteractorStyle : public vtkInteractorStyleTrackballCamera
{
public:
static KeyPressInteractorStyle* New();
vtkSmartPointer<vtkPolyData> data;
std::vector<vtkRenderer*> renderers;
std::vector<vtkSmartPointer<vtkLocator>> trees;
std::vector<vtkSmartPointer<vtkPolyDataMapper>> mappers;
std::vector<vtkSmartPointer<vtkActor>> actors;
vtkSmartPointer<vtkPolyDataMapper> meshMapper;
vtkSmartPointer<vtkActor> meshActor;
void Initialize()
{
this->meshMapper->SetInputData(this->data);
for (unsigned int i = 0; i < 4; i++)
{
vtkSmartPointer<vtkPolyDataMapper> mapper =
vtkSmartPointer<vtkPolyDataMapper>::New();
this->mappers.push_back(mapper);
vtkSmartPointer<vtkActor> actor = vtkSmartPointer<vtkActor>::New();
actor->SetMapper(mapper);
actor->GetProperty()->SetRepresentationToWireframe();
this->actors.push_back(actor);
this->renderers[i]->AddActor(actor);
this->renderers[i]->AddActor(meshActor);
}
this->Level = 1;
std::cout << "Level = " << this->Level << std::endl;
this->ReDraw();
}
KeyPressInteractorStyle()
{
this->Level = 1;
vtkSmartPointer<vtkLocator> tree0 =
vtkSmartPointer<vtkKdTreePointLocator>::New();
this->trees.push_back(tree0);
vtkSmartPointer<vtkLocator> tree1 = vtkSmartPointer<vtkOBBTree>::New();
this->trees.push_back(tree1);
vtkSmartPointer<vtkLocator> tree2 =
vtkSmartPointer<vtkOctreePointLocator>::New();
this->trees.push_back(tree2);
vtkSmartPointer<vtkLocator> tree3 =
vtkSmartPointer<vtkModifiedBSPTree>::New();
this->trees.push_back(tree3);
this->meshMapper = vtkSmartPointer<vtkPolyDataMapper>::New();
this->meshActor = vtkSmartPointer<vtkActor>::New();
this->meshActor->SetMapper(this->meshMapper);
}
virtual void OnChar()
{
char ch = this->Interactor->GetKeyCode();
switch (ch)
{
case 'n':
this->Level++;
break;
case 'p':
if (this->Level > 1)
{
this->Level--;
}
break;
default:
std::cout << "An unhandled key was pressed." << std::endl;
break;
}
this->ReDraw();
// Forward events.
if (ch != 'p') // Don't forward the "pick" command.
{
vtkInteractorStyleTrackballCamera::OnChar();
}
}
void ReDraw()
{
std::cout << "Level " << this->Level << std::endl;
for (unsigned i = 0; i < 4; i++)
{
vtkSmartPointer<vtkLocator> tree = this->trees[i];
tree->SetDataSet(data);
tree->BuildLocator();
vtkSmartPointer<vtkPolyData> polydata =
vtkSmartPointer<vtkPolyData>::New();
std::cout << "Tree " << i << " has " << tree->GetLevel() << " levels."
<< std::endl;
if (this->Level > tree->GetLevel())
{
tree->GenerateRepresentation(tree->GetLevel(), polydata);
}
else
{
tree->GenerateRepresentation(this->Level, polydata);
}
this->mappers[i]->SetInputData(polydata);
}
this->Interactor->GetRenderWindow()->Render();
}
private:
int Level;
};
vtkStandardNewMacro(KeyPressInteractorStyle);
} // namespace
DataStructureComparison::DataStructureComparison(vtkRenderWindow* renderWin)
{
vtkNew<vtkPolyData> originalMesh;
if (true) // If a file name is specified, open and use the file.
{
vtkNew<vtkPLYReader> reader;
reader->SetFileName("../resources/Armadillo.ply");
reader->Update();
originalMesh->ShallowCopy(reader->GetOutput());
}
else // If a file name is not specified, create a random cloud of points.
{
vtkNew<vtkPointSource> sphereSource;
sphereSource->SetNumberOfPoints(1000);
sphereSource->Update();
originalMesh->ShallowCopy(sphereSource->GetOutput());
}
double numberOfViewports = 4.;
vtkNew<vtkRenderWindow> renderWindow;
renderWindow->SetSize(200 * numberOfViewports, 200); //(width, height)
renderWindow->SetWindowName("DataStructureComparison");
vtkNew<KeyPressInteractorStyle> style;
style->data = originalMesh;
vtkNew<vtkCamera> camera;
for (unsigned int i = 0; i < 4; i++)
{
vtkNew<vtkRenderer> renderer;
renderWindow->AddRenderer(renderer);
style->renderers.push_back(renderer);
renderer->SetViewport(static_cast<double>(i) / numberOfViewports, 0, static_cast<double>(i + 1) / numberOfViewports, 1);
renderer->SetBackground(0.2, 0.3, 0.4);
renderer->SetActiveCamera(camera);
}
vtkNew<vtkRenderWindowInteractor> renderWindowInteractor;
renderWindowInteractor->SetRenderWindow(renderWindow);
renderWindowInteractor->SetInteractorStyle(style);
style->Initialize();
style->renderers[0]->ResetCamera();
renderWindowInteractor->Start();
}
从左到右依次是：vtkKdTreePointLocator，vtkOBBTree，vtkOctreePointLocator，vtkModifiedBSPTree。
1.2 获取线段与数据集的交点
可以通过vtkOBBTree获取一条线段与数据集的所有交点和Cells。
核心代码如下：
#include <vtkExtractCells.h>
#include <vtkIdList.h>
#include <vtkLineSource.h>
#include <vtkNamedColors.h>
#include <vtkNew.h>
#include <vtkOBBTree.h>
#include <vtkPoints.h>
#include <vtkSphereSource.h>
#include <vtkActor.h>
#include <vtkDataSetMapper.h>
#include <vtkPolyDataMapper.h>
#include <vtkProperty.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkRenderer.h>
OBBTreeExtractCells::OBBTreeExtractCells(vtkRenderWindow* renderWin)
{
vtkNew<vtkNamedColors> colors;
vtkNew<vtkSphereSource> sphereSource;
sphereSource->SetPhiResolution(20);
sphereSource->SetThetaResolution(20);
sphereSource->Update();
// Create the locator
vtkNew<vtkOBBTree> tree;
tree->SetDataSet(sphereSource->GetOutput());
tree->BuildLocator();
// Intersect the locator with the line
double lineP0[3] = { -0.6, -0.6, -0.6 };
double lineP1[3] = { .6, .6, .6 };
vtkNew<vtkPoints> intersectPoints;
vtkNew<vtkIdList> intersectCells;
double tol = 1.e-8;
tree->SetTolerance(tol);
tree->IntersectWithLine(lineP0, lineP1, intersectPoints, intersectCells);
std::cout << "NumPoints: " << intersectPoints->GetNumberOfPoints()
<< std::endl;
// Display list of intersections.
double intersection[3];
for (int i = 0; i < intersectPoints->GetNumberOfPoints(); i++)
{
intersectPoints->GetPoint(i, intersection);
std::cout << "\tPoint Intersection " << i << ": " << intersection[0] << ", "
<< intersection[1] << ", " << intersection[2] << std::endl;
}
std::cout << "NumCells: " << intersectCells->GetNumberOfIds() << std::endl;
vtkIdType cellId;
for (int i = 0; i < intersectCells->GetNumberOfIds(); i++)
{
cellId = intersectCells->GetId(i);
std::cout << "\tCellId " << i << ": " << cellId << std::endl;
}
// Render the line, sphere and intersected cells.
vtkNew<vtkLineSource> lineSource;
lineSource->SetPoint1(lineP0);
lineSource->SetPoint2(lineP1);
vtkNew<vtkPolyDataMapper> lineMapper;
lineMapper->SetInputConnection(lineSource->GetOutputPort());
vtkNew<vtkActor> lineActor;
lineActor->SetMapper(lineMapper);
vtkNew<vtkPolyDataMapper> sphereMapper;
sphereMapper->SetInputConnection(sphereSource->GetOutputPort());
vtkNew<vtkActor> sphereActor;
sphereActor->SetMapper(sphereMapper);
sphereActor->GetProperty()->SetRepresentationToWireframe();
sphereActor->GetProperty()->SetColor(colors->GetColor3d("Gold").GetData());
vtkNew<vtkExtractCells> cellSource;
cellSource->SetInputConnection(sphereSource->GetOutputPort());
cellSource->SetCellList(intersectCells);
vtkNew<vtkDataSetMapper> cellMapper;
cellMapper->SetInputConnection(cellSource->GetOutputPort());
vtkNew<vtkActor> cellActor;
cellActor->SetMapper(cellMapper);
cellActor->GetProperty()->SetColor(colors->GetColor3d("Tomato").GetData());
vtkNew<vtkRenderer> renderer;
renderWin->AddRenderer(renderer);
renderer->AddActor(lineActor);
renderer->AddActor(sphereActor);
renderer->AddActor(cellActor);
renderer->SetBackground(colors->GetColor3d("CornflowerBlue").GetData());
renderWin->Render();
}
线段与数据球体有两个交点和相交的两个Cell，分别是：
NumPoints: 2
Point Intersection 0: -0.285696, -0.285696, -0.285696
Point Intersection 1: 0.285696, 0.285696, 0.285696
NumCells: 2
CellId 0: 473
CellId 1: 116
1.3 OBB树可视化
#include <vtkActor.h>
#include <vtkCommand.h>
#include <vtkInteractorStyleTrackballCamera.h>
#include <vtkMath.h>
#include <vtkNamedColors.h>
#include <vtkNew.h>
#include <vtkOBBTree.h>
#include <vtkPolyData.h>
#include <vtkPolyDataMapper.h>
#include <vtkProperty.h>
#include <vtkProperty2D.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkRenderer.h>
#include <vtkSliderRepresentation2D.h>
#include <vtkSliderWidget.h>
#include <vtkSmartPointer.h>
#include <vtkSphereSource.h>
#include <vtkTextProperty.h>
#include <vtkBYUReader.h>
#include <vtkOBJReader.h>
#include <vtkPLYReader.h>
#include <vtkPolyDataReader.h>
#include <vtkSTLReader.h>
#include <vtkXMLPolyDataReader.h>
#include <vtksys/SystemTools.hxx>
namespace {
vtkSmartPointer<vtkPolyData> ReadPolyData(const char* fileName);
}
namespace {
class vtkSliderCallback : public vtkCommand
{
public:
static vtkSliderCallback* New()
{
return new vtkSliderCallback;
}
vtkSliderCallback() : OBBTree(0), Level(0), PolyData(0), Renderer(0)
{
}
virtual void Execute(vtkObject* caller, unsigned long, void*)
{
vtkSliderWidget* sliderWidget = reinterpret_cast<vtkSliderWidget*>(caller);
this->Level = vtkMath::Round(
static_cast<vtkSliderRepresentation*>(sliderWidget->GetRepresentation())
->GetValue());
this->OBBTree->GenerateRepresentation(this->Level, this->PolyData);
this->Renderer->Render();
}
vtkOBBTree* OBBTree;
int Level;
vtkPolyData* PolyData;
vtkRenderer* Renderer;
};
} // namespace
VisualizeOBBTree::VisualizeOBBTree(vtkRenderWindow* renderWin)
{
auto polyData = ReadPolyData("../resources/Armadillo.ply");
;
vtkNew<vtkNamedColors> colors;
vtkNew<vtkPolyDataMapper> pointsMapper;
pointsMapper->SetInputData(polyData);
pointsMapper->ScalarVisibilityOff();
vtkNew<vtkActor> pointsActor;
pointsActor->SetMapper(pointsMapper);
pointsActor->GetProperty()->SetInterpolationToFlat();
pointsActor->GetProperty()->SetColor(colors->GetColor4d("Yellow").GetData());
pointsActor->GetProperty()->SetOpacity(.3);
// Create the tree
vtkNew<vtkOBBTree> obbTree;
obbTree->SetDataSet(polyData);
obbTree->BuildLocator();
double corner[3] = { 0.0, 0.0, 0.0 };
double max[3] = { 0.0, 0.0, 0.0 };
double mid[3] = { 0.0, 0.0, 0.0 };
double min[3] = { 0.0, 0.0, 0.0 };
double size[3] = { 0.0, 0.0, 0.0 };
obbTree->ComputeOBB(polyData, corner, max, mid, min, size);
std::cout << "Corner:\t" << corner[0] << ", " << corner[1] << ", "
<< corner[2] << std::endl
<< "Max:\t" << max[0] << ", " << max[1] << ", " << max[2]
<< std::endl
<< "Mid:\t" << mid[0] << ", " << mid[1] << ", " << mid[2]
<< std::endl
<< "Min:\t" << min[0] << ", " << min[1] << ", " << min[2]
<< std::endl
<< "Size:\t" << size[0] << ", " << size[1] << ", " << size[2]
<< std::endl;
// Initialize the representation
vtkNew<vtkPolyData> polydata;
obbTree->GenerateRepresentation(0, polydata);
vtkNew<vtkPolyDataMapper> obbtreeMapper;
obbtreeMapper->SetInputData(polydata);
vtkNew<vtkActor> obbtreeActor;
obbtreeActor->SetMapper(obbtreeMapper);
obbtreeActor->GetProperty()->SetInterpolationToFlat();
obbtreeActor->GetProperty()->SetOpacity(.5);
obbtreeActor->GetProperty()->EdgeVisibilityOn();
obbtreeActor->GetProperty()->SetColor(
colors->GetColor4d("SpringGreen").GetData());
// A renderer and render window
vtkNew<vtkRenderer> renderer;
vtkNew<vtkRenderWindow> renderWindow;
renderWindow->AddRenderer(renderer);
// An interactor
vtkNew<vtkInteractorStyleTrackballCamera> style;
vtkNew<vtkRenderWindowInteractor> renderWindowInteractor;
renderWindowInteractor->SetRenderWindow(renderWindow);
renderWindowInteractor->SetInteractorStyle(style);
// Add the actors to the scene
renderer->AddActor(pointsActor);
renderer->AddActor(obbtreeActor);
renderer->SetBackground(colors->GetColor3d("MidnightBlue").GetData());
renderer->UseHiddenLineRemovalOn();
// Render an image (lights and cameras are created automatically)
renderWindow->SetWindowName("VisualizeOBBTree");
renderWindow->SetSize(600, 600);
renderWindow->Render();
vtkNew<vtkSliderRepresentation2D> sliderRep;
sliderRep->SetMinimumValue(0);
sliderRep->SetMaximumValue(obbTree->GetLevel());
sliderRep->SetValue(obbTree->GetLevel() / 2);
sliderRep->SetTitleText("Level");
sliderRep->GetPoint1Coordinate()->SetCoordinateSystemToNormalizedDisplay();
sliderRep->GetPoint1Coordinate()->SetValue(.2, .2);
sliderRep->GetPoint2Coordinate()->SetCoordinateSystemToNormalizedDisplay();
sliderRep->GetPoint2Coordinate()->SetValue(.8, .2);
sliderRep->SetSliderLength(0.075);
sliderRep->SetSliderWidth(0.05);
sliderRep->SetEndCapLength(0.05);
sliderRep->GetTitleProperty()->SetColor(
colors->GetColor3d("Beige").GetData());
sliderRep->GetCapProperty()->SetColor(
colors->GetColor3d("MistyRose").GetData());
sliderRep->GetSliderProperty()->SetColor(
colors->GetColor3d("LightBlue").GetData());
sliderRep->GetSelectedProperty()->SetColor(
colors->GetColor3d("Violet").GetData());
vtkNew<vtkSliderWidget> sliderWidget;
sliderWidget->SetInteractor(renderWindowInteractor);
sliderWidget->SetRepresentation(sliderRep);
sliderWidget->SetAnimationModeToAnimate();
sliderWidget->EnabledOn();
vtkNew<vtkSliderCallback> callback;
callback->OBBTree = obbTree;
callback->PolyData = polydata;
callback->Renderer = renderer;
callback->Execute(sliderWidget, 0, 0);
sliderWidget->AddObserver(vtkCommand::InteractionEvent, callback);
renderWindowInteractor->Initialize();
renderWindow->Render();
renderWindowInteractor->Start();
}
通过滑块调整包围盒Level可以得到不同包围细粒度的包围盒。Level有个最大值，可以通过obbTree->SetMaxLevel(maxLevel);设置，可以不不设置，ObbTree会自动设置最大Level，通过obbTree->GetLevel();函数可以获得最大Level。
1.4 对齐两个数据集
AlignTwoPolyDatas::AlignTwoPolyDatas(vtkRenderWindow* renderWin)
{
// Vis Pipeline
vtkNew<vtkNamedColors> colors;
vtkNew<vtkRenderer> renderer;
vtkNew<vtkRenderWindow> renderWindow;
renderWindow->SetSize(640, 480);
renderWindow->AddRenderer(renderer);
vtkNew<vtkRenderWindowInteractor> interactor;
interactor->SetRenderWindow(renderWindow);
renderer->SetBackground(colors->GetColor3d("sea_green_light").GetData());
renderer->UseHiddenLineRemovalOn();
auto sourcePolyData = ReadPolyData("../resources/Grey_Nurse_Shark.stl");
// Save the source polydata in case the align does not improve
// segmentation.
vtkNew<vtkPolyData> originalSourcePolyData;
originalSourcePolyData->DeepCopy(sourcePolyData);
//std::cout << "Loading target: " << argv[2] << std::endl;
auto targetPolyData = ReadPolyData("../resources/greatWhite.stl");
// If the target orientation is markedly different, you may need to apply a
// transform to orient the target with the source.
// For example, when using Grey_Nurse_Shark.stl as the
// source and thingiverse/greatWhite.stl as the target,
// you need to transform the target.
/*auto sourceFound =
std::string(argv[1]).find("Grey_Nurse_Shark.stl") != std::string::npos;
auto targetFound =
std::string(argv[2]).find("greatWhite.stl") != std::string::npos;*/
vtkNew<vtkTransform> trnf;
//if (sourceFound && targetFound)
{
trnf->RotateY(90);
}
vtkNew<vtkTransformPolyDataFilter> tpd;
tpd->SetTransform(trnf);
tpd->SetInputData(targetPolyData);
tpd->Update();
vtkNew<vtkHausdorffDistancePointSetFilter> distance;
distance->SetInputData(0, tpd->GetOutput());
distance->SetInputData(1, sourcePolyData);
distance->Update();
double distanceBeforeAlign = static_cast<vtkPointSet*>(distance->GetOutput(0))
->GetFieldData()
->GetArray("HausdorffDistance")
->GetComponent(0, 0);
// Get initial alignment using oriented bounding boxes.
AlignBoundingBoxes(sourcePolyData, tpd->GetOutput());
distance->SetInputData(0, tpd->GetOutput());
distance->SetInputData(1, sourcePolyData);
distance->Modified();
distance->Update();
double distanceAfterAlign = static_cast<vtkPointSet*>(distance->GetOutput(0))
->GetFieldData()
->GetArray("HausdorffDistance")
->GetComponent(0, 0);
double bestDistance = std::min(distanceBeforeAlign, distanceAfterAlign);
if (distanceAfterAlign > distanceBeforeAlign)
{
sourcePolyData->DeepCopy(originalSourcePolyData);
}
//Refine the alignment using IterativeClosestPoint.
vtkNew<vtkIterativeClosestPointTransform> icp;
icp->SetSource(sourcePolyData);
icp->SetTarget(tpd->GetOutput());
icp->GetLandmarkTransform()->SetModeToRigidBody();
icp->SetMaximumNumberOfLandmarks(100);
icp->SetMaximumMeanDistance(.00001);
icp->SetMaximumNumberOfIterations(500);
icp->CheckMeanDistanceOn();
icp->StartByMatchingCentroidsOn();
icp->Update();
auto icpMeanDistance = icp->GetMeanDistance();
// icp->Print(std::cout);
auto lmTransform = icp->GetLandmarkTransform();
vtkNew<vtkTransformPolyDataFilter> transform;
transform->SetInputData(sourcePolyData);
transform->SetTransform(lmTransform);
transform->SetTransform(icp);
transform->Update();
distance->SetInputData(0, tpd->GetOutput());
distance->SetInputData(1, transform->GetOutput());
distance->Update();
// Note: If there is an error extracting eigenfunctions, then this will be
// zero.
double distanceAfterICP = static_cast<vtkPointSet*>(distance->GetOutput(0))
->GetFieldData()
->GetArray("HausdorffDistance")
->GetComponent(0, 0);
if (!(std::isnan(icpMeanDistance) || std::isinf(icpMeanDistance)))
{
if (distanceAfterICP < bestDistance)
{
bestDistance = distanceAfterICP;
}
}
std::cout << "Distances:" << std::endl;
std::cout << " Before aligning: "
<< distanceBeforeAlign << std::endl;
std::cout << " Aligning using oriented bounding boxes: "
<< distanceAfterAlign << std::endl;
std::cout << " Aligning using IterativeClosestPoint: " << distanceAfterICP
<< std::endl;
std::cout << " Best distance: " << bestDistance
<< std::endl;
// Select the source to use.
vtkNew<vtkDataSetMapper> sourceMapper;
if (bestDistance == distanceBeforeAlign)
{
sourceMapper->SetInputData(originalSourcePolyData);
std::cout << "Using original alignment" << std::endl;
}
else if (bestDistance == distanceAfterAlign)
{
sourceMapper->SetInputData(sourcePolyData);
std::cout << "Using alignment by OBB" << std::endl;
}
else
{
sourceMapper->SetInputConnection(transform->GetOutputPort());
std::cout << "Using alignment by ICP" << std::endl;
}
sourceMapper->ScalarVisibilityOff();
vtkNew<vtkActor> sourceActor;
sourceActor->SetMapper(sourceMapper);
sourceActor->GetProperty()->SetOpacity(0.6);
sourceActor->GetProperty()->SetDiffuseColor(
colors->GetColor3d("White").GetData());
renderer->AddActor(sourceActor);
vtkNew<vtkDataSetMapper> targetMapper;
targetMapper->SetInputData(tpd->GetOutput());// tpd->GetOutput());
targetMapper->ScalarVisibilityOff();
vtkNew<vtkActor> targetActor;
targetActor->SetMapper(targetMapper);
targetActor->GetProperty()->SetDiffuseColor(
colors->GetColor3d("Tomato").GetData());
renderer->AddActor(targetActor);
renderWindow->AddRenderer(renderer);
renderWindow->SetWindowName("AlignTwoPolyDatas");
#if VTK_HAS_COW
vtkNew<vtkCameraOrientationWidget> camOrientManipulator;
camOrientManipulator->SetParentRenderer(renderer);
// Enable the widget.
camOrientManipulator->On();
#else
vtkNew<vtkAxesActor> axes;
vtkNew<vtkOrientationMarkerWidget> widget;
double rgba[4]{ 0.0, 0.0, 0.0, 0.0 };
colors->GetColor("Carrot", rgba);
widget->SetOutlineColor(rgba[0], rgba[1], rgba[2]);
widget->SetOrientationMarker(axes);
widget->SetInteractor(interactor);
widget->SetViewport(0.0, 0.0, 0.2, 0.2);
widget->EnabledOn();
widget->InteractiveOn();
#endif
renderWindow->Render();
interactor->Start();
}
void AlignBoundingBoxes(vtkPolyData* source, vtkPolyData* target)
{
// Use OBBTree to create an oriented bounding box for target and source
vtkNew<vtkOBBTree> sourceOBBTree;
sourceOBBTree->SetDataSet(source);
sourceOBBTree->SetMaxLevel(1);
sourceOBBTree->BuildLocator();
vtkNew<vtkOBBTree> targetOBBTree;
targetOBBTree->SetDataSet(target);
targetOBBTree->SetMaxLevel(1);
targetOBBTree->BuildLocator();
vtkNew<vtkPolyData> sourceLandmarks;
sourceOBBTree->GenerateRepresentation(0, sourceLandmarks);
vtkNew<vtkPolyData> targetLandmarks;
targetOBBTree->GenerateRepresentation(0, targetLandmarks);
vtkNew<vtkLandmarkTransform> lmTransform;
lmTransform->SetModeToSimilarity();
lmTransform->SetTargetLandmarks(targetLandmarks->GetPoints());
// vtkNew<vtkTransformPolyDataFilter> lmTransformPD;
double bestDistance = VTK_DOUBLE_MAX;
vtkNew<vtkPoints> bestPoints;
BestBoundingBox("X", target, source, targetLandmarks, sourceLandmarks,
bestDistance, bestPoints);
BestBoundingBox("Y", target, source, targetLandmarks, sourceLandmarks,
bestDistance, bestPoints);
BestBoundingBox("Z", target, source, targetLandmarks, sourceLandmarks,
bestDistance, bestPoints);
lmTransform->SetSourceLandmarks(bestPoints);
lmTransform->Modified();
vtkNew<vtkTransformPolyDataFilter> transformPD;
transformPD->SetInputData(source);
transformPD->SetTransform(lmTransform);
transformPD->Update();
source->DeepCopy(transformPD->GetOutput());
}
void BestBoundingBox(std::string const& axis, vtkPolyData* target,
vtkPolyData* source, vtkPolyData* targetLandmarks,
vtkPolyData* sourceLandmarks, double& bestDistance,
vtkPoints* bestPoints)
{
vtkNew<vtkHausdorffDistancePointSetFilter> distance;
vtkNew<vtkTransform> testTransform;
vtkNew<vtkTransformPolyDataFilter> testTransformPD;
vtkNew<vtkLandmarkTransform> lmTransform;
vtkNew<vtkTransformPolyDataFilter> lmTransformPD;
lmTransform->SetModeToSimilarity();
lmTransform->SetTargetLandmarks(targetLandmarks->GetPoints());
double sourceCenter[3];
sourceLandmarks->GetCenter(sourceCenter);
auto delta = 90.0;
for (auto i = 0; i < 4; ++i)
{
auto angle = delta * i;
// Rotate about center
testTransform->Identity();
testTransform->Translate(sourceCenter[0], sourceCenter[1], sourceCenter[2]);
if (axis == "X")
{
testTransform->RotateX(angle);
}
else if (axis == "Y")
{
testTransform->RotateY(angle);
}
else
{
testTransform->RotateZ(angle);
}
testTransform->Translate(-sourceCenter[0], -sourceCenter[1],
-sourceCenter[2]);
testTransformPD->SetTransform(testTransform);
testTransformPD->SetInputData(sourceLandmarks);
testTransformPD->Update();
lmTransform->SetSourceLandmarks(testTransformPD->GetOutput()->GetPoints());
lmTransform->Modified();
lmTransformPD->SetInputData(source);
lmTransformPD->SetTransform(lmTransform);
lmTransformPD->Update();
distance->SetInputData(0, target);
distance->SetInputData(1, lmTransformPD->GetOutput());
distance->Update();
double testDistance = static_cast<vtkPointSet*>(distance->GetOutput(0))
->GetFieldData()
->GetArray("HausdorffDistance")
->GetComponent(0, 0);
if (testDistance < bestDistance)
{
bestDistance = testDistance;
bestPoints->DeepCopy(testTransformPD->GetOutput()->GetPoints());
}
}
return;
}
对齐之前的两个模型：
对齐之后的两个模型：
对齐结果如下：
Distances:
Before aligning: 137.484
Aligning using oriented bounding boxes: 2.77965
Aligning using IterativeClosestPoint: 2.55923
Best distance: 2.55923
Using alignment by ICP
可以看到obbTree的方法跟ICP相差不大。
1.5 圆柱形有向包围盒
#include <vtkActor.h>
#include <vtkCleanPolyData.h>
#include <vtkExtractEnclosedPoints.h>
#include <vtkLineSource.h>
#include <vtkNamedColors.h>
#include <vtkNew.h>
#include <vtkOBBTree.h>
#include <vtkPoints.h>
#include <vtkPolyDataMapper.h>
#include <vtkProperty.h>
#include <vtkQuad.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkRenderer.h>
#include <vtkSmartPointer.h>
#include <vtkTubeFilter.h>
// Readers
#include <vtkBYUReader.h>
#include <vtkOBJReader.h>
#include <vtkPLYReader.h>
#include <vtkPolyDataReader.h>
#include <vtkSTLReader.h>
#include <vtkXMLPolyDataReader.h>
#include <vtkPolyData.h>
#include <vtkSphereSource.h>
#include <array>
#include <cctype> // For to_lower
#include <cstdlib>
#include <iostream>
#include <string> // For find_last_of()
namespace {
vtkSmartPointer<vtkPolyData> ReadPolyData(std::string const& fileName);
double MakeAQuad(std::vector<std::array<double, 3>>&, std::array<double, 3>&);
} // namespace
OrientedBoundingCylinder::OrientedBoundingCylinder(vtkRenderWindow* renderWin)
{
auto polyData = ReadPolyData("../resources/Armadillo.ply");
// Get bounds of polydata.
std::array<double, 6> bounds;
polyData->GetBounds(bounds.data());
// Create the tree.
vtkNew<vtkOBBTree> obbTree;
obbTree->SetDataSet(polyData);
obbTree->SetMaxLevel(1);
obbTree->BuildLocator();
// Get the PolyData for the OBB.
vtkNew<vtkPolyData> obbPolydata;
obbTree->GenerateRepresentation(0, obbPolydata);
// Get the points of the OBB.
vtkNew<vtkPoints> obbPoints;
obbPoints->DeepCopy(obbPolydata->GetPoints());
// Use a quad to find centers of OBB faces.
vtkNew<vtkQuad> aQuad;
std::vector<std::array<double, 3>> facePoints(4);
std::vector<std::array<double, 3>> centers(3);
std::vector<std::array<double, 3>> endPoints(3);
std::array<double, 3> center;
std::array<double, 3> endPoint;
std::array<double, 3> point0, point1, point2, point3, point4, point5, point6,
point7;
std::array<double, 3> radii;
std::array<double, 3> lengths;
// Transfer the points to std::array's.
obbPoints->GetPoint(0, point0.data());
obbPoints->GetPoint(1, point1.data());
obbPoints->GetPoint(2, point2.data());
obbPoints->GetPoint(3, point3.data());
obbPoints->GetPoint(4, point4.data());
obbPoints->GetPoint(5, point5.data());
obbPoints->GetPoint(6, point6.data());
obbPoints->GetPoint(7, point7.data());
// x face.
// ids[0] = 2; ids[1] = 3; ids[2] = 7; ids[3] = 6;
facePoints[0] = point2;
facePoints[1] = point3;
facePoints[2] = point7;
facePoints[3] = point6;
radii[0] = MakeAQuad(facePoints, centers[0]);
MakeAQuad(facePoints, centers[0]);
// ids[0] = 0; ids[1] = 4; ids[2] = 5; ids[3] = 1;
facePoints[0] = point0;
facePoints[1] = point4;
facePoints[2] = point5;
facePoints[3] = point1;
MakeAQuad(facePoints, endPoints[0]);
lengths[0] = std::sqrt(vtkMath::Distance2BetweenPoints(centers[0].data(),
endPoints[0].data())) /
2.0;
// y face.
// ids[0] = 0; ids[1] = 1; ids[2] = 2; ids[3] = 3;
facePoints[0] = point0;
facePoints[1] = point1;
facePoints[2] = point2;
facePoints[3] = point3;
radii[1] = MakeAQuad(facePoints, centers[1]);
// ids[0] = 4; ids[1] = 6; ids[2] = 7; ids[3] = 5;
facePoints[0] = point4;
facePoints[1] = point6;
facePoints[2] = point7;
facePoints[3] = point5;
MakeAQuad(facePoints, endPoints[1]);
lengths[1] = std::sqrt(vtkMath::Distance2BetweenPoints(centers[1].data(),
endPoints[1].data())) /
2.0;
// z face.
// ids[0] = 0; ids[1] = 2; ids[2] = 6; ids[3] = 4;
facePoints[0] = point0;
facePoints[1] = point2;
facePoints[2] = point6;
facePoints[3] = point4;
MakeAQuad(facePoints, centers[2]);
radii[2] =
std::sqrt(vtkMath::Distance2BetweenPoints(point0.data(), point2.data())) /
2.0;
double outerRadius =
std::sqrt(vtkMath::Distance2BetweenPoints(point0.data(), point6.data())) /
2.0;
// ids[0] = 1; ids[1] = 3; ids[2] = 7; ids[3] = 5;
facePoints[0] = point1;
facePoints[1] = point5;
facePoints[2] = point7;
facePoints[3] = point3;
MakeAQuad(facePoints, endPoints[2]);
lengths[2] = std::sqrt(vtkMath::Distance2BetweenPoints(centers[2].data(),
endPoints[2].data())) /
2.0;
// Find long axis.
int longAxis = -1;
double length = VTK_DOUBLE_MIN;
for (auto i = 0; i < 3; ++i)
{
std::cout << "length: " << lengths[i] << std::endl;
if (lengths[i] > length)
{
length = lengths[i];
longAxis = i;
}
}
if (longAxis < 0)
{
std::cout << "There is no long axis" << std::endl;
return;
}
std::cout << "longAxis: " << longAxis << std::endl;
std::cout << "radii: ";
double radius = radii[longAxis];
for (const auto& a : radii)
{
std::cout << a << ", ";
}
std::cout << std::endl;
std::cout << "radius: " << radius << std::endl;
std::cout << "outerRadius: " << outerRadius << std::endl;
center = centers[longAxis];
endPoint = endPoints[longAxis];
vtkNew<vtkNamedColors> colors;
vtkNew<vtkLineSource> lineSource;
lineSource->SetPoint1(center.data());
lineSource->SetPoint2(endPoint.data());
vtkNew<vtkTubeFilter> tube;
tube->SetInputConnection(lineSource->GetOutputPort());
tube->SetRadius(radius);
tube->SetNumberOfSides(51);
tube->CappingOn();
tube->Update();
// See if all points lie inside cylinder.
vtkNew<vtkCleanPolyData> clean;
clean->SetInputData(tube->GetOutput());
clean->Update();
vtkNew<vtkExtractEnclosedPoints> enclosedPoints;
enclosedPoints->SetSurfaceData(clean->GetOutput());
enclosedPoints->SetInputData(polyData);
enclosedPoints->SetTolerance(.0001);
enclosedPoints->GenerateOutliersOn();
enclosedPoints->CheckSurfaceOn();
enclosedPoints->Update();
std::cout << "polyData points: " << polyData->GetPoints()->GetNumberOfPoints()
<< " excluded points: "
<< enclosedPoints->GetOutput(1)->GetPoints()->GetNumberOfPoints()
<< std::endl;
vtkNew<vtkPolyDataMapper> repMapper;
repMapper->SetInputData(obbPolydata);
vtkNew<vtkActor> repActor;
repActor->SetMapper(repMapper);
repActor->GetProperty()->SetColor(colors->GetColor3d("peacock").GetData());
repActor->GetProperty()->SetOpacity(.6);
// Create a mapper and actor for the cylinder.
vtkNew<vtkPolyDataMapper> cylinderMapper;
cylinderMapper->SetInputConnection(tube->GetOutputPort());
vtkNew<vtkActor> cylinderActor;
cylinderActor->SetMapper(cylinderMapper);
cylinderActor->GetProperty()->SetColor(
colors->GetColor3d("banana").GetData());
cylinderActor->GetProperty()->SetOpacity(.5);
vtkNew<vtkPolyDataMapper> originalMapper;
originalMapper->SetInputData(polyData);
vtkNew<vtkActor> originalActor;
originalActor->SetMapper(originalMapper);
originalActor->GetProperty()->SetColor(
colors->GetColor3d("tomato").GetData());
// Create a renderer, render window, and interactor.
vtkNew<vtkRenderer> renderer;
renderer->UseHiddenLineRemovalOn();
// Display all centers and endpoints.
std::vector<vtkColor3d> cs;
cs.push_back(colors->GetColor3d("red"));
cs.push_back(colors->GetColor3d("green"));
cs.push_back(colors->GetColor3d("blue"));
for (auto i = 0; i < 3; ++i)
{
vtkNew<vtkSphereSource> ps1;
ps1->SetCenter(centers[i].data());
ps1->SetRadius(length * .04);
ps1->SetPhiResolution(21);
ps1->SetThetaResolution(41);
vtkNew<vtkPolyDataMapper> pm1;
pm1->SetInputConnection(ps1->GetOutputPort());
vtkNew<vtkActor> pa1;
pa1->GetProperty()->SetColor(cs[i].GetData());
pa1->GetProperty()->SetSpecularPower(50);
pa1->GetProperty()->SetSpecular(.4);
pa1->GetProperty()->SetDiffuse(.6);
pa1->SetMapper(pm1);
renderer->AddActor(pa1);
vtkNew<vtkSphereSource> ps2;
ps2->SetCenter(endPoints[i].data());
ps2->SetRadius(length * .04);
ps2->SetPhiResolution(21);
ps2->SetThetaResolution(41);
vtkNew<vtkPolyDataMapper> pm2;
pm2->SetInputConnection(ps2->GetOutputPort());
vtkNew<vtkActor> pa2;
pa2->GetProperty()->SetColor(cs[i].GetData());
pa2->SetMapper(pm2);
renderer->AddActor(pa2);
}
vtkNew<vtkRenderWindow> renderWindow;
renderWindow->AddRenderer(renderer);
renderWindow->SetWindowName("OrientedBoundingCylinder");
renderWindow->SetSize(640, 480);
vtkNew<vtkRenderWindowInteractor> renderWindowInteractor;
renderWindowInteractor->SetRenderWindow(renderWindow);
// Add the actors to the scene.
renderer->AddActor(originalActor);
renderer->AddActor(cylinderActor);
// renderer->AddActor(repActor);
renderer->GradientBackgroundOn();
renderer->SetBackground2(colors->GetColor3d("SkyBlue").GetData());
renderer->SetBackground(colors->GetColor3d("LightSeaGreen").GetData());
// double adjustedRadius = radius;
double adjustedIncr = (outerRadius - radius) / 20.0;
if (enclosedPoints->GetOutput(1)->GetPoints()->GetNumberOfPoints() > 4)
{
std::cout << "improving..." << std::endl;
for (double r = radius;
enclosedPoints->GetOutput(1)->GetPoints()->GetNumberOfPoints() > 4;
r += adjustedIncr)
{
tube->SetRadius(r);
tube->Update();
clean->Update();
enclosedPoints->Update();
if (enclosedPoints->GetOutput(1)->GetPoints() != nullptr)
{
std::cout << "r: " << r << std::endl;
std::cout
<< " excluded points: "
<< enclosedPoints->GetOutput(1)->GetPoints()->GetNumberOfPoints()
<< std::endl;
renderWindow->Render();
}
else
{
break;
}
}
}
// Render and interact.
renderWindowInteractor->Start();
}
首先建立OBB树，然后找到上下前后左右每个面的中心和对应的半径，然后找到中心距离最远的轴线作为圆柱的轴线，通过该轴线和对应的半径形成一个圆柱，然后判断是否所有点都在圆柱内部，如果有点没被包含进包围盒则递增圆柱的半径，直到所有点都在包围盒内。
扫一扫
810
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
