I want to change the camera’s position and focal point to focus on a small region on the surface of the sphere. The parameters viewUp direction, focal point, position, and parallel scale are important for our target.

Our scene comes from a previous article VTK – Find Partial Region And Its Normal, let’s create a sphere and find a red region firstly.

Here is a difference when creating a red region in our current exploring.

    // =============== find a red area ===============
    vtkSPtrNew( visitedCellIds, vtkIdList );
    int firstCellId = pd->GetNumberOfPoints() - 1;
    visitedCellIds->InsertNextId( firstCellId );
    std::vector<vtkIdType> layerCellIds;
    layerCellIds.push_back( firstCellId );
    int iterCount = 3;

After getting the normal and center of the red area, we can adjust the camera. The parallel scale is related to the size of the red area. The viewUp direction can be created by the focal plane. All implementation details are showed in the following code snippet.

void ComputeBoundsForPtsCloud(double *bounds, vtkSmartPointer<vtkPoints> points)
{
    bounds[0] = bounds[2] = bounds[4] = VTK_DOUBLE_MAX;
    bounds[1] = bounds[3] = bounds[5] = VTK_DOUBLE_MIN;
    for( int i = 0; i < points->GetNumberOfPoints(); ++i )
    {
        double *pt = points->GetPoint( i );
        if( bounds[0] > pt[0] )
        {
            bounds[0] = pt[0];
        }
        if( bounds[1] < pt[0] )
        {
            bounds[1] = pt[0];
        }


        if( bounds[2] > pt[1] )
        {
            bounds[2] = pt[1];
        }
        if( bounds[3] < pt[1] )
        {
            bounds[3] = pt[1];
        }


        if( bounds[4] > pt[2] )
        {
            bounds[4] = pt[2];
        }
        if( bounds[5] < pt[2] )
        {
            bounds[5] = pt[2];
        }
    }
}

/*
*   viewProject is start from focalPt to camera position
*/
PointStruct ComputeViewUp( PointStruct focalPt, PointStruct viewProject )
{
    vtkSPtrNew( focalPlane, vtkPlane );
    focalPlane->SetOrigin( focalPt.point );
    focalPlane->SetNormal( viewProject.point );

    PointStruct ptOnViewLine = focalPt + viewProject;
    PointStruct tmpPt = ptOnViewLine + PointStruct(1, 0, 0);
    // tmpPt and ptOnViewLine are on the view project line.
    if( fabs( (tmpPt - ptOnViewLine).Dot( viewProject ) - 1 ) < 1e-6 )
    {
        tmpPt = ptOnViewLine + PointStruct(0, 1, 0);
    }
    focalPlane->ProjectPoint( tmpPt.point, tmpPt.point );

    PointStruct result = tmpPt - focalPt;
    result.Unit();
    return result;
}

void AdjustCameraToFocus( PointStruct origin, PointStruct normal,
                          vtkSPtr<vtkRenderer> renderer, vtkSPtr<vtkPoints> points )
{
    vtkCamera *camera = renderer->GetActiveCamera();
    PointStruct curCamPos( camera->GetPosition() );
    PointStruct curFocalPt( camera->GetFocalPoint() );
    PointStruct viewDirection;
    vtkMath::Subtract(curCamPos.point, curFocalPt.point, viewDirection.point);

    double bds[6];
    ComputeBoundsForPtsCloud( bds, points );

    double w1 = bds[1] - bds[0];
    double w2 = bds[3] - bds[2];
    double w3 = bds[5] - bds[4];
    w1 *= w1;
    w2 *= w2;
    w3 *= w3;
    double radius = w1 + w2 + w3;

    // If we have just a single point, pick a radius of 1.0
    radius = (fabs(radius) < 1e-6 ) ? (1.0) : (radius);

    // compute the radius of the enclosing sphere
    radius = sqrt( radius )*0.5;
    double parallelScale = radius;

    double aspect[2];
    renderer->GetAspect(aspect);
    std::cout << "aspect[0]: " << aspect[0] << std::endl;
    if( aspect[0] < 1.0 )
    {
        parallelScale = parallelScale/aspect[0];
    }

    camera->SetFocalPoint( origin.point );
    PointStruct pos = origin + normal*viewDirection.Length();
    camera->SetPosition( pos.point );
    PointStruct viewUpDir = ComputeViewUp( origin, normal );
    camera->SetViewUp( viewUpDir.point );
    std::cout << "parallelScale: " << parallelScale << std::endl;
    camera->SetParallelScale( parallelScale );
}

int main()
{
//...

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

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

    renderer->GetActiveCamera()->SetParallelProjection( 1 );
    renderer->ResetCamera();

    // =============== start to adjust camera ==============
    vtkSPtrNew( bdsPoints, vtkPoints );
    for( int i = 0; i < boundaryPtIds->GetNumberOfIds(); ++i )
    {
        PointStruct pt( boundaryPd->GetPoint( boundaryPtIds->GetId( i ) ) );
        bdsPoints->InsertNextPoint( pt.point );
    }
    AdjustCameraToFocus( origin, normal, renderer, bdsPoints );
//...

The complete code of the program.

#include <iostream>
#include <vtkSmartPointer.h>
#include <vtkSphereSource.h>
#include <vtkActor.h>
#include <vtkConeSource.h>
#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
#include <vtkPolyDataMapper.h>
#include <vtkLineSource.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkAxesActor.h>
#include <vtkProperty.h>
#include <vtkXMLPolyDataReader.h>
#include <vtkCharArray.h>
#include <vtkPointData.h>
#include <vtkConeSource.h>
#include <vtkTransform.h>
#include <vtkColorTransferFunction.h>
#include <vtkCellData.h>
#include <vtkPoints.h>
#include <vtkCell.h>
#include <vtkCellArray.h>
#include <vtkTriangle.h>
#include <vtkLineSource.h>
#include <vtkTriangleFilter.h>
#include <vtkCamera.h>
#include <vtkPlane.h>

#include "../share/tool.h"

using namespace std;

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

const int BackCell = 0;
const int ChosenCell = 1;

void FindIndependentEdges( vtkSPtr<vtkCellArray> edges, vtkSPtr<vtkIdList> redCells, vtkSPtr<vtkPolyData> originPolyData )
{
    vtkIdType cellId;
    vtkIdType npts;
    vtkIdType p1, p2;
    vtkIdType *pts = nullptr;
    edges->Initialize();
    if( originPolyData->NeedToBuildCells() )
    {
        originPolyData->BuildLinks();
        originPolyData->BuildCells();
    }

    vtkSPtrNew( neighbors, vtkIdList );
    neighbors->Allocate(VTK_CELL_SIZE);
    vtkDataArray *cellScalars = originPolyData->GetCellData()->GetScalars();
    for( int index = 0; index < redCells->GetNumberOfIds(); ++index )
    {
        cellId = redCells->GetId( index );
        originPolyData->GetCellPoints( cellId, npts, pts );

        for (int i = 0; i < npts; i++ )
        {
            p1 = pts[ i ];
            p2 = pts[ (i+1) %npts ];

            originPolyData->GetCellEdgeNeighbors( cellId, p1, p2, neighbors );
            bool noNeighbor = true;
            for( int j = 0; j < neighbors->GetNumberOfIds(); ++j )
            {
                vtkIdType neiCellId = neighbors->GetId( j );
                double value = cellScalars->GetTuple1( neiCellId );
                if( fabs(value - ChosenCell) < 1e-6 )
                {
                    noNeighbor = false;
                    break;
                }
            }
            if( noNeighbor )
            {
                vtkIdType edgeIds[2] = { p1, p2 };
                edges->InsertNextCell( 2, edgeIds );
            }
        }
    }
}

PointStruct CalculateNormal(vtkSPtr<vtkIdList> redCells, vtkSPtr<vtkPolyData> originPolyData)
{
    PointStruct totalNormal(0, 0, 0);
    for (vtkIdType i = 0; i < redCells->GetNumberOfIds(); i++)
    {
        vtkIdType cellId = redCells->GetId(i);
        vtkIdType nPoints;
        vtkIdType* pts;
        originPolyData->GetCellPoints(cellId, nPoints, pts);
        if (nPoints != 3)
            continue;

        PointStruct pointPos[3];
        for (vtkIdType j = 0; j < nPoints; j++)
        {
            originPolyData->GetPoint(pts[j], pointPos[j].point);
        }
        PointStruct vec1 = pointPos[1] - pointPos[0];
        PointStruct vec2 = pointPos[2] - pointPos[1];
        PointStruct normal = vec1 ^ vec2;
        double area = vtkTriangle::TriangleArea(pointPos[0].point, pointPos[1].point, pointPos[2].point);
        normal.Unit();
        totalNormal += normal*area;
    }
    if( totalNormal.Length() > 0 )
    {
        totalNormal.Unit();
    }
    return totalNormal;
}

void ComputeBoundsForPtsCloud(double *bounds, vtkSmartPointer<vtkPoints> points)
{
    bounds[0] = bounds[2] = bounds[4] = VTK_DOUBLE_MAX;
    bounds[1] = bounds[3] = bounds[5] = VTK_DOUBLE_MIN;
    for( int i = 0; i < points->GetNumberOfPoints(); ++i )
    {
        double *pt = points->GetPoint( i );
        if( bounds[0] > pt[0] )
        {
            bounds[0] = pt[0];
        }
        if( bounds[1] < pt[0] )
        {
            bounds[1] = pt[0];
        }


        if( bounds[2] > pt[1] )
        {
            bounds[2] = pt[1];
        }
        if( bounds[3] < pt[1] )
        {
            bounds[3] = pt[1];
        }


        if( bounds[4] > pt[2] )
        {
            bounds[4] = pt[2];
        }
        if( bounds[5] < pt[2] )
        {
            bounds[5] = pt[2];
        }
    }
}

/*
*   viewProject is start from focalPt to camera position
*/
PointStruct ComputeViewUp( PointStruct focalPt, PointStruct viewProject )
{
    vtkSPtrNew( focalPlane, vtkPlane );
    focalPlane->SetOrigin( focalPt.point );
    focalPlane->SetNormal( viewProject.point );

    PointStruct ptOnViewLine = focalPt + viewProject;
    PointStruct tmpPt = ptOnViewLine + PointStruct(1, 0, 0);
    // tmpPt and ptOnViewLine are on the view project line.
    if( fabs( (tmpPt - ptOnViewLine).Dot( viewProject ) - 1 ) < 1e-6 )
    {
        tmpPt = ptOnViewLine + PointStruct(0, 1, 0);
    }
    focalPlane->ProjectPoint( tmpPt.point, tmpPt.point );

    PointStruct result = tmpPt - focalPt;
    result.Unit();
    return result;
}

void AdjustCameraToFocus( PointStruct origin, PointStruct normal,
                          vtkSPtr<vtkRenderer> renderer, vtkSPtr<vtkPoints> points )
{
    vtkCamera *camera = renderer->GetActiveCamera();
    PointStruct curCamPos( camera->GetPosition() );
    PointStruct curFocalPt( camera->GetFocalPoint() );
    PointStruct viewDirection;
    vtkMath::Subtract(curCamPos.point, curFocalPt.point, viewDirection.point);

    double bds[6];
    ComputeBoundsForPtsCloud( bds, points );

    double w1 = bds[1] - bds[0];
    double w2 = bds[3] - bds[2];
    double w3 = bds[5] - bds[4];
    w1 *= w1;
    w2 *= w2;
    w3 *= w3;
    double radius = w1 + w2 + w3;

    // If we have just a single point, pick a radius of 1.0
    radius = (fabs(radius) < 1e-6 ) ? (1.0) : (radius);

    // compute the radius of the enclosing sphere
    radius = sqrt( radius )*0.5;
    double parallelScale = radius;

    double aspect[2];
    renderer->GetAspect(aspect);
    std::cout << "aspect[0]: " << aspect[0] << std::endl;
    if( aspect[0] < 1.0 )
    {
        parallelScale = parallelScale/aspect[0];
    }

    camera->SetFocalPoint( origin.point );
    PointStruct pos = origin + normal*viewDirection.Length();
    camera->SetPosition( pos.point );
    PointStruct viewUpDir = ComputeViewUp( origin, normal );
    camera->SetViewUp( viewUpDir.point );
    std::cout << "parallelScale: " << parallelScale << std::endl;
    camera->SetParallelScale( parallelScale );
}

int main()
{
    vtkSPtrNew( sphereSource, vtkSphereSource );
    sphereSource->SetThetaResolution( 50 );
    sphereSource->SetPhiResolution( 50 );
    sphereSource->SetRadius( 10 );
    sphereSource->Update();

    auto *pd = sphereSource->GetOutput();
    pd->BuildCells();
    pd->BuildLinks();

    // =============== find a red area ===============
    vtkSPtrNew( visitedCellIds, vtkIdList );
    int firstCellId = pd->GetNumberOfPoints() - 1;
    visitedCellIds->InsertNextId( firstCellId );
    std::vector<vtkIdType> layerCellIds;
    layerCellIds.push_back( firstCellId );
    int iterCount = 3;
    while( iterCount-- )
    {
        std::vector<vtkIdType> currentCellIds;
        for( int i = 0; i < layerCellIds.size(); ++i )
        {
            vtkIdType cellId = layerCellIds[i];
            vtkIdType nPoints;
            vtkIdType *pts;
            pd->GetCellPoints(cellId, nPoints, pts);

            for( int j = 0; j < nPoints; ++j )
            {
                vtkSPtrNew( neiCellIds, vtkIdList );
                pd->GetCellEdgeNeighbors( cellId, pts[j], pts[(j+1)%nPoints], neiCellIds );

                for( int k = 0; k < neiCellIds->GetNumberOfIds(); ++k )
                {
                    cellId = neiCellIds->GetId( k );
                    if( -1 == visitedCellIds->IsId( cellId ) )
                    {
                        visitedCellIds->InsertUniqueId( cellId );
                        currentCellIds.push_back( cellId );
                    }
                }
            }
        }
        layerCellIds.swap( currentCellIds );
    }
    cout << "visitedCellIds count: " << visitedCellIds->GetNumberOfIds() << std::endl;

    // ================= use cell scalar to draw it ===================
    vtkSPtrNew( cellTypes, vtkCharArray );
    cellTypes->SetNumberOfComponents( 1 );
    for( int i = 0; i < pd->GetNumberOfCells(); ++i )
    {
        cellTypes->InsertNextValue( BackCell );
    }
    for( int i = 0; i < visitedCellIds->GetNumberOfIds(); ++i )
    {
        auto cellId = visitedCellIds->GetId( i );
        cellTypes->InsertValue( cellId, ChosenCell );
    }

    pd->GetCellData()->SetScalars( cellTypes );
    pd->GetCellData()->Modified();

    vtkSPtrNew( sphereMapper, vtkPolyDataMapper );
    sphereMapper->SetInputData( pd );
    sphereMapper->SetScalarModeToUseCellData();
    sphereMapper->SetColorModeToMapScalars();

    vtkSPtrNew( lut, vtkColorTransferFunction );
    lut->AddRGBPoint( BackCell, 1, 1, 1 );
    lut->AddRGBPoint( ChosenCell, 1, 0, 0 );
    sphereMapper->SetLookupTable( lut );
    // ================= finish: draw =====================

    vtkSPtrNew( actor, vtkActor );
    actor->SetMapper( sphereMapper );

    // ================= show boundary ===================
    vtkSPtrNew( edges, vtkCellArray );
    FindIndependentEdges( edges, visitedCellIds, pd );
    vtkSPtrNew( boundaryPd, vtkPolyData );
    boundaryPd->SetPoints( pd->GetPoints() );
    boundaryPd->SetLines( edges );
    //boundaryPd->SetPolys( edges );
    vtkSPtrNew( boundaryMapper, vtkPolyDataMapper );
    boundaryMapper->SetInputData( boundaryPd );
    vtkSPtrNew( boundrayActor, vtkActor );
    boundrayActor->SetMapper( boundaryMapper );
    boundrayActor->GetProperty()->SetColor( 0, 0, 1 );
    boundrayActor->GetProperty()->SetLineWidth( 5 );

    // ================= show normal ===================
    PointStruct normal = CalculateNormal( visitedCellIds, pd );
    vtkSPtrNew( boundaryPtIds, vtkIdList );
    boundaryPd->BuildCells();
    boundaryPd->BuildLinks();
    for( int i = 0; i < boundaryPd->GetNumberOfCells(); ++i )
    {
        vtkIdType nPts;
        vtkIdType *pts;
        boundaryPd->GetCellPoints( i, nPts, pts );
        for( int j = 0; j < nPts; ++j )
        {
            if( -1 == boundaryPtIds->IsId( pts[j] ) )
            {
                boundaryPtIds->InsertUniqueId( pts[j] );
            }
        }
    }
    PointStruct origin( 0, 0, 0 );

    for( int i = 0; i < boundaryPtIds->GetNumberOfIds(); ++i )
    {
        PointStruct pt( boundaryPd->GetPoint( boundaryPtIds->GetId( i ) ) );
        origin = origin + pt;
    }
    origin /= boundaryPtIds->GetNumberOfIds();

    vtkSPtrNew( line, vtkLineSource );
    line->SetPoint1( origin.point );
    line->SetPoint2( (origin + normal).point );
    line->Update();
    vtkSPtrNew( lineMapper, vtkPolyDataMapper );
    lineMapper->SetInputData( line->GetOutput() );
    vtkSPtrNew( lineActor, vtkActor );
    lineActor->SetMapper( lineMapper );
    lineActor->GetProperty()->SetColor( 0, 1, 0 );
    lineActor->GetProperty()->SetLineWidth( 5 );

    vtkSPtrNew( renderer, vtkRenderer );
    renderer->AddActor( actor );
    renderer->AddActor( boundrayActor );
    renderer->AddActor( lineActor );
    renderer->SetBackground( 0, 0, 0 );

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

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

    renderer->GetActiveCamera()->SetParallelProjection( 1 );
    renderer->ResetCamera();

    // =============== start to adjust camera ==============
    vtkSPtrNew( bdsPoints, vtkPoints );
    for( int i = 0; i < boundaryPtIds->GetNumberOfIds(); ++i )
    {
        PointStruct pt( boundaryPd->GetPoint( boundaryPtIds->GetId( i ) ) );
        bdsPoints->InsertNextPoint( pt.point );
    }
    AdjustCameraToFocus( origin, normal, renderer, bdsPoints );

    renderWindow->Render();
    renderWindowInteractor->Start();
    return 0;
}
Categories: VTK

0 0 vote
Article Rating
Subscribe
Notify of
guest
0 Comments
Inline Feedbacks
View all comments
0
Would love your thoughts, please comment.x
()
x