[vtkusers] VTKPolyData to 3D image

[David Welch]

Okay, I'm in the home stretch here... I've changed the code slightly to
output a XMLImageData file and take an input surface.  CMake runs fine on
it, but it dies in the make process with this error:

CMakeFiles/SurfaceToImage.dir/SurfaceToImage.cxx.o: In function
`vtkSmartPointer<vtkPolyDataToImageStencil>::New()':
/usr/local/include/vtk-5.7/vtkSmartPointer.h:113: undefined reference to
`vtkPolyDataToImageStencil::New()'

The smart pointer line is fundamentally no different than any other line of
code.  ???  Here's my code:
____________________________________________________________________________________________________________________
#include <vtkSmartPointer.h>
#include <vtkPolyData.h>
#include <vtkImageData.h>
#include <vtkPolyDataReader.h>
#include <vtkXMLImageDataWriter.h>
#include <vtkPolyDataToImageStencil.h>
#include <vtkImageStencil.h>
#include <vtkPointData.h>

/**
 * This program generates a sphere (closed surface, vtkPolyData) and
converts it into volume
 * representation (vtkImageData) where the foreground voxels are 1 and the
background voxels are
 * 0. Internally vtkPolyDataToImageStencil is utilized. The resultant image
is saved to disk
 * in XML image data format.
 */
int main (int argc, char *argv[])
{
  if( argc < 3 )
  {
    std::cerr << "Usage: " << argv[0];
    std::cerr << " inputImageFile outputImageFile " << std::endl;
    return EXIT_FAILURE;
  }

  vtkSmartPointer<vtkPolyDataReader> surfaceSource =
vtkSmartPointer<vtkPolyDataReader>::New();
  surfaceSource->SetFileName( argv[1] );
  vtkSmartPointer<vtkPolyData> pd = surfaceSource->GetOutput();
  surfaceSource->Update();

  vtkSmartPointer<vtkImageData> whiteImage =
vtkSmartPointer<vtkImageData>::New();
  double bounds[6];
  pd->GetBounds(bounds);
  double spacing[3]; // desired volume spacing
  spacing[0] = 0.5;
  spacing[1] = 0.5;
  spacing[2] = 0.5;
  whiteImage->SetSpacing(spacing);

  // compute dimensions
  int dim[3];
  for (int i = 0; i < 3; i++)
  {
    dim[i] = static_cast<int>(ceil((bounds[i * 2 + 1] - bounds[i * 2]) /
spacing[i]));
  }
  whiteImage->SetDimensions(dim);
  whiteImage->SetExtent(0, dim[0] - 1, 0, dim[1] - 1, 0, dim[2] - 1);

  double origin[3];
  // NOTE: I am not sure whether or not we had to add some offset!
  origin[0] = bounds[0];// + spacing[0] / 2;
  origin[1] = bounds[2];// + spacing[1] / 2;
  origin[2] = bounds[4];// + spacing[2] / 2;
  whiteImage->SetOrigin(origin);

  whiteImage->SetScalarTypeToUnsignedChar();
  whiteImage->AllocateScalars();

  // fill the image with foreground voxels:
  unsigned char inval = 255;
  unsigned char outval = 0;
  vtkIdType count = whiteImage->GetNumberOfPoints();
  for (vtkIdType i = 0; i < count; ++i)
  {
    whiteImage->GetPointData()->GetScalars()->SetTuple1(i, inval);
  }

  // polygonal data --> image stencil:
  vtkSmartPointer<vtkPolyDataToImageStencil> pol2stenc =
vtkSmartPointer<vtkPolyDataToImageStencil>::New(); //<---- This is the
'problem'???
  pol2stenc->SetInput(pd);
  pol2stenc->SetOutputOrigin(origin);
  pol2stenc->SetOutputSpacing(spacing);
  pol2stenc->SetOutputWholeExtent(whiteImage->GetExtent());
  pol2stenc->Update();

  // cut the corresponding white image and set the background:
  vtkSmartPointer<vtkImageStencil> imgstenc =
vtkSmartPointer<vtkImageStencil>::New();
  imgstenc->SetInput(whiteImage);
  imgstenc->SetStencil(pol2stenc->GetOutput());
  imgstenc->ReverseStencilOff();
  imgstenc->SetBackgroundValue(outval);
  imgstenc->Update();

  vtkSmartPointer<vtkXMLImageDataWriter> writer =
vtkSmartPointer<vtkXMLImageDataWriter>::New();
  writer->SetFileName( argv[2] ); //"SphereVolume.mhd");
  writer->SetInput(imgstenc->GetOutput());
  writer->Write();

  return EXIT_SUCCESS;
}
____________________________________________________________________________________________________________________



On Tue, Apr 20, 2010 at 6:37 AM, David Doria
<daviddoria+vtk at gmail.com<daviddoria%2Bvtk at gmail.com>
> wrote:

> On Tue, Apr 20, 2010 at 3:34 AM, Lars Friedrich Lars <
> lars-friedrich at gmx.net> wrote:
>
>> Hello David Welch,
>>
>> referring to your original question "how to voxelize a simple surface" I
>> added a VTK wiki example:
>>
>> http://www.cmake.org/Wiki/VTK/Examples/PolyDataToImageData
>>
>> It shows how we could incorporate vtkPolyDataToImageStencil into
>> voxelizing a *closed* surface.
>>
>> Please NOTE: as I wrote on the example page I'm not sure about the image
>> origin. The visual agreement in the sphere-image-overlay (paraview) is not
>> perfect and I'm not sure whether or not we had to introduce some
>> spacing-dependent image origin offset to compensate that. If you use that
>> code and find out something in this direction it would be great if you
>> contributed it to that example!
>>
>> HTH,
>>
>> lars
>>
>>
> David W -
>
> Yes, the normals are definitely necessary. With a mesh, you can simply run
> vtkPolyDataNormals. With a point set, you'd have to use my
> vtkPointSetNormalEstimation and vtkPointSetNormalOrientation filters (part
> of the same surface reconstruction journal submission). Unfortunately I
> don't know what it means by the missing ImageData element.
>
> Lars -
>
> Looks great, thanks! Hopefully David W will play around with it and solve
> the mystery of the origin problem.
>
> Thanks,
>
> David
>



-- 
David Welch
Graduate Student
Dept. of Biomedical Engineering
University of Iowa
Lab: (319) 335-5279
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.vtk.org/pipermail/vtkusers/attachments/20100420/65e41ed5/attachment.htm>