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 <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->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 );

}

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 )
{
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.

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->SetBackground( 0, 0, 0 );

// --------------- start to close the curve --------------------
curvePd->BuildCells();
/*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;
}

vtkIdType line0[2] = { 0, curvePtCount-1+1 };
for( int i = 1; i < newPtsCount; i++ )
{
vtkIdType line[2] = { curvePtCount-1+i, curvePtCount-1+i+1 };
}
vtkIdType line1[2] = { curvePtCount-1+newPtsCount, curvePtCount-1 };
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 );

vtkSPtrNew( renderWindowInteractor, vtkRenderWindowInteractor );
renderWindowInteractor->SetRenderWindow( renderWindow );

renderer->ResetCamera();
renderWindow->Render();

renderWindowInteractor->Start();
return 0;
}

Result:

Categories: VTK

Article Rating
Subscribe
Notify of