[vtkusers] PolyData-Simplification after vtkContourFilter

[Lars Friedrich Lars]

hello,
I wrote a little application which reads a binary 3D-image, extracts the 3D-surface (vtkContourFilter) and then cuts it with a vtkPlane (vtkCutter) so that the corresponding 2D-contour in plane results.
In order to reduce the resulting vtkPolyData I incorporated an algorithm posted by dean inglis (http://www.vtk.org/pipermail/vtk-developers/2008-January/004913.html). the result looks strange: the polydata after simplification criss-crosses.
I don't think that this is a problem with the algorithm, it's rather a numbering problem of the contour filter's output.
Has anyone else observed this phenomenon? Is there a workaround/fix? Or am I misunderstanding some essential paradigms?
I attachted the code of the application below.
regards,
lars


#include "vtkRenderer.h"
#include "vtkRenderWindow.h"
#include "vtkRenderWindowInteractor.h"
#include "vtkPolyDataMapper.h"
#include "vtkActor.h"
#include "vtkProperty.h"
#include "vtkContourFilter.h"
#include "vtkMetaImageReader.h"
#include "vtkImageData.h"
#include "vtkDataSetMapper.h"
#include "vtkCutter.h"
#include "vtkPlane.h"

#include "vtkCellArray.h"
#include "vtkMath.h"





//---------------------------------------------------------------------------
// guarantees that only the minimum required number of points
// is used for open or closed line polydata
// by dean inglis
//
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();
}


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;
}




int main(int argc, char *argv[])
{

  vtkRenderer *aRenderer = vtkRenderer::New();
  vtkRenderWindow *renWin = vtkRenderWindow::New();
  
  renWin->AddRenderer(aRenderer);
  
  vtkRenderWindowInteractor *iren = vtkRenderWindowInteractor::New();
    iren->SetRenderWindow(renWin);

  vtkMetaImageReader *reader = vtkMetaImageReader::New();
  reader->SetFileName("/data/image-stuff/labelex/labelex.mhd");
  reader->Update();

  vtkContourFilter *iso = vtkContourFilter::New();
  int level = 1.0;

  iso->SetInput(reader->GetOutput());
  iso->SetValue(0, level);
  iso->UseScalarTreeOn();
  
  vtkPolyDataMapper *isoMapper = vtkPolyDataMapper::New();
  
  isoMapper->SetInput(iso->GetOutput());
  isoMapper->ScalarVisibilityOff();

  vtkActor *isoActor = vtkActor::New();
  
  isoActor->SetMapper(isoMapper);
  isoActor->GetProperty()->SetColor(0.8, 0.0, 0.0);
  isoActor->GetProperty()->SetOpacity(0.5);
  
  
  vtkPlane *cutPlane = vtkPlane::New();

  double *center = reader->GetOutput()->GetCenter();
  std::cout << "volume-center: " << center[0] << "," << center[1] << "," << center[2] << std::endl;
  double *planeOrigin = new double[3];
  double *off = new double[3];
  off[0] = 0;
  off[1] = 0;
  off[2] = 0;
  for (int i = 0; i < 3; ++i)
    planeOrigin[i] = center[i] + off[i];
  cutPlane->SetOrigin(planeOrigin);
  std::cout << "-> plane-origin: " << planeOrigin[0] << "," << planeOrigin[1] << "," << planeOrigin[2] << std::endl;
  double n1 = 0;
  double n2 = 0;
  double n3 = 1;
  cutPlane->SetNormal(n1, n2, n3);

  vtkCutter *cutter = vtkCutter::New();

  cutter->SetInput(iso->GetOutput());
  cutter->SetCutFunction(cutPlane);
  cutter->SetGenerateCutScalars(false);


  // --- point reduction test

  // ACTIVATE/DEACTIVATE the polydata-simplification
  int reductionOnOff = 1;

  vtkPolyData* originalContour = NULL;
  vtkPolyData* reducedContour = vtkPolyData::New();   
  int closed = 1;
  int reductionSuccess;

  cutter->Update();
  originalContour = cutter->GetOutput();

  if (ReducePolyData2D(originalContour, reducedContour, closed))
  { std::cout << "REDUCTION SUCCEEDED!" << std::endl;
    reductionSuccess = 1;
  }
  else
  { std::cout << "REDUCTION FAILED!" << std::endl;
    reductionSuccess = 0;
  }

  reductionSuccess = reductionSuccess && reductionOnOff;

  // --- \ point reduction test

  vtkDataSetMapper *cutMapper = vtkDataSetMapper::New();

  if (reductionSuccess)
    cutMapper->SetInput(reducedContour);
  else
    cutMapper->SetInput(originalContour);
  cutMapper->ScalarVisibilityOff();

  vtkActor *cutActor = vtkActor::New();

  cutActor->SetMapper(cutMapper);
  cutActor->GetProperty()->SetColor(0.0, 8.0, 0.0);
  cutActor->GetProperty()->SetLineWidth(2.0);
  cutActor->GetProperty()->SetOpacity(1.0);



  aRenderer->AddActor(isoActor);
  aRenderer->AddActor(cutActor);
  aRenderer->ResetCamera ();
  aRenderer->SetBackground(1,1,1);
  renWin->SetSize(1024, 768);
  aRenderer->ResetCameraClippingRange ();


  iren->Initialize();
  iren->Start(); 
  

  delete[] planeOrigin;
  cutActor->Delete();
  cutMapper->Delete();
  cutPlane->Delete();
  cutter->Delete();
  isoActor->Delete();
  isoMapper->Delete();
  iso->Delete();
  reader->Delete();
  iren->Delete();
  renWin->Delete();
  aRenderer->Delete();

  return 0;
}

-- 
Psssst! Schon vom neuen GMX MultiMessenger gehört? Der kann`s mit allen: http://www.gmx.net/de/go/multimessenger