Find The Center Of Discrete Points & Get Projected Vector On Direction

I want to find the central point, not mass center, between a few points like the following image.
The points are on the same plane.

Find the longest axis and the orthogonal axis for the points set firstly. Compute the positive farthest point and negative farthest point for every axis. Finally, we can adjust the mass center of points by the four vectors, the new center is what we want.

#include <iostream>
#include <vector>

#include <iostream>
#include <vtkPolyData.h>
#include <vtkProperty.h>
#include <vtkPolyDataMapper.h>
#include <vtkActor.h>
#include <vtkPlane.h>
#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkPoints.h>
#include <vtkInteractorStyleTrackballCamera.h>
#include <vtkOBBTree.h>

#include "./point.hpp"
#include <vtkPlane.h>
#include <vtkPlaneSource.h>
#include <vtkTransform.h>
#include <vtkDelaunay3D.h>
#include <vtkDataSetSurfaceFilter.h>

#define vtkSPtr vtkSmartPointer
#define vtkSPtrNew(Var, Type) vtkSPtr<Type> Var = vtkSPtr<Type>::New();

Point GetProjectVecOnDir(Point pt, Point dir, Point startPt)
{
    dir.Unit();
    Point ptVector = pt - startPt;
    double cosTheta = ptVector.Dot( dir );
    return (dir * cosTheta);
}

Point GetCenterOfPoints(std::vector<Point> pts, Point standardVec0, Point standardVec1)
{
    standardVec0.Unit();
    standardVec1.Unit();
    Point massCenter( 0, 0, 0 );
    for( auto pt: pts )
    {
        massCenter = massCenter + pt;
    }
    massCenter /= pts.size();

    // ------------------- calculate move vector in standardVec0 -------------------
    Point maxVec;
    Point minVec;
    double maxVecValue = VTK_DOUBLE_MIN;
    double minVecValue = VTK_DOUBLE_MAX;
    for( auto pt: pts )
    {
        auto vec = GetProjectVecOnDir( pt, standardVec0, massCenter );
        auto value = vec.Dot( standardVec0 );
        if( maxVecValue < value )
        {
            maxVecValue = value;
            maxVec = vec;
        }
        if( minVecValue > value )
        {
            minVecValue = value;
            minVec = vec;
        }
    }
    Point adjustVec0 = maxVec + minVec;
    adjustVec0 /= 2;

    // ------------------- calculate move vector in standardVec1 -------------------
    maxVecValue = VTK_DOUBLE_MIN;
    minVecValue = VTK_DOUBLE_MAX;
    for( auto pt: pts )
    {
        auto vec = GetProjectVecOnDir( pt, standardVec1, massCenter );
        auto value = vec.Dot( standardVec1 );
        if( maxVecValue < value )
        {
            maxVecValue = value;
            maxVec = vec;
        }
        if( minVecValue > value )
        {
            minVecValue = value;
            minVec = vec;
        }
    }
    Point adjustVec1 = maxVec + minVec;
    adjustVec1 /= 2;

    return (massCenter + adjustVec0 + adjustVec1);
}

int main() {
    Point normal( -0.0170471, 0.98981, -0.141371 );
    std::vector<Point> samplePts{ {-25.8979, -4.12537, -13.2241},
                                  {-19.1172, -4.24538, -14.882},
                                  {-25.9613, -3.41124, -8.21646},
                                  {-19.7642, -3.57102, -10.0824},
                                  {-23.3579, -2.62887, -3.05267},
                                  {-17.4944, -2.82134, -5.10724},
                                  {-22.9702, -1.96459, 1.55158},
                                  {-17.1261, -2.20013, -0.802314},
                                  {-17.8337, -0.932447, 8.15871},
                                  {-19.5861, -1.61249, 3.60869},
                                  {-22.3959, -1.17335, 7.02217},
                                  {-15.4835, -1.32562, 5.12249},
                                  {-14.2595, -0.103165, 13.5339},
                                  {-17.1115, -0.729489, 9.49264},
                                  {-18.1866, -0.195735, 13.3593},
                                  {-13.1472, -0.472227, 10.8158},
                                  {-12.6442, 1.07117, 21.5613},
                                  {-15.602, 0.497715, 17.9029},
                                  {-14.9051, 0.909782, 20.7039},
                                  {-7.35411, 1.78261, 25.9045},
                                  {-10.9231, 1.34243, 23.253},
                                  {-9.81953, 1.66037, 25.3459},
                                  {0.460095, 2.21861, 28.0149},
                                  {-5.69037, 2.04031, 27.5081},
                                  {-2.63058, 2.27101, 28.7545},
                                  {8.11476, 2.37368, 28.1776},
                                  {2.30691, 2.21672, 27.779},
                                  {5.25961, 2.43996, 28.9859},
                                  {14.1472, 2.00728, 24.8848},
                                  {10.7012, 2.32136, 27.4994},
                                  {13.1834, 2.29874, 27.0417},
                                  {21.2462, 1.11331, 17.7697},
                                  {17.1191, 1.41558, 20.3837},
                                  {19.8277, 1.36909, 19.7316},
                                  {21.1534, -0.140178, 9.00458},
                                  {19.4958, 0.4626, 13.4248},
                                  {22.9479, 0.425237, 12.7469},
                                  {17.9374, 0.131773, 11.2964},
                                  {32.7624, -2.96318, -12.1605},
                                  {26.8703, -3.32433, -13.9786},
                                  {31.6866, -2.23754, -6.9502},
                                  {26.3174, -2.72705, -9.73007}};
    // =============== draw points =============
    vtkSPtrNew( polyData, vtkPolyData );
    vtkSPtrNew( points, vtkPoints );
    vtkSPtrNew( cells, vtkCellArray );
    for( int i = 0; i < samplePts.size(); ++i )
    {
        points->InsertNextPoint( samplePts[i].point );
        vtkIdType ids[] = { i };
        cells->InsertNextCell( 1, ids );
    }
    polyData->SetPoints( points );
    polyData->SetVerts( cells );
    polyData->Modified();

    vtkSPtrNew( ptsMapper, vtkPolyDataMapper );
    ptsMapper->SetInputData( polyData );

    vtkSPtrNew( ptsActor, vtkActor );
    ptsActor->SetMapper( ptsMapper );
    ptsActor->GetProperty()->SetPointSize( 2 );
    // ==========================================

    // ------------------------- compute obb -------------------------
    vtkSPtrNew( virtualData, vtkPolyData );
    virtualData->DeepCopy( polyData );
    Point tmpPt( samplePts[0] );
    tmpPt = tmpPt + (normal + (samplePts[1] - samplePts[0]).Unit() )*7;
    virtualData->GetPoints()->InsertNextPoint( tmpPt.point ); // insert new point to form 3D object

    tmpPt = Point( samplePts[samplePts.size()-1] );
    tmpPt = tmpPt + (normal + (samplePts[samplePts.size()-1] - samplePts[0]).Unit() )*7;
    virtualData->GetPoints()->InsertNextPoint( tmpPt.point );

    vtkSmartPointer<vtkDelaunay3D> delaunay = vtkSmartPointer<vtkDelaunay3D>::New();
    delaunay->SetInputData( virtualData );
    delaunay->Update();

    vtkSmartPointer<vtkDataSetSurfaceFilter> surfaceFilter = vtkSmartPointer<vtkDataSetSurfaceFilter>::New();
    surfaceFilter->SetInputConnection(delaunay->GetOutputPort());
    surfaceFilter->Update();

    vtkSPtrNew( obbTree, vtkOBBTree );
    Point corner, sizeVec[3], size;
    obbTree->ComputeOBB( surfaceFilter->GetOutput(), corner.point, sizeVec[0].point, sizeVec[1].point, sizeVec[2].point, size.point );
    std::sort( sizeVec, sizeVec+3, [](Point vec1, Point vec2){ return vec1.Length() < vec2.Length(); } );

    // -------------- make vectors orthogonal ------------------
    Point crossVec = normal^sizeVec[2];
    crossVec.Unit();
    if( crossVec.Dot( sizeVec[1] ) < 0 )
    {
        crossVec = -crossVec;
    }
    sizeVec[1] = crossVec*sizeVec[1].Length();

    // --------------- adjust center --------------------
    auto center = GetCenterOfPoints( samplePts, sizeVec[2], sizeVec[1] );


    Point p1 = center + sizeVec[2]*0.5 - sizeVec[1]*0.5;
    Point p2 = center - sizeVec[2]*0.5 + sizeVec[1]*0.5;
    Point p0 = center - sizeVec[2]*0.5 - sizeVec[1]*0.5;

    vtkSmartPointer<vtkPlaneSource> planeSource =
            vtkSmartPointer<vtkPlaneSource>::New();
    planeSource->SetOrigin( p0.point );
    planeSource->SetPoint1( p1.point );
    planeSource->SetPoint2( p2.point );
    planeSource->Update();

    vtkSmartPointer<vtkPolyDataMapper> planeMapper =
            vtkSmartPointer<vtkPolyDataMapper>::New();
    planeMapper->SetInputData( planeSource->GetOutput() );

    vtkSmartPointer<vtkActor> planeActor =
            vtkSmartPointer<vtkActor>::New();
    planeActor->SetMapper( planeMapper );
    planeActor->GetProperty()->SetColor( 0, 0.5, 0.5 );
    planeActor->GetProperty()->SetOpacity( 0.5 );

    vtkSPtrNew( renderer, vtkRenderer );
    renderer->AddActor( ptsActor );
    renderer->AddActor( planeActor );
    renderer->SetBackground( 0, 0, 0 );

    vtkSPtrNew( renderWindow, vtkRenderWindow );
    renderWindow->AddRenderer( renderer );

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

    vtkSPtrNew(style, vtkInteractorStyleTrackballCamera);
    renderWindowInteractor->SetInteractorStyle(style);

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

    return 0;
}