[Insight-users] Reduce the number of control points in vtkContourWidget

Xiaopeng Yang yxp233 at postech.ac.kr
Mon Mar 21 21:00:45 EDT 2011


Hi John, Thanks very much. I implemented Kent's code, but my program got crashed after I ran it. I have no idea which part is wrong. If possible, could you please check my code briefly?



vtkSmartPointer<vtkContourWidget> ContourWidget[157];

class vtkSliderCallback2 : public vtkCommand
{
public:
  static vtkSliderCallback2 *New() 
    { return new vtkSliderCallback2; }
  void SetImageViewer(vtkImageViewer2 *viewer)
    { this->Viewer =  viewer; }
  virtual void Execute(vtkObject *caller, unsigned long , void* )
    {
      vtkSliderWidget *slider = static_cast<vtkSliderWidget *>(caller);
      vtkSliderRepresentation *sliderRepres = static_cast<vtkSliderRepresentation *>(slider->GetRepresentation());
      int pos = static_cast<int>(sliderRepres->GetValue());


	for (int j = 0; j < 157; j++)
      {
		  if (pos == j)
		  {
			  ContourWidget[j]->SetEnabled(true);
			  //ContourWidget2->SetEnabled(true);

		  }
		  else
		  {
			  ContourWidget[j]->SetEnabled(false);
			  //ContourWidget2->SetEnabled(false);
		  }
	  }
	  
    this->Viewer->SetSlice(pos);
    }
protected:
  vtkImageViewer2 *Viewer;
};

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

namespace
{
double PointsArea(vtkPoints *points)
{
  int       numPoints = points->GetNumberOfPoints();
  vtkIdType *ids = new vtkIdType[numPoints + 1];
  int       i;

  for ( i = 0; i < numPoints; i++ )
    {
    ids[i] = i;
    }
  ids[i] = 0;
  double normal[3];
  double rval( vtkPolygon::ComputeArea(points, numPoints, ids, normal) );
  delete[] ids;
  return rval;
}

double PolyDataArea(vtkPolyData *pd)
{
  assert(pd->GetNumberOfLines() == 1);
  vtkIdType npts(0);
  vtkIdType *pts(0);
  double    normal[3];
  pd->GetLines()->InitTraversal();
  pd->GetLines()->GetNextCell(npts, pts);
  return vtkPolygon::ComputeArea(pd->GetPoints(),
                                 npts,
                                 pts,
                                 normal);
}

inline
void
IdsMinusThisPoint(vtkIdType *ids, int PointCount, int i)
{
  int k = 0;

  for ( int j = 0; j < PointCount; j++ )
    {
    if ( j != i )
      {
      ids[k] = j;
      k++;
      }
    }
  // leaving off first point? 1 is index of first point
  // in the reduced polygon;
  ids[k] = i == 0 ? 1 : 0;
}

void
Cull(vtkPolyData *in, vtkPolyData *out)
{
  //
  // Algorithm:
  // Error = 0
  // Create BTPolygon from points in input. Compute Original Area.
  //
  // while Error < MaxError
  //
  //    foreach vertex in input
  //      Create Polygon from input, minus this vertex
  //      Compute area.
  //      subtract area from original area.
  //    endfor
  //    find minimum error vertex, remove it.
  //    error = minimum error
  //

  //

  // do an initial point count reduction based on
  // curvature through vertices of the polygons.
  vtkPolyData *in2 = vtkPolyData::New();

  ReducePolyData2D(in, in2, 1);

  double originalArea( PolyDataArea(in2) );

  //
  // SWAG numbers -- accept
  // area change of 0.5%,
  // regard 0.005% as the same as zero
  double maxError = originalArea * 0.005;
  double errEpsilon = maxError * 0.001;

  vtkPoints *curPoints = vtkPoints::New();
  curPoints->DeepCopy( in2->GetPoints() );
  int       PointCount = curPoints->GetNumberOfPoints();
  vtkIdType *ids = new vtkIdType[PointCount];
  vtkIdType minErrorPointID = -1;

  for (;; )
    {
    double minError = 10000000.00;
    //
    // remove each point, one at a time and find the minimum error
    for ( int i = 0; i < PointCount; i++ )
      {
      // build id list, minus the current point;
      IdsMinusThisPoint(ids, PointCount, i);
      double normal[3];
      double curArea = vtkPolygon::ComputeArea(curPoints,
                                               PointCount - 1,
                                               ids,
                                               normal);
      double thisError = fabs(originalArea - curArea);
      if ( thisError < minError )
        {
        minError = thisError;
        minErrorPointID = i;
        // if the area error is absurdly low, just get rid of
        // this point and move on.
        if ( thisError < errEpsilon )
          {
          break;
          }
        }
      }
    //
    // if we have a new winner for least important point
    if ( minError <= maxError )
      {
      vtkPoints *newPoints = vtkPoints::New();

      for ( int i = 0; i < PointCount; i++ )
        {
        if ( i == minErrorPointID )
          {
          continue;
          }
        double point[3];
        curPoints->GetPoint(i, point);
        newPoints->InsertNextPoint(point);
        }
      curPoints->Delete();
      curPoints = newPoints;
      --PointCount;
      }
    else
      {
      break;
      }
    }
  ConvertPointSequenceToPolyData(curPoints, 1, out);
  curPoints->Delete();
  in2->Delete();
  delete[] ids;
}
}

int main (int argc, char *argv[])
{
  typedef    float    InputPixelType;
  typedef    float    OutputPixelType;
  typedef itk::Image< InputPixelType,  3 >   InputImageType;
  typedef itk::Image< OutputPixelType, 3 >   OutputImageType;
  typedef itk::Image< InputPixelType,       2 >   SliceImageType;

  typedef itk::ExtractImageFilter<InputImageType, SliceImageType> ExtractFilterType;

  // Load DICOM files
  typedef itk::ImageSeriesReader< InputImageType >     ReaderType;
  ReaderType::Pointer reader1 = ReaderType::New();
  ReaderType::Pointer reader2 = ReaderType::New();

  typedef itk::GDCMImageIO                        ImageIOType;
  typedef itk::GDCMSeriesFileNames                NamesGeneratorType;

  ImageIOType::Pointer gdcmIO = ImageIOType::New();
  //load original CT
  NamesGeneratorType::Pointer namesGenerator1 = NamesGeneratorType::New();
  namesGenerator1->SetInputDirectory( "C:/Users/User/Desktop/CT data/SHS/Diffusion_filter" );
  const ReaderType::FileNamesContainer & filenames1 = 
                            namesGenerator1->GetInputFileNames();
  reader1->SetImageIO( gdcmIO );
  reader1->SetFileNames( filenames1 );
  reader1->Update();

  typedef itk::ImageToVTKImageFilter<InputImageType>FilterType;
  FilterType::Pointer connector1 = FilterType::New();
  
  connector1->SetInput(reader1->GetOutput());

  vtkSmartPointer<vtkLookupTable> table1 = vtkSmartPointer<vtkLookupTable>::New();
  table1->SetRange(0, 255); // image intensity range
  table1->SetValueRange(0.0, 1.0); // from black to white
  table1->SetSaturationRange(0.0, 0.0); // no color saturation 
  table1->SetRampToLinear();
  table1->Build();

  vtkSmartPointer<vtkImageMapToColors> color1 = vtkSmartPointer<vtkImageMapToColors>::New();
  color1->SetLookupTable(table1);
  color1->SetInput(connector1->GetOutput());
  color1->Update();
 
  //load segmented CT
  NamesGeneratorType::Pointer namesGenerator2 = NamesGeneratorType::New();
  namesGenerator2->SetInputDirectory( "C:/Users/User/Desktop/CT data/SHS/Threshold_LS/Hole filled" );
  const ReaderType::FileNamesContainer & filenames2 = 
                            namesGenerator2->GetInputFileNames();
  reader2->SetImageIO( gdcmIO );
  reader2->SetFileNames( filenames2 );
  reader2->Update();
 
  FilterType::Pointer connector2 = FilterType::New();
  
  connector2->SetInput(reader2->GetOutput());
   
  InputImageType::Pointer inputImage = reader2->GetOutput();
  InputImageType::PointType origin = inputImage->GetOrigin();
  InputImageType::SpacingType spacing = inputImage->GetSpacing();
  InputImageType::SizeType  size = inputImage->GetLargestPossibleRegion().GetSize();
  
  InputImageType::RegionType extractRegion;

  InputImageType::SizeType extractSize(size);
  extractSize[2] = 0;

  InputImageType::IndexType       extractIndex;
  ExtractFilterType::Pointer extractFilter = ExtractFilterType::New();
  extractFilter->SetInput(inputImage);
  double minArea = ( spacing[0] * spacing[1] ) / 0.1;
  
  vtkImageViewer2 *imageViewer1 = vtkImageViewer2::New();
  imageViewer1->SetInput(color1->GetOutput());
  imageViewer1->SetColorLevel(127);
  imageViewer1->SetColorWindow(255);

  vtkRenderWindowInteractor *iren = vtkRenderWindowInteractor::New();

  imageViewer1->SetupInteractor(iren);
  imageViewer1->GetRenderWindow()->SetMultiSamples(0);
  imageViewer1->GetRenderWindow()->SetSize(600, 600);

  vtkSliderRepresentation2D *SliderRepres1 = vtkSliderRepresentation2D::New();
  int min = imageViewer1->GetSliceMin();
  int max = imageViewer1->GetSliceMax();
  SliderRepres1->SetMinimumValue(min);
  SliderRepres1->SetMaximumValue(max);
  SliderRepres1->SetValue(71);//static_cast<int>((min + max) / 2)
  //SliderRepres->SetTitleText("Slice");
  SliderRepres1->GetPoint1Coordinate()->SetCoordinateSystemToNormalizedDisplay();
  SliderRepres1->GetPoint1Coordinate()->SetValue(0.3, 0.05);
  SliderRepres1->GetPoint2Coordinate()->SetCoordinateSystemToNormalizedDisplay();
  SliderRepres1->GetPoint2Coordinate()->SetValue(0.7, 0.05);
  SliderRepres1->SetSliderLength(0.02);
  SliderRepres1->SetSliderWidth(0.03);
  SliderRepres1->SetEndCapLength(0.01);
  SliderRepres1->SetEndCapWidth(0.03);
  SliderRepres1->SetTubeWidth(0.005);
  SliderRepres1->SetLabelFormat("%3.0lf");
  SliderRepres1->SetTitleHeight(0.02);
  SliderRepres1->SetLabelHeight(0.02);

  vtkSliderWidget *SliderWidget1 = vtkSliderWidget::New();
  SliderWidget1->SetInteractor(iren);
  SliderWidget1->SetRepresentation(SliderRepres1);
  SliderWidget1->KeyPressActivationOff();
  SliderWidget1->SetAnimationModeToAnimate();
  SliderWidget1->SetEnabled(true);
  
  vtkSliderCallback2 *SliderCb1 = vtkSliderCallback2::New();
  SliderCb1->SetImageViewer(imageViewer1);
  SliderWidget1->AddObserver(vtkCommand::InteractionEvent, SliderCb1);  

  imageViewer1->SetSlice(static_cast<int>(SliderRepres1->GetValue()));
  imageViewer1->SetSliceOrientationToXY();

  vtkSmartPointer<vtkOrientedGlyphContourRepresentation> rep[157];

  for ( unsigned i = 0; i < 157; i++ )//size[2]
  {
	   
	   vtkExtractVOI * extract = vtkExtractVOI::New();
       extract->SetInput(connector2->GetOutput());
       extract->SetVOI(0,512,0,512,i,i);
  
       vtkContourFilter * contour = vtkContourFilter::New();
       contour->SetInput( extract->GetOutput() );
       contour->SetValue(0, 0.5);

	   vtkStripper *stripper = vtkStripper::New(); 
	   stripper->SetInput(contour->GetOutput()); 
	   stripper->Update(); 

	   vtkPolyData * pd = vtkPolyData::New();
	   pd = stripper->GetOutput();
   	   
	   vtkPolyData *pd3 = vtkPolyData::New();
	   if ( pd->GetPoints()->GetNumberOfPoints() > 0 )
       {
		   vtkPoints *pts = pd->GetPoints();
		   double    zPos = static_cast<double>( origin[2] )
			 + ( static_cast<double>( i ) * spacing[2] );
		   // std::cout << "Z Position " << zPos << std::endl;
		   for ( int j = 0; j < pd->GetNumberOfPoints(); j++ )
			 {
				 double point[3];
				 pts->GetPoint(j, point);
				 point[2] = zPos;
				 pts->SetPoint(j, point);
			 }

		   vtkCellArray *lines = pd->GetLines();
		   vtkIdType    numPoints;
		   vtkIdType    *points;
		  
		   while ( lines->GetNextCell(numPoints, points) != 0 )
			 {
				 vtkPoints *tmpPoints = vtkPoints::New();
				 for ( int j = 0; j < numPoints; j++ )
				   {
					   double point[3];
					   pts->GetPoint(points[j], point);
					   tmpPoints->InsertNextPoint(point);
				   }
				 if ( PointsArea(tmpPoints) < minArea )
				   {
					   tmpPoints->Delete();
					   continue;
				   }
				 vtkPolyData *pd2 = vtkPolyData::New();
				 ConvertPointSequenceToPolyData(tmpPoints, 1, pd2);
				 Cull(pd2, pd3);
			 }
	   }

	   ContourWidget[i] = vtkSmartPointer<vtkContourWidget>::New(); 
	   rep[i] = vtkSmartPointer<vtkOrientedGlyphContourRepresentation>::New(); 

	   ContourWidget[i]->SetInteractor(iren);
	   ContourWidget[i]->SetRepresentation(rep[i]);

	   rep[i]->GetProperty()->SetColor(1,0,1);
	   rep[i]->GetProperty()->SetPointSize(5);
	   rep[i]->GetLinesProperty()->SetColor(0,1,1);
	   rep[i]->GetLinesProperty()->SetLineWidth(3);

	   ContourWidget[i]->On();
	   ContourWidget[i]->Initialize(pd3);
	   ContourWidget[i]->SetEnabled(false);
  }

  imageViewer1->Render();
  imageViewer1->GetRenderer()->ResetCamera();

  iren->Initialize();
  iren->Start();
 
  return 0;
}






-----邮件原件-----
发件人: John Drescher [mailto:drescherjm at gmail.com] 
发送时间: 2011년 3월 22일 화요일 오전 3:28
收件人: Xiaopeng Yang
抄送: Karthik Krishnan; insight-users at itk.org; vtk
主题: Re: [Insight-users] [vtkusers] Reduce the number of control points in vtkContourWidget

> It works, but after decimation, the control points are not evenly
> distributed, which causes the change of the contour shape due to no control
> points at that part. So my question is how to make the points to be
> distributed evenly after decimation?
>

I experienced the same with the decimation for my usage (editing the
segmentation result for lungct). So I ended up using the node
reduction technique in the following example that Kent Williams (on
this list from time to time) worked on:
http://www.nitrc.org/projects/brainstracer/

The method is basically pulling out a node 1 by one and determining
the change in area of the filled polygon caused by the node removal.

John









More information about the Insight-users mailing list