VTK – Find Partial Region And Its Normal

Published by StephenWei on

Find a local partial region on the surface of the sphere, mark it by red color. Finally, we will compute the normal of the partial region and show it by a green line. The whole scene is just like the following image.

Here are a few steps about the process. Please read them and you can get all the implementation details.

Find A Local Partial Area

  1. Take a cell which id is 0 firstly.
  2. Search all neighbor cells of it and collect them.
  3. Start to find the new neighbor cells of the collected cells in step 2.

Iterate the third step until we find enough cells and store them in a cell id list.

    vtkSPtrNew( visitedCellIds, vtkIdList );
    visitedCellIds->InsertNextId( 0 );
    std::vector<vtkIdType> layerCellIds;
    layerCellIds.push_back( 0 );
    int iterCount = 8;
    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;

Find The Boundary Of The Red Region

Iterate over all of the red cells, collect all the edges which are used by only red cells. Store all edge points in an array, we will use them in the later step.

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

Calculate The Region Normal

Compute every cell’s normal and use them to get an area-weighted normal. It means the bigger cell has more influence on the final result.

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

The complete program code is in the following part.
We can switch to wireframe representation to look at the origin point which is computed from the blue boundary.

#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 "../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;
}

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

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

    // =============== find a red area ===============
    vtkSPtrNew( visitedCellIds, vtkIdList );
    visitedCellIds->InsertNextId( 0 );
    std::vector<vtkIdType> layerCellIds;
    layerCellIds.push_back( 0 );
    int iterCount = 8;
    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 );

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

    renderer->ResetCamera();
    renderWindow->Render();
    renderWindowInteractor->Start();
    return 0;
}
Categories: AlgorithmVTK
Tags:
0 0 votes
Article Rating
Subscribe
Notify of
guest

1 Comment
Oldest
Newest Most Voted
Inline Feedbacks
View all comments
trackback
Adjust Camera To Focus On Small Region | weiy
3 years ago

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

0
Reply

Related Posts

VTK

Size Of VTK Render Window Will Not Match The External Canvas

I found an interesting issue about vtkSDL2OpenGLRenderWindow::SetSize when changed canvas for VTK WASM 3D window.The related post on the website is Change Canvas For VTK WASM Project. The interface SetSize can work if it was used before m_RenderWindowInteractor->SetRenderWindow( newWin );.It means Read more…

VTK

VTK WASM Error – Uncaught TypeError: Cannot read properties of null (reading ‘bindFramebuffer’) at _glBindFramebuffer

There are some errors and warnings after I changed canvas for VTK wasm. The related post: https://www.weiy.city/2023/11/change-canvas-for-vtk-wasm-project/ We need to free up graphics resources with the old window before we render again. worker.cpp:

VTK

Pick Operation Of vtkPropPicker On vtkSDL2OpenGLRenderWindow Is Not Working

Environment:VTK Branch_v9.1.0emcc (Emscripten gcc/clang-like replacement + linker emulating GNU ld) 3.1.39 I tested two interfaces for pick operation of vtkPropPicker on web page. The result value is awalys 0. Invoking code: There are warning on Read more…

Content Summary
 : Input your strings, the tool can get a brief summary of the content for you.

Got it!
X
1
0
Would love your thoughts, please comment.x
()
x