[vtkusers] using stencils to render stl objects

[Richard Beare]

Hi,

I have an stl model generated by solidworks that I'm trying to render
into a 3D image. I've successfully created various views but run into
issues when trying to consolidate the various world coordinate
information, and I experience segmentation faults depending on
transformation of the polydata.

Anyway, my pipeline is

1) Load a template image using itk. This image is derived from the MRI
scan which I'll eventually register the rendered model to

2) transfer the itk image to vtk.

3) Load the stl using the vtkSTLReader

4) Apply a series of transformations, at the moment done individually
so I can figure out what's going on.

5) Generate a stencil and apply it.

My questions are: should the vtkPolyDataToImageStencil have problems
if the stl and image bounds don't overlap? I've attempted to use a
clipping filter to eliminate potential problems, but no luck. And,
should I be doing this a different way? The code is as follows:

template <class PixType, int dim>
void createPhantom(const CmdLineType &CmdLineObj)
{
    vtkSTLReader *stlReader = vtkSTLReader::New();
    vtkTransformPolyDataFilter *stlTransform1 =
vtkTransformPolyDataFilter::New();
    vtkTransformPolyDataFilter *stlTransform2 =
vtkTransformPolyDataFilter::New();
    vtkTransformPolyDataFilter *stlTransform3 =
vtkTransformPolyDataFilter::New();
    vtkTransform * tf1 = vtkTransform::New();
    vtkTransform * tf2 = vtkTransform::New();
    vtkTransform * tf3 = vtkTransform::New();

    // Load up the sample image using itk
    typedef typename itk::Image<PixType, dim> ImageType;
    typename ImageType::Pointer TempIm =
readIm<ImageType>(CmdLineObj.templateIm);

    // figure out world coordinates of image centre
    typedef typename itk::ContinuousIndex<double,
ImageType::ImageDimension> ContinuousIndexType;

    typename ImageType::SizeType size =
TempIm->GetLargestPossibleRegion().GetSize();
    ContinuousIndexType Middle;
    typename ImageType::PointType WorldMiddle, ModelMiddle;
    for (unsigned i = 0; i < ImageType::ImageDimension; i++)
    {
        Middle[i] = ((float)size[i])/2;
    }
    TempIm->TransformContinuousIndexToPhysicalPoint(Middle, WorldMiddle);
    TempIm->FillBuffer(100);
    // import the image to vtk
    typedef typename itk::ImageToVTKImageFilter<ImageType> itk2vtkType;
    typename itk2vtkType::Pointer itkexporter = itk2vtkType::New();
    itkexporter->SetInput(TempIm);
    itkexporter->Update();

    typename ImageType::PointType origin = TempIm->GetOrigin();
    std::cout << origin << std::endl;
    stlReader->SetFileName(CmdLineObj.StlFile.c_str());
    stlReader->Update();
    std::cout << "Loaded" << std::endl;

    printBounds(stlReader->GetOutput());
    // apply the transforms one at a time - inefficient but clear
    // first centre the cad model
    computeMiddle(stlReader->GetOutput(), ModelMiddle);
    tf1->Translate(-ModelMiddle[0], -ModelMiddle[1], -ModelMiddle[2]);

    std::cout << ModelMiddle;
    stlTransform1->SetInput(stlReader->GetOutput());
    stlTransform1->SetTransform(tf1);
    stlTransform1->Update();

    printBounds(stlTransform1->GetOutput());

    tf2->RotateX(90);
    stlTransform2->SetTransform(tf2);
    stlTransform2->SetInput(stlTransform1->GetOutput());
    stlTransform2->Update();

    computeMiddle(stlTransform2->GetOutput(), ModelMiddle);

    // we want the centre of the phantom to be at the image centre

    stlTransform3->SetInput(stlTransform2->GetOutput());
    tf3->Translate(WorldMiddle[0] - ModelMiddle[0],
                   WorldMiddle[1] - ModelMiddle[1],
                   WorldMiddle[2] - ModelMiddle[2]);

    stlTransform3->SetTransform(tf3);
    stlTransform3->Update();

    std::cout << "Bounds of stl object are: ";
//    double bounds[6];
//    stlReader->GetOutput()->GetBounds(bounds);
    double bounds[6];
    stlTransform3->GetOutput()->ComputeBounds();
    stlTransform3->GetOutput()->GetBounds(bounds);
    {
    for (unsigned i = 0; i<6; i++)
      {
      std::cout << bounds[i] << ",";
      }
    std::cout << std::endl << "Make sure the template image is
appropriate" << std::endl;
    }

    // clip the polydata against the image bounds
    double imagebounds[6];
    vtkBox *imagebox = vtkBox::New();
    vtkClipPolyData *boxclipper = vtkClipPolyData::New();

    itkexporter->GetOutput()->GetBounds(imagebounds);
    imagebox->SetBounds(imagebounds);
    boxclipper->SetClipFunction(imagebox);
    boxclipper->SetInput(stlTransform3->GetOutput());


    vtkPolyDataToImageStencil *surfaceConverter =
vtkPolyDataToImageStencil::New();
//    surfaceConverter->SetInput(stlReader->GetOutput());
    surfaceConverter->SetInput(boxclipper->GetOutput());
    surfaceConverter->SetInformationInput(itkexporter->GetOutput());
    surfaceConverter->SetTolerance( 0.0 );

    // set up the stencil
    vtkImageStencil * stencil = vtkImageStencil::New();
    stencil->SetInput( itkexporter->GetOutput() );
    stencil->SetStencil( surfaceConverter->GetOutput() );
    stencil->SetBackgroundValue( 200 );
    stencil->ReverseStencilOn();
    stencil->Update();

    // send the result back to itk for saving
    typedef typename itk::VTKImageToImageFilter<ImageType> vtk2itkType;
    typename vtk2itkType::Pointer itkimporter = vtk2itkType::New();
    itkimporter->SetInput(stencil->GetOutput());
    itkimporter->Update();

    typedef typename itk::MaskImageFilter<ImageType, ImageType,
ImageType> MaskType;
    typename MaskType::Pointer maskfilt = MaskType::New();
    maskfilt->SetInput1(itkimporter->GetOutput());
    maskfilt->SetInput2(TempIm);

    typedef typename itk::ImageFileWriter<ImageType> WriterType;
    typename WriterType::Pointer writer = WriterType::New();
    writer->SetInput(itkimporter->GetOutput());
    writer->SetFileName(CmdLineObj.OutputIm.c_str());
    writer->Update();

}