Here is a curve in the 3D scene, the sphere represents its center.
I want to find the axes of the 3D model just as vtkOBBTree, relative article:
Explore The Orientation Of 3D Object.
But all points on the curve are in the same plane, so neither vtkOBBTree nor vtkHull can help us.
The points look just like
470.903412 875.383850 0.000000
471.292023 875.435852 0.000000
471.680634 875.487122 0.000000
472.069305 875.537659 0.000000
472.457977 875.587402 0.000000
I closed the curve and put a line which goes through the central sphere on the scene, there are two intersecting points on the curve, then we can compute a distance between them. Add another line which is orthogonal to the first line. Finally, we can get four intersect points A, B, C, and D as the following image.
Lets’ rotate the line around the center point until the difference between AB’s length and CD’s length becomes the biggest.
The result shows height and width vectors as the above image. We have to project the points on the curve to the width vector and find the end points if we want to get the real width.
Here are all the implementation details about the above logic.
#include <iostream>
#include <vtkSmartPointer.h>
#include <vtkSphereSource.h>
#include <vtkActor.h>
#include <vtkConeSource.h>
#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
#include <vtkPolyDataMapper.h>
#include <vtkXMLPolyDataReader.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkXMLPolyDataWriter.h>
#include <vtkHull.h>
#include <vtkProperty.h>
#include "Point.hpp"
#define vtkSPtr vtkSmartPointer
#define vtkSPtrNew(Var, Type) vtkSPtr<Type> Var = vtkSPtr<Type>::New();
using namespace std;
vtkPolyData *curvePd = nullptr;
vtkSmartPointer<vtkRenderer> renderer;
void ShowPoint( Point pt, vtkRenderer *renderer, double *color )
{
vtkSmartPointer<vtkSphereSource> sphere =
vtkSmartPointer<vtkSphereSource>::New();
sphere->SetThetaResolution( 100 );
sphere->SetPhiResolution( 50 );
sphere->SetRadius( 10 );
sphere->SetCenter( pt.point );
sphere->Update();
vtkSmartPointer<vtkPolyDataMapper> sphereMapper =
vtkSmartPointer<vtkPolyDataMapper>::New();
sphereMapper->SetInputData( sphere->GetOutput() );
vtkSmartPointer<vtkActor> actor = vtkSmartPointer<vtkActor>::New();
actor->SetMapper( sphereMapper );
actor->GetProperty()->SetColor( color );
renderer->AddActor( actor );
}
double CalculateDisFromPtToLine(Point pt, Point l1, Point l2)
{
Point line = l1 - l2;
line.Unit();
Point vec = pt - l1;
double dotValue = vec.Dot( line );
Point projectedVec = line*dotValue;
Point final = -projectedVec + vec;
return final.Length();
}
double GetDistanceForLineIntersect( Point center, Point dir, Point &resultClosetPt1, Point &resultClosetPt2 )
{
int projectLen = 10;
Point linePts[2] = { center+projectLen*dir, center-projectLen*dir };
Point negDir = linePts[1] - center;
Point closedPt1( 0, 0, 0 ), closedPt2( 0, 0, 0 );
double closedDis1 = VTK_DOUBLE_MAX, closedDis2 = VTK_DOUBLE_MAX;
for( int i = 0; i < curvePd->GetNumberOfPoints(); ++i )
{
Point pt( curvePd->GetPoint( i ) );
Point vec = pt - center;
double dis = CalculateDisFromPtToLine( pt, linePts[0], linePts[1] );
if( vec.Dot( negDir ) < 0 )
{
if( closedDis1 > dis )
{
closedDis1 = dis;
closedPt1 = pt;
}
}
else {
if( closedDis2 > dis )
{
closedDis2 = dis;
closedPt2 = pt;
}
}
}
cout << closedDis1 << " vs " << closedDis2 << endl;
resultClosetPt1 = closedPt1;
resultClosetPt2 = closedPt2;
return (closedPt1 - closedPt2).Length();
}
Point GetProjectPtOnDir(Point pt, Point dir, Point startPt)
{
dir.Unit();
Point ptVector = pt - startPt;
double cosTheta = ptVector.Dot( dir );
Point resultPt = startPt + dir * cosTheta;
return resultPt;
}
void FindSizeInfo( Point center )
{
double maxDis = VTK_DOUBLE_MIN, minDis = VTK_DOUBLE_MAX;
Point endPt[4];
double maxDiff = 0;
for( double angle = 0; angle < 180; angle = angle+0.5 )
{
double radians1 = vtkMath::RadiansFromDegrees( angle );
double radians2 = vtkMath::RadiansFromDegrees( angle+90 );
Point dir1( std::cos( radians1 ), std::sin( radians1 ), 0 );
Point dir2( std::cos( radians2 ), std::sin( radians2 ), 0 );
Point resultClosetPt1[2], resultClosetPt2[2];
double dis1 = GetDistanceForLineIntersect( center, dir1, resultClosetPt1[0], resultClosetPt1[1] );
double dis2 = GetDistanceForLineIntersect( center, dir2, resultClosetPt2[0], resultClosetPt2[1] );
double diff = fabs(dis1 - dis2);
if( maxDiff < diff )
{
endPt[0] = resultClosetPt1[0];
endPt[1] = resultClosetPt1[1];
endPt[2] = resultClosetPt2[0];
endPt[3] = resultClosetPt2[1];
}
}
// ------------ debug -------------
/*double color[3] = { 1, 0, 0 };
ShowPoint( endPt[0], renderer, color );
ShowPoint( endPt[1], renderer, color );
color[2] = 1;
ShowPoint( endPt[2], renderer, color );
ShowPoint( endPt[3], renderer, color ); */
// ------------- calculate size -----------
Point dir[2] = { endPt[1] - endPt[0], endPt[3] - endPt[2] };
double negVecLen = 0, posVecLen = 0;
Point edge1[2], edge2[2];
for( int i = 0; i < curvePd->GetNumberOfPoints(); ++i )
{
Point pt( curvePd->GetPoint( i ) );
Point projectPt = GetProjectPtOnDir( pt, dir[0], center );
double dotValue = (projectPt - center).Dot( dir[0] );
double vecLen = (projectPt - center).Length();
if( dotValue > 0 )
{
if( posVecLen < vecLen )
{
posVecLen = vecLen;
edge1[0] = projectPt;
}
}
else
{
if( negVecLen < vecLen )
{
negVecLen = vecLen;
edge1[1] = projectPt;
}
}
}
negVecLen = posVecLen = 0;
for( int i = 0; i < curvePd->GetNumberOfPoints(); ++i )
{
Point pt( curvePd->GetPoint( i ) );
Point projectPt = GetProjectPtOnDir( pt, dir[1], center );
double dotValue = (projectPt - center).Dot( dir[0] );
double vecLen = (projectPt - center).Length();
if( dotValue > 0 )
{
if( posVecLen < vecLen )
{
posVecLen = vecLen;
edge2[0] = projectPt;
}
}
else
{
if( negVecLen < vecLen )
{
negVecLen = vecLen;
edge2[1] = projectPt;
}
}
}
double color[3] = { 1, 0, 0 };
ShowPoint( edge1[0], renderer, color );
ShowPoint( edge1[1], renderer, color );
color[2] = 1;
ShowPoint( edge2[0], renderer, color );
ShowPoint( edge2[1], renderer, color );
}
int main()
{
// we can't use vtkOBBTree or vtkHull to analyze the 2D curve through it is represented as 3D model.
vtkSPtrNew( xmlReader, vtkXMLPolyDataReader );
xmlReader->SetFileName( "./scurve.vtp" );
xmlReader->Update();
curvePd = xmlReader->GetOutput();
std::cout << curvePd->GetNumberOfCells() << std::endl;
std::cout << curvePd->GetNumberOfPoints() << std::endl;
vtkSPtrNew( mapper, vtkPolyDataMapper );
mapper->SetInputData( curvePd );
vtkSPtrNew( actor, vtkActor );
actor->SetMapper( mapper );
renderer = vtkSmartPointer<vtkRenderer>::New();
renderer->AddActor( actor );
renderer->SetBackground( 0, 0, 0 );
// --------------- start to close the curve --------------------
curvePd->BuildCells();
curvePd->BuildLinks();
/*for( int i = 0; i < curvePd->GetNumberOfPoints(); ++i )
{
Point pt( curvePd->GetPoint( i ) );
printf( "%lf\t%lf\t%lf\n", pt[0], pt[1], pt[2] );
} */
int curvePtCount = curvePd->GetNumberOfPoints();
Point firstPt( curvePd->GetPoint( 0 ) );
Point lastPt( curvePd->GetPoint( curvePtCount - 1 ) );
Point dir = lastPt - firstPt;
double dis = dir.Length();
dir.Unit();
// insert points and connect them.
int newPtsCount = 99;
for( int i = 0; i < newPtsCount; i++ )
{
double len = 1.0*(i+1)/(newPtsCount+1)*dis;
Point tmpPt = firstPt + dir * len;
curvePd->InsertNextLinkedPoint( tmpPt.point, 1 );
}
vtkIdType line0[2] = { 0, curvePtCount-1+1 };
curvePd->InsertNextLinkedCell( VTK_POLY_LINE, 2, line0 );
for( int i = 1; i < newPtsCount; i++ )
{
vtkIdType line[2] = { curvePtCount-1+i, curvePtCount-1+i+1 };
curvePd->InsertNextLinkedCell( VTK_POLY_LINE, 2, line );
}
vtkIdType line1[2] = { curvePtCount-1+newPtsCount, curvePtCount-1 };
curvePd->InsertNextLinkedCell( VTK_POLY_LINE, 2, line1 );
curvePd->Modified();
// ----------------- explore the line intersections -------------------
Point center( actor->GetCenter() );
double color[3] = { 0, 1, 1 };
ShowPoint( center, renderer, color );
FindSizeInfo( center );
vtkSPtrNew( renderWindow, vtkRenderWindow );
renderWindow->AddRenderer( renderer );
vtkSPtrNew( renderWindowInteractor, vtkRenderWindowInteractor );
renderWindowInteractor->SetRenderWindow( renderWindow );
renderer->ResetCamera();
renderWindow->Render();
renderWindowInteractor->Start();
return 0;
}
Result:
