[vtkusers] using stencils to render stl objects

[David Gobbi]

Hi Richard,

I'll see if I can answer some of your questions (some of them are
easier than others).  Note that, if you would be satisfied with a pure
ITK solution, there is an itkTriangleMeshToBinaryImageFilter class.
It doesn't use the same algorithm as PolyDataToImageStencil but it is
meant to do essentially the same thing.

Notes:

1) make sure that you are using VTK 5.2 or VTK cvs, since
PolyDataToImageStencil had some bugs in VTK 5.0

2) you don't have to clip the polydata, in fact it is best if you don't

3) try shrinking the polydata with a transform so that its bounds fit
completely within the image bounds, just as a test (also try replacing
the STL polydata with a vtkSphereSource to see if that works).

Of course, things are not supposed to crash... so is there any chance
that you could send me the STL file and the spacing/origin/size of the
image?  If you can, then I'll email you information about how to
connect to my ftp.

   David




On Tue, Aug 5, 2008 at 6:06 AM, Richard Beare <richard.beare at gmail.com> wrote:
> 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();
>
> }
> _______________________________________________
> This is the private VTK discussion list.
> Please keep messages on-topic. Check the FAQ at: http://www.vtk.org/Wiki/VTK_FAQ
> Follow this link to subscribe/unsubscribe:
> http://www.vtk.org/mailman/listinfo/vtkusers
>