I want to find the boundary of all projected points which are computed from the points on the 3D model. The red arrow in the following image represents the projected vector. I will project the original points on the 3D model to a plane.

Scan the projected points with a series of straight lines in the plane to find all points on the boundary.
Find the two farthest points to generate vector Vec0, then we can calculate the cross vector of Vec0 and project vector, I name the result vector Vec1. Move the vector Vec0 along Vec1 to find all intersecting points.

    \[P=(x_0, y_0, z_0)\]

    \[Vec_0 = (a, b, c)\]

    \[\frac{ x - x_0 }{ a } = \frac{ y - y_0 }{ b } = \frac{ z - z_0 }{ c }\]

We can convert the equation to the multiplication form.

    \[bc (x - x_0) = ac ( y - y_0 ) = ab (z - z_0)\]

If a projected point coordinate value satisfies this equation, we can regard it as an intersecting point.
Because the projected points have discrete distribution, so I allow computing error in the comparison process.

#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 <vtkSTLReader.h>
#include <vtkPlane.h>
#include <vtkProperty.h>
#include "point.hpp"
#define vtkSPtr vtkSmartPointer
#define vtkSPtrNew(Var, Type) vtkSPtr<Type> Var = vtkSPtr<Type>::New();
using namespace std;
int main() {
    vtkSPtrNew( reader, vtkSTLReader );
    reader->SetFileName( "C:\\Users\\Administrator\\Desktop\\bottomTest.stl" );
    reader->Update();
    vtkPolyData *data = reader->GetOutput();
    vtkSPtrNew( mapper, vtkPolyDataMapper );
    mapper->SetInputData( reader->GetOutput() );
    vtkSPtrNew( actor, vtkActor );
    actor->SetMapper( mapper );
    vtkSPtrNew( renderer, vtkRenderer );
    renderer->AddActor(actor);
    renderer->SetBackground( 0, 0, 0 );
    vtkSPtrNew( renderWindow, vtkRenderWindow );
    renderWindow->AddRenderer( renderer );
    vtkSPtrNew( renderWindowInteractor, vtkRenderWindowInteractor );
    renderWindowInteractor->SetRenderWindow( renderWindow );
    // ============ start to project points =================
    Point proOrigin( 0, 0, 5 );
    Point proNormal( 0.0593944, 8.18427e-05, -0.997451 );
    vtkSPtrNew( proPlane, vtkPlane );
    proPlane->SetOrigin( proOrigin.point );
    proPlane->SetNormal( proNormal.point );
    vtkSPtrNew( proPoints, vtkPoints );
    for( int i = 0; i < data->GetNumberOfPoints(); ++i )
    {
        Point pt( data->GetPoint( i ) );
        proPlane->ProjectPoint( pt.point, pt.point );
        proPoints->InsertNextPoint( pt.point );
    }
    // ==================== finished ========================
    // ============ generate the fixed point and move vector =================
    Point comVec;
    Point pt0( proPoints->GetPoint(0) );
    Point pt1 = pt0;
    Point fixedPoint = pt0;
    double maxDis = VTK_DOUBLE_MIN;
    for( int i = 0; i < proPoints->GetNumberOfPoints(); ++i )
    {
        for( int j = i+1; j < proPoints->GetNumberOfPoints(); ++j )
        {
            Point p0( proPoints->GetPoint( i ) );
            Point p1( proPoints->GetPoint( j ) );
            double tmpDis = (p0 - p1).Length();
            if( tmpDis > maxDis )
            {
                maxDis = tmpDis;
                comVec = p0 - p1;
                fixedPoint = p0;
            }
        }
    }
    comVec.Unit();
    Point moveVec = comVec ^ proNormal;
    moveVec.Unit();
    // ==================== finished ========================
    std::vector<Point> boundaryPts;
    // ==================== positive search ==================
    double limit = 1e-2;
    double delta = 0.001;
    int count = 10000;
    int step = -count;
    Point comPt;
    std::vector<Point> linePts;
    Point midPt( 0, 0, 0 );
    while ( step < count )
    {
        // find next line
        comPt = fixedPoint + moveVec * delta * (step++);
        linePts.clear();
        midPt = Point(0, 0, 0);
        for( int i = 0; i < proPoints->GetNumberOfPoints(); ++i )
        {
            Point tmpPt( proPoints->GetPoint( i ) );
            double v0 = (tmpPt[0] - comPt[0]) * comVec[1] *comVec[2];
            double v1 = (tmpPt[1] - comPt[1]) * comVec[0] *comVec[2];
            double v2 = (tmpPt[2] - comPt[2]) * comVec[0] *comVec[1];
            if( fabs( v0 - v1 ) < limit && fabs( v1 - v2 ) < limit )
            {
                linePts.push_back( tmpPt );
                midPt = midPt + tmpPt;
                // {
                //     vtkSPtrNew( sphereSource, vtkSphereSource );
                //     sphereSource->SetCenter( tmpPt.point );
                //     sphereSource->SetRadius( 0.05 );
                //     sphereSource->Update();
                //     vtkSPtrNew( mapper, vtkPolyDataMapper );
                //     mapper->SetInputData( sphereSource->GetOutput() );
                //     mapper->Update();
                //     vtkSPtrNew( ptActor, vtkActor );
                //     ptActor->SetMapper( mapper );
                //     ptActor->GetProperty()->SetColor( 1, 0, 0 );
                //     renderer->AddActor( ptActor );
                // }
            }
        }
        cout << "linePts size: " << linePts.size() << endl;
        if( linePts.size() > 0 )
        {
            midPt /= linePts.size();
            double maxValue = VTK_DOUBLE_MIN;
            double minValue = VTK_DOUBLE_MAX;
            Point endPts[2];
            for( int i = 0; i < linePts.size(); ++i )
            {
                Point vec = linePts[i] - midPt;
                double dotValue = vec.Dot( comVec );
                if( maxValue < dotValue )
                {
                    maxValue = dotValue;
                    endPts[1] = linePts[i];
                }
                if( minValue > dotValue )
                {
                    minValue = dotValue;
                    endPts[0] = linePts[i];
                }
            }
            boundaryPts.push_back( endPts[0] );
            boundaryPts.push_back( endPts[1] );
            if( endPts[0] == endPts[1] )
            {
                boundaryPts.pop_back();
            }
        }
    }
    // ==================== finished =======================
    cout << "boundaryPts.size: " << boundaryPts.size() << endl;
    for( int i = 0; i < boundaryPts.size(); ++i )
    {
        cout << boundaryPts[i][0] << ", " << boundaryPts[i][1] << ", " << boundaryPts[i][2] << endl;
        vtkSPtrNew( sphereSource, vtkSphereSource );
        sphereSource->SetCenter( boundaryPts[i].point );
        sphereSource->SetRadius( 0.01 );
        sphereSource->Update();
        vtkSPtrNew( mapper, vtkPolyDataMapper );
        mapper->SetInputData( sphereSource->GetOutput() );
        mapper->Update();
        vtkSPtrNew( ptActor, vtkActor );
        ptActor->SetMapper( mapper );
        ptActor->GetProperty()->SetColor( 1, 1, 0 );
        renderer->AddActor( ptActor );
    }
    renderer->ResetCamera();
    renderWindow->Render();
    renderWindowInteractor->Start();
    return 0;
}

Categories: VTK

0 0 votes
Article Rating
Subscribe
Notify of
guest
0 Comments
Inline Feedbacks
View all comments

3D Model Viewer: add grid plane and convex hull.

X
A prohibited operation
0
Would love your thoughts, please comment.x
()
x