[vtkusers] using vtkitertiveclosestpointtransform

[Jonathan Bailleul]

> Moti Freiman wrote:
> 
> Hello All!
> I'm trying to develop registartion system using vtk. i found the
> vtkitertiveclosestpointtransform as implementation of icp algorithm.
> but to acclert the algorithm i want to use k-d-tree.
> a. i found that someone develop k-d-tree for vtk 5.0. can someone send
> me that code?
> b. i'm not sure about using the vkitertiveclosestpointtransform well.
> can someone send me c++ code example?
> 
> thanks to all
> moti

Maybe not the best example, but that works.


#include "vtkRenderer.h"
#include "vtkRenderWindow.h"
#include "vtkRenderWindowInteractor.h"
#include "vtkPolyData.h"
#include "vtkPoints.h"
#include "vtkCellArray.h"
#include "vtkFloatArray.h"
#include "vtkPolyDataMapper.h"
#include "vtkActor.h"
#include "vtkCamera.h"

#include "vtkProcrustesAlignmentFilter.h"
#include "vtkIterativeClosestPointTransform.h"
#include "vtkTransformPolyDataFilter.h"
#include "vtkPolyDataReader.h"
#include "vtkPolyDataWriter.h"
#include "vtkPolyDataMapper.h"
#include "vtkLandmarkTransform.h"
#include "vtkTransform.h"

#include "vtkPointSetToPointSetFilter.h"
#include "vtkDataSetWriter.h"

#include <assert.h>
#include "Vector.h"




/* original author   Allan Reinhold Kildeby */

/* GLOS (Graphic Library in Open Source), an ANSI Common Lisp OpenGL
subset.
Copyright (C) 2000 the GLOS development team
(http://glos.sourceforge.net) */

      

static void
usage(int argc, char **argv) 
{
  if (argc < 6)
    {
      printf("Aligns given 3d landmark files using ICP.\n");
      printf("Usage: %s <output files radix> <reference file> <auto
resize nonref files?> <0: similarity; 1: rigid-body> <to align files>
\n", argv[0]);

      //       printf("assumes file names of form radixNN only!\n");
      exit(1);
    }  
}



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


  int i;
  char* filename_out;

  char* out_radix = argv[1];
  char* reference_vol = argv[2];
  BOOL resize_mode = atoi(argv[3]);
  int rigid_mode = atoi(argv[4]);
  char** volumes = argv + 5;

  int nLandmarks = 66; ///100
  int nIterations = 100; 
  double dTolerance = 0.000001; /// 100;
  

  /* number of volumes besides reference */
  int nb_volumes = argc - 5;


  
  assert(filename_out = (char *)calloc(strlen(out_radix) + 10,
sizeof(char)));


  vtkPolyDataReader *ref_reader = vtkPolyDataReader::New();
  ref_reader -> SetFileName(reference_vol);
  ref_reader -> Update();

  vtkPolyDataReader *vol_reader = vtkPolyDataReader::New();
  
  vtkDataSetWriter *writer = vtkDataSetWriter::New(); 

  
  // Prepare mappers for getting bounds of objects
  vtkPolyDataMapper* ref_mapper = vtkPolyDataMapper::New();
  ref_mapper -> SetInput(ref_reader -> GetOutput());
  float ref_bounds[6]; ref_mapper -> GetBounds(ref_bounds); 
  vtkPolyDataMapper* vol_mapper = vtkPolyDataMapper::New();
  float vol_bounds[6];

  // retrieves norm of the diagonal of the bounding box of reference
  double ref_diagsize = Vector<float>(Point<float>(ref_bounds[0],
ref_bounds[2], ref_bounds[4]),
				      Point<float>(ref_bounds[1], ref_bounds[3],
ref_bounds[5])).norm();  
  double vol_diagsize, vol_factor;

  // scaling preparation
  vtkTransform* scaling = vtkTransform::New();
  vtkTransformPolyDataFilter* scaler =
vtkTransformPolyDataFilter::New();
  

  for (i = 0; i < nb_volumes; i++)
    {

      // compute a transformation associating current mesh with
reference one
      vtkIterativeClosestPointTransform *ICP =
vtkIterativeClosestPointTransform::New();
      vol_reader -> SetFileName(volumes[i]);
      vol_reader -> Update();

      
      if (resize_mode)
	{
	  vol_mapper -> SetInput(vol_reader -> GetOutput()); vol_mapper ->
Update();
	  vol_mapper -> GetBounds(vol_bounds); 

	  vol_diagsize = Vector<float>(Point<float>(vol_bounds[0],
vol_bounds[2], vol_bounds[4]),
				       Point<float>(vol_bounds[1], vol_bounds[3],
vol_bounds[5])).norm();  
	  vol_factor = (1.00 * ref_diagsize) / vol_diagsize; 
	  cout << endl << "Diagonals size:" << "ref: " << 1.0 * ref_diagsize 
<< "vol: " << vol_diagsize << "factor: " << vol_factor; 
	  
	  scaling -> Scale(vol_factor, vol_factor, vol_factor);
	  scaling -> Update();
	  scaler -> SetTransform(scaling);
	  scaler -> SetInput(vol_reader -> GetOutput());
	  scaler -> Update();
	  
	  ICP -> SetSource(scaler -> GetOutput());
	  //writer -> SetInput(scaler -> GetOutput());  writer ->
SetFileName("FECHIE.vtk"); writer -> Update();
	}
      else 
	{ICP -> SetSource(vol_reader -> GetOutput());}
      

      /* ICP parameters settings */
      ICP -> SetTarget(ref_reader -> GetOutput());
      
      if (rigid_mode)
	ICP -> GetLandmarkTransform() -> SetModeToRigidBody();
      else
	ICP -> GetLandmarkTransform() -> SetModeToSimilarity();
      
      ICP -> CheckMeanDistanceOn();
      // transform->SetMaximumMeanDistance(dTolerance);
      ICP -> SetMaximumNumberOfIterations(nIterations);
      ICP -> SetMaximumNumberOfLandmarks(nLandmarks);
      ICP -> StartByMatchingCentroidsOn();

      ICP -> SetMeanDistanceModeToRMS(); /// 

      ICP -> Update();



      // apply this ICP to current volume 
      vtkTransformPolyDataFilter *filter =
vtkTransformPolyDataFilter::New();
      filter -> SetTransform(ICP);
      if (resize_mode)
	filter -> SetInput(scaler -> GetOutput());
      else
	filter -> SetInput(vol_reader -> GetOutput());
      filter -> Update();

      // save transformed current volume      
      sprintf(filename_out, "%s%s%d.vtk", out_radix, ((i < 10) ? "0" :
""), i);

      writer -> SetInput(filter->GetOutput());
      writer -> SetFileName(filename_out);
      writer -> Update();

      filter -> Delete();
      ICP -> Delete();

	
    }

  return 0;
}


--
-----------------------------------
Jonathan BAILLEUL, Doctorant
GREYC Image - Université de Caen 
http://www.greyc.ismra.fr/~bailleul