Programming for Image Processing/Analysis and Visualization using The Visualization Toolkit

Programming for Image Processing/Analysis and Visualization using The Visualization Toolkit Week 6: Case Study 1 The Iterative Closest Point Algorithm Xenios Papademetris xenios.papademetris@yale.edu TAC (ex CAB) N119 5-6148 http://noodle.med.yale.edu/seminar/seminar.html

Schedule – Part 2 • Review of Part 1 and Course Overview • C++ Pointers/Classes, Object Oriented Programming • VTK – an object library • Adding new VTK Commands/Cmake • Image-to-image filters • Case Study I -- Iterative Closest Point surface matching • Case Study II – A Simple Segmentation Algorithm

Talk Outline • A Review of last few lectures • A review of the VTK Transformation Classes • Key Pieces • vtkLandmarkTransform • vtkCellLocator • The ICP Algorithm • The VTK implementation of ICP • vtkIterativeClosestPointTransform

Variables and Arrays • C++ has basic data types such as short, int, float, char, double etc. • The statements • int a; • float b[10]; • double c[5][5]; • define a single integer a, an one-dimensional array of • floats b: b[0] .. b[9] and a two-dimensional array of doubles c[0][0] .. c[4][4] (All array indices start at 0.) • Both a,b and c are implicitly allocated and will be de-allocated when the function containing them exits. Their sizes are fixed at compile time.

Dynamic Memory Allocation • Dynamic allocation allows the creation of arrays whose size is determined at runtime (e.g. loading an image whose size can vary). • It is one of the key to writing memory efficient programs. • It is, arguable, also the biggest source of problems and crashes in most C++ code. Most of the problems are caused by: • Accessing/Deleting arrays/objects before they are allocated and initialized. • Accessing/Deleting arrays/objects after they have been already deleted • Neglecting to delete arrays/objects

Allocating/De-allocating Pointers • Memory allocation is performed by the newcommand e.g. • int* a=new int; • a is NOT an integer. It is a pointer which contains a memory location. To access the value stored at the memory location pointed by a use the de-referencing operator * e.g. • *a = 1 • //This sets the value of the integer whose location is stored in a to 1 • When a is no longer needed the memory can be released using the delete command e.g. • delete a;

Allocating/De-allocating Arrays of Pointers • As before allocation is performed by the newcommand e.g. • int* a=new int[10]; // Allocate an array of 10 integers • a is a pointer which contains the memory location of the first element of the array. To access the values stored at the memory locations use: • a[0] = 1 ; a[3]=2; etc. • When a is no longer needed the memory can be released using the delete [] command e.g. • delete [] a; • Do not use the delete operator to delete arrays, it causes lots of problems use delete [] !!!

Class Design • Classes consists of: • Data members -- just like C structures • Methods – special functions which are embedded in the class • Members and methods can be declared as public which makes them accessible from outside the class or protected/private which makes them inaccessible from outside the class methods.

VTK and Macros • VTK uses certain macros to help with the parsing of the C++ header files and the creation of wrapper code for TCL/Java/Python. • The most common are: • vtkSetMacro(member,type); vtkGetMacro(member,type); • e.g. vtkSetMacro(Length,float)results in • void SetLength(float f) { Length=f; } • and vtkGetMacro(Length,float)results in • float GetLength() { return Length; } • vtkSetObjectMacro(object,type) and vtkGetObjectMacro(object,type) • are used to access/modify pointer objects derived from vtkObject

The vtkDataObject Branch

VTK Pipeline File Output Sources Filters Mappers Props vtkDataSet vtkDataSet • Two key filter methods: • Execute Information() • Execute()

Talk Outline • A Review of last few lectures • A review of the VTK Transformation Classes • Key Pieces • vtkLandmarkTransform • vtkCellLocator • The ICP algorithm • The VTK implementation of ICP • vtkIterativeClosestPointTransform

The vtkAbstractTransform Branch

vtkAbstractTransform • Abstract parentclass of all VTK transformations • Defines interface (but not the implementation) for most common operations e.g. • TransformPoint() • Identity() • Inverse() • TranformDerivative() • TransformNormalAtPoint() • vtkHomogeneousTransform provides a generic interface for homogeneous transformations, i.e. transformations which can be represented by multiplying a 4x4 matrix with a homogeneous coordinate.

The vtkHomogeneousTransform Branch

vtkIterativeClosestPoint Transform • Derived from vtkLinearTransform • It computes an “optimal” 4x4 matrix which defines the transformation between two data sets (stored as vtkDataSet) by iteratively estimating the correspondence between the point sets • From the VTK man page • “Match two surfaces using the iterative closest point (ICP) algorithm. The core of the algorithm is to match each vertex in one surface with the closest surface point on the other, then apply the transformation that modify one surface to best match the other (in a least square sense). This has to be iterated to get proper convergence of the surfaces. • Attention: • Use vtkTransformPolyDataFilter to apply the resulting ICP transform to your data. This class makes use of vtkLandmarkTransform internally to compute the best fit. Use the GetLandmarkTransform member to get a pointer to that transform and set its parameters. You might, for example, constrain the number of degrees of freedom of the solution (i.e. rigid body, similarity, etc.) by checking the vtkLandmarkTransform documentation for its SetMode member.”

vtkLandmarkTransform • Derived from vtkLinearTransform • It computes an “optimal” 4x4 matrix which defines the transformation between two sets of point (stored as vtkPoints) which have the same number of points and are assumed to be in correspondence. • From the VTK man page: • “A vtkLandmarkTransform is defined by two sets of landmarks, the transform computed gives the best fit mapping one onto the other, in a least squares sense. The indices are taken to correspond, so point 1 in the first set will get mapped close to point 1 in the second set, etc. Call SetSourceLandmarks and SetTargetLandmarks to specify the two sets of landmarks, ensure they have the same number of points. • Warning: • Whenever you add, subtract, or set points you must call Modified() on the vtkPoints object, or the transformation might not update.”

vtkLandmarkTransform II • Key Methods: • void SetSourceLandmarks (vtkPoints *points) • void SetTargetLandmarks (vtkPoints *points) • void SetModeToRigidBody () • void SetModeToSimilarity () • void SetModeToAffine () • From vtkHomogeneousTransform • vtkMatrix4x4 * GetMatrix() • void Update()

vtkLocator Tree From the man page: “vtkLocator is an abstract base class for spatial search objects, or locators. The principle behind locators is that they divide 3-space into small pieces (or "buckets") that can be quickly found in response to queries like point location, line intersection, or object-object intersection.”

vtkCellLocator • From the man page “vtkCellLocator is a spatial search object to quickly locate cells in 3D. vtkCellLocator uses a uniform-level octree subdivision, where each octant (an octant is also referred to as a bucket) carries an indication of whether it is empty or not, and each leaf octant carries a list of the cells inside of it. (An octant is not empty if it has one or more cells inside of it.) Typical operations are intersection with a line to return candidate cells, or intersection with another vtkCellLocator to return candidate cells.” • A cell is a triangle/voxel etc depending on the dataset. • Key methods: • void SetDataSet (vtkDataSet *) • void BuildLocator () • void FindClosestPoint (float x[3], float closestPoint[3], vtkIdType &cellId, int &subId, float &dist2)

Algorithm • Initialization • Take two points sets A and B • Initialize the transformation series with T0 as the identity. • Iteration k: • For each point in set A pj find the closest point in B to it qj • Estimate the best transformation Tkin a least squares sense that minimizes the distance between (pi and qi) • Compute the average distance d between the set of (pi and qi) If less than threshold stop. • Let pi=qi • Final transformation T is the concatenation of all Tks i.e. T=T3T2T1T0

vtkIterativeClosestPointTransform.h --1 • Standard stuff, some default constants and specification of inheritance etc • #ifndef __vtkIterativeClosestPointTransform_h • #define __vtkIterativeClosestPointTransform_h • #include "vtkLinearTransform.h" • #define VTK_ICP_MODE_RMS 0 • #define VTK_ICP_MODE_AV 1 • class vtkCellLocator; • class vtkLandmarkTransform; • class vtkDataSet; • class vtkIterativeClosestPointTransform : public vtkLinearTransform • { • public: • static vtkIterativeClosestPointTransform *New(); • vtkTypeMacro(vtkIterativeClosestPointTransform,vtkLinearTransform); • void PrintSelf(ostream& os, vtkIndent indent);

vtkIterativeClosestPointTransform.h --2 • Public Interface: all standard Get() macros ommitted. • // Description: Specify the source and target data sets. • void SetSource(vtkDataSet *source); • void SetTarget(vtkDataSet *target); • // Description: Set/Get a spatial locator for speeding up the search process. • // An instance of vtkCellLocator is used by default. • void SetLocator(vtkCellLocator *locator); • // Description: Set/Get the maximum number of iterations • vtkSetMacro(MaximumNumberOfIterations, int); • // Description: Force the algorithm to check the mean distance between two iteration. • vtkSetMacro(CheckMeanDistance, int); • // Description: Specify the mean distance mode. • void SetMeanDistanceModeToRMS() • void SetMeanDistanceModeToAbsoluteValue() • // Description: Set/Get the maximum mean distance between two iteration. If the mean distance is lower than this, the convergence stops. • vtkSetMacro(MaximumMeanDistance, float); • vtkGetMacro(MaximumMeanDistance, float);

vtkIterativeClosestPointTransform.h --3 • Public Interface continued : all standard Get() macros ommitted. • // Description: Set/Get the maximum number of landmarks sampled in your dataset. If your dataset is dense, then you will typically not need all the points to compute the ICP transform. • vtkSetMacro(MaximumNumberOfLandmarks, int); • // Description: Starts the process by translating source centroid to target centroid. • vtkSetMacro(StartByMatchingCentroids, int); • vtkBooleanMacro(StartByMatchingCentroids, int); • // Description: Get the internal landmark transform. Use it to constrain the number of degrees of freedom of the solution (i.e. rigid body, similarity, etc.). • vtkGetObjectMacro(LandmarkTransform,vtkLandmarkTransform); • // Description: Invert the transformation. This is done by switching the source and target. • void Inverse(); • // Description: Make another transform of the same type. • vtkAbstractTransform *MakeTransform();

vtkIterativeClosestPointTransform.h --4 • Protected Interface • // Description: Release source and target and locator (reference counting) • void ReleaseSource(void); • void ReleaseTarget(void); • void ReleaseLocator(void); • // Description Create default locator. Used to create one when none is specified. • void CreateDefaultLocator(void); • // Description: Get the MTime of this object also considering the locator. This is the Update test! • unsigned long int GetMTime(); • vtkIterativeClosestPointTransform(); • ~vtkIterativeClosestPointTransform(); • void InternalUpdate(); • // Description: This method does no type checking, use DeepCopy instead. • void InternalDeepCopy(vtkAbstractTransform *transform);

vtkIterativeClosestPointTransform.h --5 • Protected Interface continued • // Data Members • vtkDataSet* Source; • vtkDataSet* Target; • vtkCellLocator *Locator; • int MaximumNumberOfIterations; • int CheckMeanDistance; • int MeanDistanceMode; • float MaximumMeanDistance; • int MaximumNumberOfLandmarks; • int StartByMatchingCentroids; • int NumberOfIterations; • float MeanDistance; • vtkLandmarkTransform *LandmarkTransform; • // Stuff to keep off limits for users (copy constructor and = operator). // Not implemented. • private: • vtkIterativeClosestPointTransform(const vtkIterativeClosestPointTransform&); • void operator=(const vtkIterativeClosestPointTransform&); • }; • #endif

vtkIterativeClosestPointTransform.cxx • We will examine in detail the following methods • Constructor and Destructor • SetSource() and ReleaseSource() which are examples of obtaining external VTK objects and modifying the number of references to them to ensure that • The Source object does not get deleted accidentally (SetSource increases its reference number by 1) • The Source object does not stay in memory when it is no longer used (ReleaseSource decreases its reference number by 1) • Appropriate reference modification when the source object is changed! • The pairs SetTaget()/ GetTarget() SetLocator/ ReleaseLocator() are very similar to the SetSource()/ ReleaseSource() • CreateDefaultLocator() for creating the Cell Locator if no external one is specified! • InternalDeepCopy() for copying the object • And finally: • InternalUpdate() where the registration is computed.

vtkIterativeClosestPointTransform.cxx – 1Constructor and Destructor // The constructor sets all the default values vtkIterativeClosestPointTransform::vtkIterativeClosestPointTransform() : vtkLinearTransform() { this->Source = NULL; this->Target = NULL; this->Locator = NULL; this->LandmarkTransform = vtkLandmarkTransform::New(); this->MaximumNumberOfIterations = 50; this->CheckMeanDistance = 0; this->MeanDistanceMode = VTK_ICP_MODE_RMS; this->MaximumMeanDistance = 0.01; this->MaximumNumberOfLandmarks = 200; this->StartByMatchingCentroids = 0; this->NumberOfIterations = 0; this->MeanDistance = 0.0; } // The destructor releases all external inputs (Source,Target, Locator) and deletes the landmark transform object which was allocated in the constructor vtkIterativeClosestPointTransform::~vtkIterativeClosestPointTransform() { ReleaseSource(); ReleaseTarget(); ReleaseLocator(); this->LandmarkTransform->Delete(); }

vtkIterativeClosestPointTransform.cxx – 2SetSource() / ReleaseSource() void vtkIterativeClosestPointTransform::SetSource(vtkDataSet *source) { // Case 1 – the current source is the same as the one the method is trying to set: do nothing if (this->Source == source) return; // Otherwise first let go of the current source object (if it is not NULL) before taking a new source // if (this->Source) is equivalent to if (this->Source !=NULL) if (this->Source) this->ReleaseSource(); // If the input source is not NULL register it first and then set this->Source=source if (source) source->Register(this); this->Source = source; // The Modified() command tells the object that it needs to update its output this->Modified(); } void vtkIterativeClosestPointTransform::ReleaseSource(void) { if (this->Source) { this->Source->UnRegister(this); this->Source = NULL; } }

vtkIterativeClosestPointTransform.cxx – 3CreateDefaultLocator() & CreateInternalCopy() void vtkIterativeClosestPointTransform::CreateDefaultLocator() { if (this->Locator) { this->ReleaseLocator(); } this->Locator = vtkCellLocator::New(); } void vtkIterativeClosestPointTransform::InternalDeepCopy(vtkAbstractTransform *transform) { vtkIterativeClosestPointTransform *t = (vtkIterativeClosestPointTransform *)transform; this->SetSource(t->GetSource()); this->SetTarget(t->GetTarget()); this->SetLocator(t->GetLocator()); this->SetMaximumNumberOfIterations(t->GetMaximumNumberOfIterations()); this->SetCheckMeanDistance(t->GetCheckMeanDistance()); this->SetMeanDistanceMode(t->GetMeanDistanceMode()); this->SetMaximumMeanDistance(t->GetMaximumMeanDistance()); this->SetMaximumNumberOfLandmarks(t->GetMaximumNumberOfLandmarks()); this->Modified(); }

vtkIterativeClosestPointTransform.cxx – 4InternalUpdate() • InternalUpdate() gets called by the Update() function • The Update() function defined in a parent class checks for • The last modification time for this class • The last modification time for its inputs/parameters • It determines whether InternalUpdate() needs to be called. • To force execution use: • icp->Modified(); • icp->Update();

vtkIterativeClosestPointTransform.cxx – 4InternalUpdate() continued void vtkIterativeClosestPointTransform::InternalUpdate() { // Check source, target if (this->Source == NULL || !this->Source->GetNumberOfPoints()) { vtkErrorMacro(<<"Can't execute with NULL or empty input"); return; } if (this->Target == NULL || !this->Target->GetNumberOfPoints()) { vtkErrorMacro(<<"Can't execute with NULL or empty target"); return; } // Create locator this->CreateDefaultLocator(); this->Locator->SetDataSet(this->Target); this->Locator->SetNumberOfCellsPerBucket(1); this->Locator->BuildLocator();

vtkIterativeClosestPointTransform.cxx – 5InternalUpdate() continued // Create two sets of points to handle iteration step=sampling rate int step = 1; if (this->Source->GetNumberOfPoints() > this->MaximumNumberOfLandmarks) { step = this->Source->GetNumberOfPoints() / this->MaximumNumberOfLandmarks; } // Allocate some points. vtkIdType nb_points = this->Source->GetNumberOfPoints() / step; vtkPoints *points1 = vtkPoints::New(); points1->SetNumberOfPoints(nb_points); vtkPoints *closestp = vtkPoints::New(); closestp->SetNumberOfPoints(nb_points); vtkPoints *points2 = vtkPoints::New(); points2->SetNumberOfPoints(nb_points); // Fill with initial positions (sample dataset using step) vtkTransform *accumulate = vtkTransform::New(); accumulate->PostMultiply();

vtkIterativeClosestPointTransform.cxx – 6InternalUpdate() continued vtkIdType i; int j; float *p1, *p2; if (StartByMatchingCentroids) { float source_centroid[3] = {0,0,0}; for (i = 0; i < this->Source->GetNumberOfPoints(); i++) { p1 = this->Source->GetPoint(i); source_centroid[0] += p1[0]; source_centroid[1] += p1[1]; source_centroid[2] += p1[2]; } source_centroid[0] /= this->Source->GetNumberOfPoints(); source_centroid[1] /= this->Source->GetNumberOfPoints(); source_centroid[2] /= this->Source->GetNumberOfPoints(); float target_centroid[3] = {0,0,0}; for (i = 0; i < this->Target->GetNumberOfPoints(); i++) { p2 = this->Target->GetPoint(i); target_centroid[0] += p2[0]; target_centroid[1] += p2[1]; target_centroid[2] += p2[2]; } target_centroid[0] /= this->Target->GetNumberOfPoints(); target_centroid[1] /= this->Target->GetNumberOfPoints(); target_centroid[2] /= this->Target->GetNumberOfPoints(); accumulate->Translate(target_centroid[0] - source_centroid[0], target_centroid[1] - source_centroid[1], target_centroid[2] - source_centroid[2]); accumulate->Update(); for (i = 0, j = 0; i < nb_points; i++, j += step) accumulate->InternalTransformPoint(this->Source->GetPoint(j), points1->GetPoint(i)); } else { for (i = 0, j = 0; i < nb_points; i++, j += step) points1->SetPoint(i, this->Source->GetPoint(j)); }

Algorithm • Initialization • Take two points sets A and B • Initialize the transformation series with T0 as the identity. • Iteration k: • For each point in set A pj find the closest point in B to it qj • Estimate the best transformation Tkin a least squares sense that minimizes the distance between (pi and qi) • Compute the average distance d between the set of (pi and qi) If less than threshold stop. • Let pi=qi • Final transformation T is the concatenation of all Tks i.e. T=T3T2T1T0

vtkIterativeClosestPointTransform.cxx – 7InternalUpdate() continued // Go vtkIdType cell_id; int sub_id; float dist2, totaldist = 0; vtkPoints *temp, *a = points1, *b = points2; this->NumberOfIterations = 0; do { // Fill points with the closest points to each vertex in input for(i = 0; i < nb_points; i++) this->Locator->FindClosestPoint(a->GetPoint(i), closestp->GetPoint(i), cell_id, sub_id, dist2); // Build the landmark transform this->LandmarkTransform->SetSourceLandmarks(a); this->LandmarkTransform->SetTargetLandmarks(closestp); this->LandmarkTransform->Update(); // Concatenate accumulate->Concatenate(this->LandmarkTransform->GetMatrix());

vtkIterativeClosestPointTransform.cxx – 8InternalUpdate() continued // Check for Maximum Number of Iterations this->NumberOfIterations++; if (this->NumberOfIterations >= this->MaximumNumberOfIterations) break;// break out of do .. while loop // Move mesh and compute mean distance if needed if (this->CheckMeanDistance) { totaldist = 0.0; for(i = 0; i < nb_points; i++) { p1 = a->GetPoint(i); p2 = b->GetPoint(i); this->LandmarkTransform->InternalTransformPoint(p1, p2); totaldist += sqrt(vtkMath::Distance2BetweenPoints(p1, p2)); } this->MeanDistance = totaldist / (float)nb_points; } if (this->MeanDistance <= this->MaximumMeanDistance) break; } temp = a; a = b; b = temp; } while (1);

vtkIterativeClosestPointTransform.cxx – 9InternalUpdate() continued // Now recover accumulated result this->Matrix->DeepCopy(accumulate->GetMatrix()); accumulate->Delete(); points1->Delete(); closestp->Delete(); points2->Delete(); }

Next Week • Final Case Study: • Choices • Image-intensity based registration • Segmentation using snakes • Levelsets? (Marcel) • The idea would be to show how the design process will work rather than focusing too much on the underlying code • Highlight key VTK tools for making the implementation of these easier

Programming for Image Processing/Analysis and Visualization using The Visualization Toolkit