Add Two vtkImageData To Avoid Connectivity Removing Object

We create a cone and a sphere individually, then combine the cone and the sphere, show it in the middle renderer.

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

unsigned char inval = 255;
unsigned char outval = 0;
double spacingVal = 0.5;

vtkSPtr<vtkImageData> ConvertPolydataToImage(vtkPolyData *polydata)
{
    double spacingVal = 0.5;
    double bounds[6];
    polydata->GetBounds(bounds);

    int dim[3];
    for (int i = 0; i < 3; i++)
    {
        dim[i] = static_cast<int>( ceil((bounds[2 * i + 1] - bounds[2 * i]) / spacingVal) );
    }
    double origin[3];
    origin[0] = bounds[0] + spacingVal / 2;
    origin[1] = bounds[2] + spacingVal / 2;
    origin[2] = bounds[4] + spacingVal / 2;

    vtkSPtrNew(imageData, vtkImageData)
    imageData->SetSpacing(spacingVal, spacingVal, spacingVal);
    imageData->SetDimensions(dim);
    imageData->SetOrigin(origin);
    imageData->AllocateScalars(VTK_UNSIGNED_CHAR, 1);
    imageData->SetExtent( 0, dim[0] - 1, 0, dim[1] - 1, 0, dim[2] - 1 );

    // fill the imageData with foreground voxels
    vtkIdType count = imageData->GetNumberOfPoints();
    for (vtkIdType i = 0; i < count; ++i)
    {
        imageData->GetPointData()->GetScalars()->SetTuple1(i, inval);
    }

    // polygonal data --> imageData stencil:
    vtkSPtrNew(pd2stenc, vtkPolyDataToImageStencil);
    pd2stenc->SetInputData(polydata);
    pd2stenc->SetOutputOrigin(origin);
    pd2stenc->SetOutputSpacing(spacingVal, spacingVal, spacingVal);
    pd2stenc->SetOutputWholeExtent(imageData->GetExtent());
    pd2stenc->Update();

    // cut the corresponding white imageData and set the background:
    vtkSPtrNew(imgstenc, vtkImageStencil);
    imgstenc->SetInputData(imageData);
    imgstenc->SetStencilConnection(pd2stenc->GetOutputPort());
    imgstenc->ReverseStencilOff();
    imgstenc->SetBackgroundValue( outval );
    imgstenc->Update();

    imageData->DeepCopy(imgstenc->GetOutput());
    return imageData;
}

int main(int, char *[])
{
    setbuf( stdout, nullptr );

    vtkSPtrNew( sphereSource, vtkSphereSource );
    sphereSource->SetRadius(20);
    sphereSource->SetPhiResolution(30);
    sphereSource->SetThetaResolution(30);
    sphereSource->Update();
    vtkSPtr<vtkPolyData> spherePd = sphereSource->GetOutput();
    vtkSPtrNew( sphereMapper, vtkPolyDataMapper );
    sphereMapper->SetInputData( spherePd );
    sphereMapper->Update();
    vtkSPtrNew( sphereActor, vtkActor );
    sphereActor->SetMapper( sphereMapper );

    vtkSPtrNew( coneSource, vtkConeSource );
    coneSource->SetRadius( 20 );
    coneSource->SetHeight( 50 );
    coneSource->SetResolution( 30 );
    coneSource->Update();
    vtkSPtr<vtkPolyData> conePd = coneSource->GetOutput();
    vtkSPtrNew( coneMapper, vtkPolyDataMapper );
    coneMapper->SetInputData( conePd );
    coneMapper->Update();
    vtkSPtrNew( coneActor, vtkActor );
    coneActor->SetMapper( coneMapper );

    vtkSPtrNew( appendFilter, vtkAppendPolyData );
    appendFilter->AddInputData( conePd );
    appendFilter->AddInputData( spherePd );
    appendFilter->Update();
    vtkSPtrNew( appendMapper, vtkPolyDataMapper );
    appendMapper->SetInputData( appendFilter->GetOutput() );
    appendMapper->Update();
    vtkSPtrNew( appendActor, vtkActor );
    appendActor->SetMapper( appendMapper );

    double leftViewport[4] = {0.0, 0.0, 0.3, 1.0};
    double midViewport[4] = {0.3, 0.0, 0.7, 1.0};
    double rightViewport[4] = {0.7, 0.0, 1.0, 1.0};

    // Setup renderers
    vtkSPtrNew( leftRenderer, vtkRenderer );
    leftRenderer->SetViewport( leftViewport );
    leftRenderer->AddActor( sphereActor );
    leftRenderer->SetBackground(.6, .5, .4);
    leftRenderer->ResetCamera();

    vtkSPtrNew( midRenderer, vtkRenderer );
    midRenderer->SetViewport(midViewport);
    midRenderer->AddActor( appendActor );
    midRenderer->ResetCamera();

    vtkSPtrNew( rightRenderer, vtkRenderer );
    rightRenderer->SetViewport(rightViewport);
    rightRenderer->AddActor( coneActor );
    rightRenderer->SetBackground(.4, .5, .6);
    rightRenderer->ResetCamera();

    // Setup render window
    vtkSPtrNew( renderWindow, vtkRenderWindow );
    renderWindow->AddRenderer( leftRenderer );
    renderWindow->AddRenderer( midRenderer );
    renderWindow->AddRenderer( rightRenderer );
    renderWindow->SetSize( 900, 300 );

    // Setup render window interactor
    vtkSPtrNew( renderWindowInteractor, vtkRenderWindowInteractor );

    // Render and start interaction
    renderWindowInteractor->SetRenderWindow(renderWindow);
    renderWindowInteractor->Start();

    return EXIT_SUCCESS;
}

We use vtkPolyDataConnectivityFilter to fetch the biggest part of the combination, the result is only sphere.

    vtkSPtrNew( connectivity, vtkPolyDataConnectivityFilter );
    connectivity->SetInputData( appendFilter->GetOutput() );
    connectivity->Update();
    vtkSPtrNew( appendMapper, vtkPolyDataMapper );
    appendMapper->SetInputData( connectivity->GetOutput() );
    appendMapper->Update();

We can create a big image data which contains cone and sphere, vtkPolyDataConnectivityFilter will regards it as a complete object.

void AddImageData( vtkSPtr<vtkImageData> finalImage, vtkSPtr<vtkImageData> image )
{
    int *dim = image->GetDimensions();
    for( int i = 0; i < dim[0]; ++i )
    {
        for( int j = 0; j < dim[1]; ++j )
        {
            for( int k = 0; k < dim[2]; ++k )
            {
                char *pixel = static_cast<char *>(image->GetScalarPointer(i,j,k));
                if( pixel[0] == outval )
                {
                    continue;
                }
                int index[3] = {i, j, k};
                int pointId = image->ComputePointId( index );
                double pos[3];
                image->GetPoint( pointId, pos );
                int finalImagePointId = finalImage->FindPoint( pos );
                if( -1 != finalImagePointId )
                {
                    finalImage->GetPointData()->GetScalars()->SetTuple1( finalImagePointId, inval );
                }
            }
        }
    }
}

vtkSPtr<vtkImageData> GetFinalImage(std::vector<vtkSPtr<vtkImageData>> images, double spacingVal)
{
    double bounds[6] = { VTK_DOUBLE_MAX, VTK_DOUBLE_MIN,
                         VTK_DOUBLE_MAX, VTK_DOUBLE_MIN,
                         VTK_DOUBLE_MAX, VTK_DOUBLE_MIN };
    for( int i = 0; i < images.size(); ++i )
    {
        double origin[3];
        int dim[3];
        images[i]->GetOrigin( origin );
        images[i]->GetDimensions( dim );
        for( int j = 0; j < 3; ++j )
        {
            if( bounds[2*j] > origin[j] )
            {
                bounds[2*j] = origin[j];
            }
        }
        for( int j = 0; j < 3; ++j )
        {
            double boundsValue = origin[j] + dim[j] * spacingVal;
            if( boundsValue > bounds[2*j + 1] )
            {
                bounds[2*j + 1] = boundsValue;
            }
        }
    }
    int dim[3];
    for( int i = 0; i < 3; ++i )
    {
        dim[i] = static_cast<int>( ceil((bounds[2*i + 1] - bounds[2*i]) / spacingVal) );
    }
    double origin[3];
    origin[0] = bounds[0] + spacingVal / 2;
    origin[1] = bounds[2] + spacingVal / 2;
    origin[2] = bounds[4] + spacingVal / 2;

    vtkSPtrNew( finalImage, vtkImageData );
    finalImage->SetSpacing( spacingVal, spacingVal, spacingVal );
    finalImage->SetDimensions( dim );
    finalImage->SetOrigin( origin );
    finalImage->AllocateScalars( VTK_UNSIGNED_CHAR, 1 );

    vtkIdType count = finalImage->GetNumberOfPoints();
    for( vtkIdType i = 0; i < count; ++i )
    {
        finalImage->GetPointData()->GetScalars()->SetTuple1( i, outval );
    }
    return finalImage;
}

int main()
{
    // ================= add image ==================
    vtkSPtr<vtkImageData> sphereImage = ConvertPolydataToImage( spherePd );
    vtkSPtr<vtkImageData> coneImage = ConvertPolydataToImage( conePd );
    std::vector<vtkSPtr<vtkImageData>> images;
    images.push_back( sphereImage );
    images.push_back( coneImage );
    vtkSPtr<vtkImageData> appendImage = GetFinalImage( images, spacingVal );
    for( auto image : images )
    {
        AddImageData( appendImage, image );
    }
    // ================= finished: add image ==================

    // ============== convert image ===============
    vtkSPtrNew( marchingCubes, vtkImageMarchingCubes );
    marchingCubes->SetInputData( appendImage );
    marchingCubes->SetValue( 0, (inval + outval)/2 );
    marchingCubes->SetNumberOfContours( 1 );
    marchingCubes->Update();

    // Fill the holes
    auto fillHoles = vtkSmartPointer<vtkFillHolesFilter>::New();
    fillHoles->SetInputConnection( marchingCubes->GetOutputPort() );
    fillHoles->SetHoleSize(1000.0);

    // Make the triangle winding order consistent
    auto normals = vtkSmartPointer<vtkPolyDataNormals>::New();
    normals->SetInputConnection(fillHoles->GetOutputPort());
    normals->ConsistencyOn();
    normals->SplittingOff();
    normals->GetOutput()->GetPointData()->SetNormals(
        marchingCubes->GetOutput()->GetPointData()->GetNormals());
    normals->Update();
    // ==============finished: convert image ===============

    vtkSPtrNew( connectivity, vtkPolyDataConnectivityFilter );
    connectivity->SetInputData( normals->GetOutput() );
    connectivity->Update();

    vtkSPtrNew( appendMapper, vtkPolyDataMapper );
    appendMapper->SetInputData( connectivity->GetOutput() );
    appendMapper->Update();

    vtkSPtrNew( appendActor, vtkActor );
    appendActor->SetMapper( appendMapper );
}

We can also subtract an image to form a new object.

void SubImageData( vtkSPtr<vtkImageData> finalImage, vtkSPtr<vtkImageData> image )
{
    int *dim = image->GetDimensions();
    for( int i = 0; i < dim[0]; ++i )
    {
        for( int j = 0; j < dim[1]; ++j )
        {
            for( int k = 0; k < dim[2]; ++k )
            {
                char *pixel = static_cast<char *>(image->GetScalarPointer(i,j,k));
                if( pixel[0] == outval )
                {
                    continue;
                }
                int index[3] = {i, j, k};
                int pointId = image->ComputePointId( index );
                double pos[3];
                image->GetPoint( pointId, pos );
                int finalImagePointId = finalImage->FindPoint( pos );
                if( -1 != finalImagePointId )
                {
                    finalImage->GetPointData()->GetScalars()->SetTuple1( finalImagePointId, outval );
                }
            }
        }
    }
}

The following code snippet contains all header files in our cpp program.

#include <vtkVersion.h>
#include <vtkImageData.h>
#include <vtkImageMapper3D.h>
#include <vtkImageStencil.h>
#include <vtkImageStencilData.h>
#include <vtkImageToImageStencil.h>
#include <vtkJPEGReader.h>
#include <vtkPointData.h>
#include <vtkSmartPointer.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkInteractorStyleImage.h>
#include <vtkRenderer.h>
#include <vtkImageActor.h>
#include <vtkPolyDataToImageStencil.h>
#include <vtkPolyData.h>
#include <vtkSTLReader.h>
#include <vtkProperty.h>
#include <vtkImageProperty.h>
#include <vtkSphereSource.h>
#include <vtkPolyDataMapper.h>
#include <vtkImageDataGeometryFilter.h>
#include <vtkImageMarchingCubes.h>
#include <vtkPolyDataNormals.h>
#include <vtkFillHolesFilter.h>
#include <vtkCleanPolyData.h>
#include <vtkTriangleFilter.h>
#include <vtkAppendPolyData.h>
#include <vtkImageMarchingCubes.h>
#include <vtkImageWeightedSum.h>
#include <vtkConeSource.h>
#include <vtkPolyDataConnectivityFilter.h>