[vtk-developers] Use vtkContour to find outline of mask in vtkImageData?

[dean.inglis at camris.ca]

Hi Kent, 

if you have extracted a set of points that are coincident with pixel
centers, you can try the routine below for reducing the number
of points...

Dean

// SAMPLE USE TO REDUCE CONTOUR POINTS COINCIDENT
// WITH PIXEL CENTERS

vtkPolyData* orig = const_cast<vtkPolyData*>( this->ContourRep->GetContourRepresentationAsPolyData() );

  vtkPolyData* temp = vtkPolyData::New();
  int closed = this->ContourRep->GetClosedLoop();

  if( vtkGeometryUtilities::ReducePolyData2D( orig, temp, closed ) )
    {
    this->ContourWidget->Initialize( temp );
    }

  temp->Delete();

// END SAMPLE USE


int vtkGeometryUtilities::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;
}

//---------------------------------------------------------------------------
// guarantees that only the minimum required number of points
// is used for open or closed line polydata
//
void vtkGeometryUtilities::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();
}
> 
> Kent:
> 
> Sorry for the long overdue reply.  Yes there is a way to do this.
> 
> The method you're looking for is
> 
>   vtkContourWidget::Initialize( polydata_containing_contour_control_points,
> StartingWidgetState )
> 
> What this does is to take a vtkPolyData (that contains 2 or more points and
> just one cell) and populate a contour from it.
> 
> StartingWidgetState as the documentation states can be
> vtkContourWidget::Define or vtkContourWidget::Manipulate. The former will
> let you go about adding new points. The latter puts it into manipulate
> state, (the same state you have on a vtkContourWidget, after you add the
> last point or close the loop).
> 
> ----------
> Now for the caveat:
> 
> I had a similar situation a while ago where I needed to go from a binary
> representation to a parametric representation, for editing purposes.
> If you run a vtkContourFilter and get a polydata out, you will have a
> polydata with several cells. You'll need to extract out each of those cells
> into seperate vtkPolydata. Take one of these cells. It will contain several
> points. If you initialize a contour widget from it, it will work, but it
> will have a few hundred control points and that doesn't look pretty. You can
> try it out :)
> 
> What you need is a downsampler, that picks every nth point, or better still
> a smart one that picks the ones that lie on high curvature.
> ----------
> 
> I just added a test that demonstrates initializing a contour widget from a
> user-supplied polydata. (a closed circle in this case).
> 
> Please let us know if you have issues.
> 
> --
> Karthik Krishnan
> R&D Engineer,
> Kitware Inc.
> 
> 
> On 1/9/08, kent williams <nkwmailinglists at gmail.com> wrote:
> >
> > I posted this to VTKUsers and got no suggestions.  Anyone here want to
> > comment?
> >
> > I need to go from a binary image, representing a 3D mask, to a set of
> > outlines on each slice of the binary image.
> >
> > It seems like something like this would work:
> >
> > Use vtkImageReslice to extract 2D slices from the 3D volume
> > For each slice, use vtkContourFilter to extract contours around areas
> > of set pixels.
> >
> > Right now I have a tracer that uses one or more vtkContourWidgets to
> > interactively generate a set of contours defining a mask.  I just
> > finished code that iterates through the contour set a slice at a time,
> > and generates a binary image volume where a non-zero pixel value is
> > inside the mask, and zero pixel value is outside the mask.  What I
> > want to do is reverse this process -- given a binary image, generate
> > the countour set.
> >
> > It looks like I should maybe use GenerateValues to generate contours,
> > but I don't understand the relationship between the output data of
> > vtkContour filter and the data i would need to initialize a
> > vtkContourWidget.  The few examples extant don't make that clear to me
> > either.
> >
> > Thanks!
> > _______________________________________________
> > vtk-developers mailing list
> > vtk-developers at vtk.org
> > http://www.vtk.org/mailman/listinfo/vtk-developers
> >
> 
>