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

[kent williams]

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;

}