[vtkusers] Once more -- how do I get the outline contours of regions in a mask/binary image

Bill Lorensen bill.lorensen at gmail.com
Thu Mar 27 15:57:35 EDT 2008


Kent,

In the past, I did this sort of thing:
extract contours from an image with vtkContourFilter with an isovalue of .5.
convert disconnected line segments into polylines using vtkStripper

Bill

On Thu, Mar 27, 2008 at 12:17 PM, kent williams
<nkwmailinglists at gmail.com> wrote:
> Here's what I want:
>
> Input: Image where pixel == 0 is outside mask pixel != 0 is inside mask
> Output: Polygon/Polyline outline of all regions of nonzero pixels in image.
>
> Reason: My program uses vtkContourWidget to trace anatomical features.
>  I save the contours in an ITK xml file. I can also save a binary
> image -- i.e. an image file the same dimensions as the image being
> traced.  If a voxel is inside the traced contour, then it has a value
> of 255, otherwise zero.
>
> What I want to do is to read a binary image, and then find all the
> contours around the set pixels.  Basically the inverse of the
> operation above -- rather than using outline contours to generate a
> binary image, I want to use a binary image to get a set of outline
> contours.
>
> I've written what I think should work, below, but it looks really,
> really wrong when I load it -- it seems to be some woven agglomeration
> of many too many conour outlines.
>
> #include "MaskToContoursCLP.h"
> #include "itkBrains2MaskImageIO.h"
> #include "itkBrains2MaskImageIOFactory.h"
> #include "itkIO.h"
> #include "vtkKWImage.h"
> #include <vtkAppendPolyData.h>
> #include <vtkContourFilter.h>
> #include <vtkXMLPolyDataWriter.h>
> #include <vtkPolyData.h>
> #include <vtkImageData.h>
> #include <vtkDataObject.h>
> #include <vtkMath.h>
> #include <vtkIdTypeArray.h>
> #include <vtkCellArray.h>
> #include <itkExtractImageFilter.h>
>
> //
> // ConverPointSequenceToPolyData
> // does what the name says -- code borrowed from someon
> // on the VTK Mailing list
> //
> // logically what you end up is a collection of lines; at some point
> // we might want something different. At which point I'll have to dive
> // in and learn what the hell is going on with vtkCells.
> void ConvertPointSequenceToPolyData(vtkPoints* inPts,
>                                    const int& closed,
>                                    vtkPolyData* outPoly)
> {
>  if ( !inPts || !outPoly )
>    {
>    return;
>    }
>
>  int npts = inPts->GetNumberOfPoints();
>
>  if ( npts < 2 )
>    {
>    return;
>    }
>
>  double p0[3];
>  double p1[3];
>  inPts->GetPoint( 0, p0 );
>  inPts->GetPoint( npts - 1, p1 );
>  if ( p0[0] == p1[0] && p0[1] == p1[1] && p0[2] == p1[2] && closed )
>    {
>    --npts;
>    }
>
>  vtkPoints* temp = vtkPoints::New();
>  temp->SetNumberOfPoints( npts );
>  for ( int i = 0; i < npts; ++i )
>    {
>    temp->SetPoint( i, inPts->GetPoint( i ) );
>    }
>
>  vtkCellArray *cells = vtkCellArray::New();
>  cells->Allocate( cells->EstimateSize( npts + closed, 2 ) );
>
>  cells->InsertNextCell( npts + closed );
>
>  for ( int i = 0; i < npts; ++i )
>    {
>    cells->InsertCellPoint( i );
>    }
>
>  if ( closed )
>    {
>    cells->InsertCellPoint( 0 );
>    }
>
>  outPoly->SetPoints( temp );
>  temp->Delete();
>  outPoly->SetLines( cells );
>  cells->Delete();
> }
> //--------------------------------------------------------------------------
> // assumes a piecwise linear polyline with points at discrete locations
> //
> int
> ReducePolyData2D(vtkPolyData* inPoly,
>                 vtkPolyData* outPoly, const int& closed )
> {
>  if ( !inPoly || !outPoly )
>    {
>    return 0;
>    }
>
>  vtkPoints* inPts = inPoly->GetPoints();
>  if ( !inPts )
>    {
>    return 0;
>    }
>  int n = inPts->GetNumberOfPoints();
>  if ( n < 3 )
>    {
>    return 0;
>    }
>
>  double p0[3];
>  inPts->GetPoint( 0, p0 );
>  double p1[3];
>  inPts->GetPoint( n - 1, p1 );
>  bool minusNth = ( p0[0] == p1[0] && p0[1] == p1[1] && p0[2] == p1[2] );
>  if ( minusNth && closed )
>    {
>    --n;
>    }
>
>  struct frenet
>  {
>    double t[3];   // unit tangent vector
>    bool   state; // state of kappa: zero or non-zero  T/F
>  };
>
>  frenet* f;
>  f = new frenet[n];
>  double tL;
>  // calculate the tangent vector by forward differences
>  for ( int i = 0; i < n; ++i )
>    {
>    inPts->GetPoint( i, p0 );
>    inPts->GetPoint( ( i + 1 ) % n, p1 );
>    tL = sqrt( vtkMath::Distance2BetweenPoints( p0, p1 ) );
>    if ( tL == 0.0 ){ tL = 1.0; }
>    for ( int j = 0 ; j < 3; ++j )
>      {
>      f[i].t[j] = (p1[j] - p0[j]) / tL;
>      }
>    }
>
>  // calculate kappa from tangent vectors by forward differences
>  // mark those points that have very low curvature
>  double eps = 1.e-10;
>
>  for ( int i = 0; i < n; ++i )
>    {
>    f[i].state = ( fabs( vtkMath::Dot( f[i].t, f[( i + 1 ) % n].t ) - 1.0 )
>                   < eps );
>    }
>
>  vtkPoints* tempPts = vtkPoints::New();
>
>  // mark keepers
>  vtkIdTypeArray* ids = vtkIdTypeArray::New();
>
>  // for now, insist on keeping the first point for closure
>  ids->InsertNextValue( 0 );
>
>  for ( int i = 1; i < n; ++i )
>    {
>    bool pre = f[( i - 1 + n ) % n].state; // means fik != 1
>    bool cur = f[i].state;  // means fik = 1
>    bool nex = f[( i + 1 ) % n].state;
>
>    // possible vertex bend patterns for keep: pre cur nex
>    // 0 0 1
>    // 0 1 1
>    // 0 0 0
>    // 0 1 0
>
>    // definite delete pattern
>    // 1 1 1
>
>    bool keep = false;
>
>    if      (  pre &&  cur &&  nex ) { keep = false; }
>    else if ( !pre && !cur &&  nex ) { keep = true; }
>    else if ( !pre &&  cur &&  nex ) { keep = true; }
>    else if ( !pre && !cur && !nex ) { keep = true; }
>    else if ( !pre &&  cur && !nex ) { keep = true; }
>
>    if ( keep  ) // not a current sure thing but the preceding delete was
>      {
>      ids->InsertNextValue( i );
>      }
>    }
>
>  for ( int i = 0; i < ids->GetNumberOfTuples(); ++i )
>    {
>    tempPts->InsertNextPoint( inPts->GetPoint( ids->GetValue( i ) ) );
>    }
>
>  if ( closed )
>    {
>    tempPts->InsertNextPoint( inPts->GetPoint( ids->GetValue( 0 ) ) );
>    }
>
>  ConvertPointSequenceToPolyData( tempPts, closed, outPoly );
>
>  ids->Delete();
>  tempPts->Delete();
>  delete [] f;
>  return 1;
> }
>
> /*This program generates a signed distance image for a binary image.
> The intensity values inside the structues will be positive and outside
> the structure will be negative in the signed distance image. */
>
> int main( int argc, char * argv[] )
> {
>  PARSE_ARGS;
>  itk::Brains2MaskImageIOFactory::RegisterOneFactory();
>
>  const     unsigned int   Dimension = 3;
>  typedef float PixelType;
>  typedef itk::Image< PixelType,  Dimension >   ImageType;
>  typedef ImageType SliceImageType;
>  //  typedef itk::Image<PixelType,2> SliceImageType;
>  typedef itk::ExtractImageFilter<ImageType,SliceImageType> ExtractFilterType;
>
>  ImageType::Pointer inputImage(itkUtil::ReadImage<ImageType>(InputImage));
>
>  ImageType::SizeType size = inputImage->GetLargestPossibleRegion().GetSize();
>  //
>  // accumulate PolyData
>  vtkAppendPolyData *accumulator = vtkAppendPolyData::New();
>
>  //
>  // for each slice
>  for(unsigned i = 0; i < size[2]; i++)
>    {
>    //
>    // extract a slice from the image.
>    ImageType::RegionType extractRegion;
>
>    ImageType::SizeType extractSize(size);
>    extractSize[2] = 1;
>
>    ImageType::IndexType extractIndex;
>    extractIndex[0] = 0;
>    extractIndex[1] = 0;
>    extractIndex[2] = i;
>
>    extractRegion.SetSize(extractSize);
>    extractRegion.SetIndex(extractIndex);
>
>    ExtractFilterType::Pointer extractFilter = ExtractFilterType::New();
>    extractFilter->SetInput(inputImage);
>    extractFilter->SetExtractionRegion(extractRegion);
>    extractFilter->Update();
>
>    SliceImageType::Pointer sliceImage(extractFilter->GetOutput());
>
>    vtkKWImage *adapter = vtkKWImage::New();
>    adapter->SetITKImageBase(sliceImage);
>    vtkContourFilter *contourFilter = vtkContourFilter::New();
>    contourFilter->SetValue(0,0.0);
>    contourFilter->SetValue(1,1.0);
>    contourFilter->SetInput(vtkDataObject::SafeDownCast(adapter->GetVTKImage()));
>    contourFilter->Update();
>    vtkPolyData *pd = contourFilter->GetOutput();
>    vtkPolyData *pd2 = vtkPolyData::New();
>    ReducePolyData2D(pd,pd2,true);
>    accumulator->AddInput(pd2);
>    adapter->Delete();
>    contourFilter->Delete();
>    }
>  accumulator->Update();
>  vtkXMLPolyDataWriter *writer = vtkXMLPolyDataWriter::New();
>  writer->SetInput(accumulator->GetOutput());
>  writer->SetDataModeToAscii();
>  writer->SetFileName(OutputContours.c_str());
>  writer->Write();
>  writer->Delete();
>
>  return EXIT_SUCCESS;
>
> }
> _______________________________________________
> 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
>



More information about the vtkusers mailing list