Create Oriented Bounding Box Based On Three Custom Axis

The post shows how to generate an oriented bounding box based on the three axis known in advance.
The 3D model is displayed by a CUModel object. We will calculate the two farthest points, positive farthest point and negative farthest point along the special direction. Then generate a polydata from the six farthest points.
Attention should be paid to the determination of the center point when we create polydata from the six points.

void GetOBBBox(CUModel *model, vtkSmartPointer<vtkPolyData> result)
{
    // we have three directions as axes. PointStruct xDir, yDir, zDir; 

    PointStruct boundPts[6];
    // -x, x, -y, y, -z, z
    CalculateFarPoints( boundPts[1], boundPts[0], xDir, model );
    CalculateFarPoints( boundPts[3], boundPts[2], yDir, model );
    CalculateFarPoints( boundPts[5], boundPts[4], zDir, model );

    // compute the polyData center
    PointStruct center( (boundPts[0]+boundPts[1])*0.5 );
    vSPNew( plane, vtkPlane );
    plane->SetOrigin( center.point );
    plane->SetNormal( xDir.point );
    PointStruct tmp = (boundPts[2]+boundPts[3])*0.5;
    plane->ProjectPoint( tmp.point, center.point ); // x mid, y mid.
    tmp = (boundPts[4]+boundPts[5])*0.5;
    vSPNew( plane1, vtkPlane );
    plane1->SetOrigin( tmp.point );
    plane1->SetNormal( zDir.point );
    plane1->ProjectPoint( center.point, center.point ); // x mid, y mid, z mid.

    // generate polyData
    PointStruct sizeVec[3] = { (boundPts[1] - boundPts[0])*0.5, (boundPts[3] - boundPts[2])*0.5, (boundPts[5] - boundPts[4])*0.5 };
    PointStruct pts[8] = {
        center - sizeVec[0] - sizeVec[1] - sizeVec[2],
        center + sizeVec[0] - sizeVec[1] - sizeVec[2],
        center + sizeVec[0] + sizeVec[1] - sizeVec[2],
        center - sizeVec[0] + sizeVec[1] - sizeVec[2],
        center - sizeVec[0] - sizeVec[1] + sizeVec[2],
        center + sizeVec[0] - sizeVec[1] + sizeVec[2],
        center + sizeVec[0] + sizeVec[1] + sizeVec[2],
        center - sizeVec[0] + sizeVec[1] + sizeVec[2],
    };
    vSPNew( points, vtkPoints );
    for( auto pt: pts )
    {
        points->InsertNextPoint( pt.point );
    }
    vSPNew( cells, vtkCellArray );
    vtkIdType bottom_ptIds1[3] = { 1, 0, 2 };
    cells->InsertNextCell( 3, bottom_ptIds1 );
    vtkIdType bottom_ptIds2[3] = { 3, 2, 0 };
    cells->InsertNextCell( 3, bottom_ptIds2 );

    vtkIdType top_ptIds1[3] = { 4, 5, 6 };
    cells->InsertNextCell( 3, top_ptIds1 );
    vtkIdType top_ptIds2[3] = { 6, 7, 4 };
    cells->InsertNextCell( 3, top_ptIds2 );

    vtkIdType left_ptIds1[3] = { 4, 3, 0 };
    cells->InsertNextCell( 3, left_ptIds1 );
    vtkIdType left_ptIds2[3] = { 3, 4, 7 };
    cells->InsertNextCell( 3, left_ptIds2 );

    vtkIdType right_ptIds1[3] = { 2, 6, 5 };
    cells->InsertNextCell( 3, right_ptIds1 );
    vtkIdType right_ptIds2[3] = { 2, 5, 1 };
    cells->InsertNextCell( 3, right_ptIds2 );

    vtkIdType front_ptIds1[3] = { 1, 4, 0 };
    cells->InsertNextCell( 3, front_ptIds1 );
    vtkIdType front_ptIds2[3] = { 5, 4, 1 };
    cells->InsertNextCell( 3, front_ptIds2 );

    vtkIdType back_ptIds1[3] = { 7, 6, 2 };
    cells->InsertNextCell( 3, back_ptIds1 );
    vtkIdType back_ptIds2[3] = { 7, 2, 3 };
    cells->InsertNextCell( 3, back_ptIds2 );

    result->SetPoints( points );
    result->SetPolys( cells );
    result->BuildCells();
    result->BuildLinks();
    result->Modified();
}

void CalculateFarPoints
    (
    PointStruct &posPt,
    PointStruct &negPt,
    PointStruct axis,
    CUModel *model
    )
{
    auto data = model->Getm_Actor()->GetInputPolyData();
    data->ComputeBounds();
    auto bds = data->GetBounds();
    PointStruct center( (bds[0]+bds[1])/2, (bds[2]+bds[3])/2, (bds[4]+bds[5])/2 );

    vSPNew( computePlane, vtkPlane );
    computePlane->SetOrigin( center.point );
    computePlane->SetNormal( axis.point );

    double maxValue = VTK_DOUBLE_MIN;
    double minValue = VTK_DOUBLE_MAX;
    PointStruct posFarPt(0, 0, 0), negFarPt(0, 0, 0);
    for( int i = 0; i < data->GetNumberOfPoints(); ++i )
    {
        PointStruct pt( data->GetPoint( i ) );
        double computeValue = computePlane->EvaluateFunction( pt.point );
        if( computeValue > maxValue )
        {
            maxValue = computeValue;
            posFarPt = pt;
        }

        if( computeValue < minValue )
        {
            minValue = computeValue;
            negFarPt = pt;
        }
    }

    LOG( INFO, "maxValue: ", maxValue, ", minValue: ", minValue );
    posPt = center + maxValue*axis;
    negPt = center + minValue*axis;
}