[vtkusers] Iterative closest point transform for point clouds - need some help
Amy Squillacote
amy.squillacote at kitware.com
Wed Nov 9 08:55:28 EST 2005
At 05:34 PM 11/8/2005, Sylvain Jaume wrote:
>Hi Vikas,
>
>You could define one vertex for every point in your input polydata.
>You may as well use the attached code.
>
>Cheers,
>Sylvain
>
>Vikas Singh wrote:
>
>>Hi Sylvain:
>>
>>I am sorry but your suggestions are not helping me solve the problem.
>>I just started on VTK yesterday, maybe I do not know VTK well enough
>>to understand what you're trying to tell me :-)
>>
>>For instance, your suggestion
>>vtkMaskPoints *maskPoints1 = vtkMaskPoints::New();
>>maskPoints1->SetInput(polyData1);
>>maskPoints1->SetMaximumNumberOfPoints(500);
>>maskPoints1->GenerateVerticesOn();
>>
>>_icpt->SetSource(maskPoints1);
Change the above line to the following.
_icpt->SetSource(maskPoints1->GetOutput());
>>does not compile. When I replace the last line by typecasting
>>maskPoints1 it gives me a seg-fault same as before. I am attaching the
>>code - its only about 50-60 lines. See if you can help me here.
>>
>>regards,
>>
>>Vikas.
>>
>>On 11/7/05, Sylvain Jaume <sylvain.jaume at kitware.com> wrote:
>>
>>
>>>Hi Vikas,
>>>
>>>You could try to create vertices with vtkMaskPoints before the ICP filter:
>>>
>>>vtkMaskPoints *maskPoints = vtkMaskPoints::New();
>>>maskPoints->SetInput(polyData1);
>>>maskPoints->SetMaximumNumberOfPoints(numPoints);
>>>maskPoints->GenerateVerticesOn();
>>>
>>>You don't need to cast polyData1. Use vtkPolyDataMaper instead of
>>>vtkDataSetMapper.
>>>
>>>Cheers,
>>>Sylvain
>>>
>>>Vikas Singh wrote:
>>>
>>>
>>>
>>>>Hi Sylvain:
>>>>
>>>>Unfortunately, that doesnt work. I mentioned that in my email
>>>>(attached below) that if I do not do a delaunay triangulation and call
>>>>icpt directly as
>>>>
>>>>_icpt->SetSource(polydata_1);
>>>>_icpt->SetTarget(polydata_2);
>>>>
>>>>it gives me a segmentation fault and the following error message:
>>>>ERROR: In /usr/local/VTK/Filtering/vtkCellLocator.cxx, line 1190
>>>>vtkCellLocator (0x8093b88): No cells to subdivide.
>>>>
>>>>Same error message even if I typecast it.
>>>>
>>>>But If I do a Delaunay, I avoid the error message and seg fault.
>>>>However, after applying the transformation on the source cloud set, I
>>>>cannot see the point set when it renders it.
>>>>See partial code below:
>>>>
>>>>// Iterative Closest Point Transform
>>>>vtkIterativeClosestPointTransform *_icpt =
>>>>vtkIterativeClosestPointTransform::New();
>>>>_icpt->GetLandmarkTransform()->SetModeToRigidBody();
>>>>_icpt->StartByMatchingCentroidsOn();
>>>>_icpt->SetMaximumNumberOfIterations(100.);
>>>>
>>>>// del_1 and del_2 are the delaunay triangulations for first and
>>>>second point clouds.
>>>>_icpt->SetSource((vtkDataSet *)del_1->GetOutput());
>>>>cout << "Calling icpt" << endl;
>>>>_icpt->SetTarget((vtkDataSet *)del_2->GetOutput());
>>>>
>>>>_icpt->Update();
>>>>
>>>>cout << "done icpt" << endl;
>>>>
>>>>// Now lets tranform the data
>>>>vtkTransformPolyDataFilter* transformer =
>>>>vtkTransformPolyDataFilter::New();
>>>>transformer->SetTransform(_icpt);
>>>>
>>>>transformer->SetInput((vtkPolyData*)polydata_1);
>>>>transformer->Update();
>>>>transformer->PrintSelf(std::cout,0);
>>>>
>>>>// map for displaying
>>>>
>>>>vtkDataSetMapper *surface_mapper_1 = vtkDataSetMapper::New();
>>>>
>>>>//surface_mapper_1->SetInput((vtkDataSet*)del_1->GetOutput());
>>>>// if i do this ^^^^^^^^^^^^^ it shows the untransformed first point cloud
>>>>// and the one below does not show anything (I want it to show the
>>>>// the transformed (aligned) point cloud
>>>>surface_mapper_1->SetInput(transformer->GetOutput());
>>>>
>>>>vtkDataSetMapper *surface_mapper_2 = vtkDataSetMapper::New();
>>>>surface_mapper_2->SetInput((vtkDataSet*)_icpt->GetTarget());
>>>>
>>>>// Create an actor for cloud 1.
>>>>vtkActor* actor_1 = vtkActor::New();
>>>>actor_1->SetMapper(surface_mapper_1);
>>>>actor_1->GetProperty()->SetRepresentationToPoints();
>>>>actor_1->GetProperty()->SetColor(1,1,1);
>>>>
>>>>// create an actor for cloud 2
>>>>vtkActor* actor_2 = vtkActor::New();
>>>>actor_2->SetMapper(surface_mapper_2);
>>>>actor_2->GetProperty()->SetRepresentationToPoints();
>>>>actor_2->GetProperty()->SetColor(0,1,1);
>>>>
>>>>......and then the rendering part.
>>>>
>>>>
>>>>Any suggestions ?
>>>>
>>>>--Vikas
>>>>
>>>>
>>>>On 11/7/05, Sylvain Jaume <sylvain.jaume at kitware.com> wrote:
>>>>
>>>>
>>>>
>>>>
>>>>>Hi Vikas,
>>>>>
>>>>>vtkIterativeClosestPointTransform gives a transform. You need to apply
>>>>>the resulting transform to your input dataset using
>>>>>vtkTransformPolyDataFilter. You don't need to use vtkDelaunay3D, nor to
>>>>>type cast your datasets.
>>>>>
>>>>>Cheers,
>>>>>Sylvain
>>>>>
>>>>>Vikas Singh wrote:
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>>Hello!
>>>>>>
>>>>>>I am a newbie to VTK and need some help with ICPT. Been trying all
>>>>>>sorts of things since today morning including the VTK User's guide and
>>>>>>the online documentation but nothing has worked so far. Seems to me
>>>>>>that the problem is simple for someone who knows VTK well.
>>>>>>
>>>>>>I am trying to randomly generate point clouds and align them
>>>>>>(as a start).
>>>>>>
>>>>>>vtkPolyData* polydata_1 = vtkPolyData::New();
>>>>>>vtkDelaunay3D* del_1 = vtkDelaunay3D::New();
>>>>>>vtkPoints* cloud_1 = vtkPoints::New();
>>>>>>
>>>>>>static float mypoint_1[3];
>>>>>>
>>>>>>for (int i = 0; i < 50; i++)
>>>>>> {
>>>>>> mypoint_1[0] = drand48()*50;
>>>>>> mypoint_1[1] = drand48()*50;
>>>>>> mypoint_1[2] = drand48()*50;
>>>>>>
>>>>>> cloud_1->InsertNextPoint(mypoint_1);
>>>>>> }
>>>>>>
>>>>>>polydata_1->SetPoints(cloud_1);
>>>>>>
>>>>>>int num_of_points_1 = polydata_1->GetPoints()->GetNumberOfPoints();
>>>>>>cout << "in all, num of points in cloud 1 is " <<
>>>>>>num_of_points_1 << endl;
>>>>>>
>>>>>>del_1->SetInput((vtkPointSet*)polydata_1);
>>>>>>del_1->SetTolerance(0.01);
>>>>>>del_1->PrintSelf(std::cout,0);
>>>>>>
>>>>>>I repeat the same for cloud 2 - variable names end with "_2".
>>>>>>
>>>>>>
>>>>>>Then I call ICPT as follows:
>>>>>>
>>>>>>vtkIterativeClosestPointTransform *_icpt =
>>>>>>vtkIterativeClosestPointTransform::New();
>>>>>>_icpt->GetLandmarkTransform()->SetModeToRigidBody();
>>>>>>_icpt->StartByMatchingCentroidsOn();
>>>>>>_icpt->SetMaximumNumberOfIterations(10);
>>>>>>
>>>>>>_icpt->SetSource((vtkDataSet *)del_1->GetOutput());
>>>>>>cout << "Calling icpt" << endl;
>>>>>>_icpt->SetTarget((vtkDataSet *)del_2->GetOutput());
>>>>>>
>>>>>>_icpt->Update();
>>>>>>
>>>>>>
>>>>>>To this it gives me an error
>>>>>>
>>>>>>ERROR: In
>>>>>>/usr/local/VTK/Hybrid/vtkIterativeClosestPointTransform.cxx, line 274
>>>>>>vtkIterativeClosestPointTransform (0x80650e0): Can't execute
>>>>>>with NULL or empty
>>>>>>input
>>>>>>
>>>>>>I think that up until the delaunay triangulation, it is working ok - I
>>>>>>rendered it and it is alright. If I do not do a delaunay triangulation
>>>>>>and typecast cloud_1 and cloud_2 directly it never returns from ICPT.
>>>>>>
>>>>>>any pointers, please?
>>>>>>
>>>>>>VS.
>>>>>>_______________________________________________
>>>>>>This is the private VTK discussion list.
>>>>>>Please keep messages on-topic. Check the FAQ at:
>>>>>>http://www.vtk.org/Wiki/VTK_FAQ
>>>>>>Follow this link to subscribe/unsubscribe:
>>>>>>http://www.vtk.org/mailman/listinfo/vtkusers
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>_______________________________________________
>>>>This is the private VTK discussion list.
>>>>Please keep messages on-topic. Check the FAQ at:
>>>>http://www.vtk.org/Wiki/VTK_FAQ
>>>>Follow this link to subscribe/unsubscribe:
>>>>http://www.vtk.org/mailman/listinfo/vtkusers
>>>>
>>>>
>>>>
>>>>
>>>------------------------------------------------------------------------
>>>
>>>// This example demonstrate the use of VTK data arrays as attribute
>>>// data as well as field data. It creates geometry (vtkPolyData) as
>>>// well as attribute data explicitly.
>>>
>>>// first include the required header files for the vtk classes we are using
>>>#include "vtkActor.h"
>>>#include "vtkCellArray.h"
>>>#include "vtkDoubleArray.h"
>>>#include "vtkFloatArray.h"
>>>#include "vtkIntArray.h"
>>>#include "vtkPointData.h"
>>>#include "vtkPoints.h"
>>>#include "vtkPointSet.h"
>>>#include "vtkPolyData.h"
>>>#include "vtkPolyDataMapper.h"
>>>#include "vtkRenderWindow.h"
>>>#include "vtkRenderWindowInteractor.h"
>>>#include "vtkRenderer.h"
>>>#include "vtkDelaunay3D.h"
>>>#include "vtkDataSetMapper.h"
>>>#include "vtkOutlineFilter.h"
>>>#include "vtkProperty.h"
>>>#include "vtkIterativeClosestPointTransform.h"
>>>#include "vtkLandmarkTransform.h"
>>>#include "vtkTransformPolyDataFilter.h"
>>>#include "vtkMaskPoints.h"
>>>
>>>using namespace std;
>>>
>>>int main( int argc, char *argv[] )
>>>{
>>> srand48(time(NULL));
>>> cout << "Program Starting" << endl;
>>>
>>> // Create the dataset. In this case, we create a vtkPolyData for
>>> first point cloud
>>> vtkPolyData* polydata_1 = vtkPolyData::New();
>>> vtkDelaunay3D* del_1 = vtkDelaunay3D::New();
>>> vtkPoints* cloud_1 = vtkPoints::New();
>>>
>>> static float mypoint_1[3];
>>>
>>> for (int i = 0; i < 500; i++)
>>> {
>>> mypoint_1[0] = drand48()*50;
>>> mypoint_1[1] = drand48()*50;
>>> mypoint_1[2] = drand48()*50;
>>>
>>> cloud_1->InsertNextPoint(mypoint_1);
>>> }
>>>
>>> polydata_1->SetPoints(cloud_1);
>>>
>>> int num_of_points_1 = polydata_1->GetPoints()->GetNumberOfPoints();
>>> cout << "in all, num of points in cloud 1 is " << num_of_points_1 << endl;
>>>
>>> del_1->SetInput((vtkPointSet*)polydata_1);
>>> del_1->SetTolerance(0.01);
>>> del_1->Update();
>>> del_1->SetAlpha(0);
>>> //del_1->BoundingTriangulationOn();
>>> //del_1->PrintSelf(std::cout,0);
>>>
>>> // finished dealing with the first point cloud
>>>
>>>
>>> // Create the dataset. In this case, we create a vtkPolyData for
>>> second point cloud
>>> vtkPolyData* polydata_2 = vtkPolyData::New();
>>> vtkDelaunay3D* del_2 = vtkDelaunay3D::New();
>>> vtkPoints* cloud_2 = vtkPoints::New();
>>>
>>> static float mypoint_2[3];
>>>
>>> for (int i = 0; i < 500; i++)
>>> {
>>> mypoint_2[0] = drand48()*50;
>>> mypoint_2[1] = drand48()*50;
>>> mypoint_2[2] = drand48()*50;
>>>
>>> cloud_2->InsertNextPoint(mypoint_2);
>>> }
>>>
>>> polydata_2->SetPoints(cloud_2);
>>>
>>> int num_of_points_2 = polydata_2->GetPoints()->GetNumberOfPoints();
>>> cout << "in all, num of points in cloud 2 is " << num_of_points_2 << endl;
>>>
>>> del_2->SetInput(polydata_2);
>>> del_2->SetTolerance(0.01);
>>> del_2->Update();
>>> del_2->SetAlpha(0);
>>> //del_1->BoundingTriangulationOn();
>>> //del_2->PrintSelf(std::cout,0);
>>>
>>> // finished dealing with the second point cloud
>>>
>>>
>>> // Iterative Closest Point Transform
>>> vtkIterativeClosestPointTransform *_icpt =
>>> vtkIterativeClosestPointTransform::New();
>>> _icpt->GetLandmarkTransform()->SetModeToRigidBody();
>>> _icpt->StartByMatchingCentroidsOn();
>>> //_icpt->SetMaximumMeanDistance(atof(MaxMeanDist->Text.c_str()));
>>> _icpt->SetMaximumNumberOfIterations(10.);
>>> //_icpt->SetMaximumNumberOfLandmarks(atof(MaxNLandmarks->Text.c_str()));
>>>
>>> vtkMaskPoints *maskPoints1 = vtkMaskPoints::New();
>>> maskPoints1->SetInput(polydata_1);
>>> maskPoints1->SetMaximumNumberOfPoints(500);
>>> maskPoints1->GenerateVerticesOn();
>>>
>>> vtkMaskPoints *maskPoints2 = vtkMaskPoints::New();
>>> maskPoints2->SetInput(polydata_2);
>>> maskPoints2->SetMaximumNumberOfPoints(500);
>>> maskPoints2->GenerateVerticesOn();
>>>
>>> //_icpt->SetSource((vtkDataSet *)del_1->GetOutput());
>>> _icpt->SetSource((vtkPolyData *)maskPoints1);
>>> cout << "Calling icpt" << endl;
>>> //_icpt->SetTarget((vtkDataSet *)del_2->GetOutput());
>>> _icpt->SetTarget((vtkPolyData *)maskPoints2);
>>>
>>> _icpt->SetMeanDistanceModeToAbsoluteValue();
>>> //_icpt->PrintSelf(std::cout, 0);
>>> _icpt->Update();
>>>
>>> cout << " Mean distance is " << _icpt->GetMeanDistance() << endl;
>>> cout << "done icpt" << endl;
>>>
>>> // Now lets tranform the data
>>> vtkTransformPolyDataFilter* transformer =
>>> vtkTransformPolyDataFilter::New();
>>> transformer->SetTransform(_icpt);
>>>
>>> //transformer->SetInput((vtkPolyData*)polydata_1);
>>> transformer->SetInput((vtkPolyData*)maskPoints1);
>>> transformer->Update();
>>> transformer->PrintSelf(std::cout,0);
>>>
>>>
>>> // map for displaying
>>>
>>> vtkDataSetMapper *surface_mapper_1 = vtkDataSetMapper::New();
>>> //surface_mapper_1->SetInput((vtkDataSet*)del_1->GetOutput());
>>>
>>> surface_mapper_1->SetInput((vtkDataSet*)_icpt->GetSource());
>>> //surface_mapper_1->SetInput(transformer->GetOutput());
>>> //surface_mapper_1->PrintSelf(std::cout,0);
>>>
>>> vtkDataSetMapper *surface_mapper_2 = vtkDataSetMapper::New();
>>> surface_mapper_2->SetInput((vtkDataSet*)_icpt->GetTarget());
>>> //surface_mapper_2->PrintSelf(std::cout,0);
>>>
>>>
>>> //cout << (_icpt->GetSource())->GetCenter() << endl;
>>>
>>> // Create an actor for cloud 1.
>>> vtkActor* actor_1 = vtkActor::New();
>>> actor_1->SetMapper(surface_mapper_1);
>>> //actor->GetProperty()->SetRepresentationToWireframe();
>>> actor_1->GetProperty()->SetRepresentationToPoints();
>>> actor_1->GetProperty()->SetColor(1,1,1);
>>>
>>> // create an actor for cloud 2
>>> vtkActor* actor_2 = vtkActor::New();
>>> actor_2->SetMapper(surface_mapper_2);
>>> //actor->GetProperty()->SetRepresentationToWireframe();
>>> actor_2->GetProperty()->SetRepresentationToPoints();
>>> actor_2->GetProperty()->SetColor(0,1,1);
>>>
>>> // Create the rendering objects.
>>> vtkRenderer* ren = vtkRenderer::New();
>>> ren->AddActor(actor_1);
>>> ren->AddActor(actor_2);
>>>
>>> vtkRenderWindow* renWin = vtkRenderWindow::New();
>>> renWin->AddRenderer(ren);
>>>
>>> vtkRenderWindowInteractor* iren = vtkRenderWindowInteractor::New();
>>> iren->SetRenderWindow(renWin);
>>> iren->Initialize();
>>> iren->Start();
>>>
>>> // pcoords->Delete();
>>> // points->Delete();
>>> // strips->Delete();
>>> // temperature->Delete();
>>> // vorticity->Delete();
>>> // polydata->Delete();
>>> // mapper->Delete();
>>> actor_1->Delete();
>>> actor_2->Delete();
>>> ren->Delete();
>>> renWin->Delete();
>>> iren->Delete();
>>>
>>> polydata_1->Delete();
>>> del_1->Delete();
>>> cloud_1->Delete();
>>>
>>> return 0;
>>>}
>>>
>>>
>>>
>
>
>#include "vtkActor.h"
>#include "vtkDoubleArray.h"
>#include "vtkPointData.h"
>#include "vtkPoints.h"
>#include "vtkPolyData.h"
>#include "vtkPolyDataMapper.h"
>#include "vtkRenderWindow.h"
>#include "vtkRenderWindowInteractor.h"
>#include "vtkRenderer.h"
>#include "vtkProperty.h"
>#include "vtkIterativeClosestPointTransform2.h"
>#include "vtkLandmarkTransform.h"
>#include "vtkTransformPolyDataFilter.h"
>#include "vtkMaskPoints.h"
>
>using namespace std;
>
>int main( int argc, char *argv[] )
>{
> srand48(time(NULL));
> cout << "Program Starting" << endl;
>
> vtkPolyData* polyData1 = vtkPolyData::New();
> vtkPoints* cloud1 = vtkPoints::New();
>
> int numPoints = 500;
> static float mypoint1[3];
>
> for (int i = 0; i < numPoints; i++)
> {
> mypoint1[0] = drand48()*50;
> mypoint1[1] = drand48()*50;
> mypoint1[2] = drand48()*50;
>
> cloud1->InsertNextPoint(mypoint1);
> }
>
> polyData1->SetPoints(cloud1);
>
> int num_of_points1 = polyData1->GetPoints()->GetNumberOfPoints();
> cout << "in all, num of points in cloud 1 is " << num_of_points1 << endl;
>
> vtkPolyData* polyData2 = vtkPolyData::New();
> vtkPoints* cloud2 = vtkPoints::New();
>
> static float mypoint2[3];
>
> for (int i = 0; i < numPoints; i++)
> {
> mypoint2[0] = drand48()*50;
> mypoint2[1] = drand48()*50;
> mypoint2[2] = drand48()*50;
>
> cloud2->InsertNextPoint(mypoint2);
> }
>
> polyData2->SetPoints(cloud2);
>
> int num_of_points2 = polyData2->GetPoints()->GetNumberOfPoints();
> cout << "in all, num of points in cloud 2 is " << num_of_points2 << endl;
>
> // Iterative Closest Point Transform
> vtkIterativeClosestPointTransform2
> *iterativeClosestPointTransform = vtkIterativeClosestPointTransform2::New();
> iterativeClosestPointTransform->SetSource(polyData1);
> iterativeClosestPointTransform->SetTarget(polyData2);
>
>iterativeClosestPointTransform->GetLandmarkTransform()->SetModeToRigidBody();
> iterativeClosestPointTransform->StartByMatchingCentroidsOn();
> iterativeClosestPointTransform->SetMaximumNumberOfIterations(10);
> iterativeClosestPointTransform->SetMeanDistanceModeToAbsoluteValue();
> iterativeClosestPointTransform->Update();
>
> cout << "Mean distance: " <<
> iterativeClosestPointTransform->GetMeanDistance() << endl;
>
> // tranform the data
> vtkTransformPolyDataFilter* transformPolyDataFilter =
> vtkTransformPolyDataFilter::New();
> transformPolyDataFilter->SetTransform(iterativeClosestPointTransform);
> transformPolyDataFilter->SetInput(polyData1);
> polyData1->Delete();
> transformPolyDataFilter->Update();
>
> vtkMaskPoints *maskPoints1 = vtkMaskPoints::New();
> maskPoints1->SetInput(transformPolyDataFilter->GetOutput());
> transformPolyDataFilter->Delete();
> maskPoints1->SetMaximumNumberOfPoints(numPoints);
> maskPoints1->GenerateVerticesOn();
> maskPoints1->Update();
>
> vtkMaskPoints *maskPoints2 = vtkMaskPoints::New();
> maskPoints2->SetInput(polyData2);
> polyData2->Delete();
> maskPoints2->SetMaximumNumberOfPoints(numPoints);
> maskPoints2->GenerateVerticesOn();
> maskPoints2->Update();
>
> vtkPolyDataMapper *polyDataMapper1 = vtkPolyDataMapper::New();
> polyDataMapper1->SetInput(maskPoints1->GetOutput());
> maskPoints1->Delete();
>
> vtkPolyDataMapper *polyDataMapper2 = vtkPolyDataMapper::New();
> polyDataMapper2->SetInput(maskPoints2->GetOutput());
> maskPoints2->Delete();
>
> // Create an actor for cloud 1.
> vtkActor* actor1 = vtkActor::New();
> actor1->SetMapper(polyDataMapper1);
> actor1->GetProperty()->SetColor(1,0,0);
>
> // create an actor for cloud 2
> vtkActor* actor2 = vtkActor::New();
> actor2->SetMapper(polyDataMapper2);
> actor2->GetProperty()->SetColor(0,1,0);
>
> // Create the rendering objects.
> vtkRenderer* ren = vtkRenderer::New();
> ren->AddActor(actor1);
> ren->AddActor(actor2);
>
> vtkRenderWindow* renWin = vtkRenderWindow::New();
> renWin->AddRenderer(ren);
>
> vtkRenderWindowInteractor* iren = vtkRenderWindowInteractor::New();
> iren->SetRenderWindow(renWin);
> iren->Initialize();
> iren->Start();
>
> actor1->Delete();
> actor2->Delete();
> ren->Delete();
> renWin->Delete();
> iren->Delete();
>
> return 0;
>}
>
>
>
>/*=========================================================================
>
> Program: Visualization Toolkit
> Module: $RCSfile: vtkIterativeClosestPointTransform2.cxx,v $
>
> Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
> All rights reserved.
> See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
>
> This software is distributed WITHOUT ANY WARRANTY; without even
> the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
> PURPOSE. See the above copyright notice for more information.
>
>=========================================================================*/
>#include "vtkIterativeClosestPointTransform2.h"
>
>#include "vtkPointLocator.h"
>#include "vtkDataSet.h"
>#include "vtkLandmarkTransform.h"
>#include "vtkMath.h"
>#include "vtkObjectFactory.h"
>#include "vtkPoints.h"
>#include "vtkTransform.h"
>
>vtkCxxRevisionMacro(vtkIterativeClosestPointTransform2, "$Revision: 1.14 $");
>vtkStandardNewMacro(vtkIterativeClosestPointTransform2);
>
>//----------------------------------------------------------------------------
>
>vtkIterativeClosestPointTransform2::vtkIterativeClosestPointTransform2()
> : 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;
>}
>
>//----------------------------------------------------------------------------
>
>const char *vtkIterativeClosestPointTransform2::GetMeanDistanceModeAsString()
>{
> if ( this->MeanDistanceMode == VTK_ICP_MODE_RMS )
> {
> return "RMS";
> }
> else
> {
> return "AbsoluteValue";
> }
>}
>
>//----------------------------------------------------------------------------
>
>vtkIterativeClosestPointTransform2::~vtkIterativeClosestPointTransform2()
>{
> ReleaseSource();
> ReleaseTarget();
> ReleaseLocator();
> this->LandmarkTransform->Delete();
>}
>
>//----------------------------------------------------------------------------
>
>void vtkIterativeClosestPointTransform2::SetSource(vtkDataSet *source)
>{
> if (this->Source == source)
> {
> return;
> }
>
> if (this->Source)
> {
> this->ReleaseSource();
> }
>
> if (source)
> {
> source->Register(this);
> }
>
> this->Source = source;
> this->Modified();
>}
>
>//----------------------------------------------------------------------------
>
>void vtkIterativeClosestPointTransform2::ReleaseSource(void) {
> if (this->Source)
> {
> this->Source->UnRegister(this);
> this->Source = NULL;
> }
>}
>
>//----------------------------------------------------------------------------
>
>void vtkIterativeClosestPointTransform2::SetTarget(vtkDataSet *target)
>{
> if (this->Target == target)
> {
> return;
> }
>
> if (this->Target)
> {
> this->ReleaseTarget();
> }
>
> if (target)
> {
> target->Register(this);
> }
>
> this->Target = target;
> this->Modified();
>}
>
>//----------------------------------------------------------------------------
>
>void vtkIterativeClosestPointTransform2::ReleaseTarget(void) {
> if (this->Target)
> {
> this->Target->UnRegister(this);
> this->Target = NULL;
> }
>}
>
>//----------------------------------------------------------------------------
>
>void vtkIterativeClosestPointTransform2::SetLocator(vtkPointLocator *locator)
>{
> if (this->Locator == locator)
> {
> return;
> }
>
> if (this->Locator)
> {
> this->ReleaseLocator();
> }
>
> if (locator)
> {
> locator->Register(this);
> }
>
> this->Locator = locator;
> this->Modified();
>}
>
>//----------------------------------------------------------------------------
>
>void vtkIterativeClosestPointTransform2::ReleaseLocator(void) {
> if (this->Locator)
> {
> this->Locator->UnRegister(this);
> this->Locator = NULL;
> }
>}
>
>//----------------------------------------------------------------------------
>
>void vtkIterativeClosestPointTransform2::CreateDefaultLocator() {
> if (this->Locator)
> {
> this->ReleaseLocator();
> }
>
> this->Locator = vtkPointLocator::New();
>}
>
>//------------------------------------------------------------------------
>
>unsigned long vtkIterativeClosestPointTransform2::GetMTime()
>{
> unsigned long result = this->vtkLinearTransform::GetMTime();
> unsigned long mtime;
>
> if (this->Source)
> {
> mtime = this->Source->GetMTime();
> if (mtime > result)
> {
> result = mtime;
> }
> }
>
> if (this->Target)
> {
> mtime = this->Target->GetMTime();
> if (mtime > result)
> {
> result = mtime;
> }
> }
>
> if (this->Locator)
> {
> mtime = this->Locator->GetMTime();
> if (mtime > result)
> {
> result = mtime;
> }
> }
>
> if (this->LandmarkTransform)
> {
> mtime = this->LandmarkTransform->GetMTime();
> if (mtime > result)
> {
> result = mtime;
> }
> }
>
> return result;
>}
>
>//----------------------------------------------------------------------------
>
>void vtkIterativeClosestPointTransform2::Inverse()
>{
> vtkDataSet *tmp1 = this->Source;
> this->Source = this->Target;
> this->Target = tmp1;
> this->Modified();
>}
>
>//----------------------------------------------------------------------------
>
>vtkAbstractTransform *vtkIterativeClosestPointTransform2::MakeTransform()
>{
> return vtkIterativeClosestPointTransform2::New();
>}
>
>//----------------------------------------------------------------------------
>
>void
>vtkIterativeClosestPointTransform2::InternalDeepCopy(vtkAbstractTransform
>*transform)
>{
> vtkIterativeClosestPointTransform2 *t =
> (vtkIterativeClosestPointTransform2 *)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();
>}
>
>//----------------------------------------------------------------------------
>
>void vtkIterativeClosestPointTransform2::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();
>
> // Create two sets of points to handle iteration
>
> int step = 1;
> if (this->Source->GetNumberOfPoints() > this->MaximumNumberOfLandmarks)
> {
> step = this->Source->GetNumberOfPoints() /
> this->MaximumNumberOfLandmarks;
> vtkDebugMacro(<< "Landmarks step is now : " << step);
> }
>
> vtkIdType nb_points = this->Source->GetNumberOfPoints() / step;
>
> // Allocate some points.
> // - closestp is used so that the internal state of
> LandmarkTransform remains
> // correct whenever the iteration process is stopped (hence its source
> // and landmark points might be used in a vtkThinPlateSplineTransform).
> // - points2 could have been avoided, but do not ask me why
> // InternalTransformPoint is not working correctly on my computer when
> // in and out are the same pointer.
>
> 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();
>
> vtkIdType i;
> int j;
> double p1[3], p2[3];
>
> if (StartByMatchingCentroids)
> {
> double source_centroid[3] = {0,0,0};
> for (i = 0; i < this->Source->GetNumberOfPoints(); i++)
> {
> this->Source->GetPoint(i, p1);
> 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();
>
> double target_centroid[3] = {0,0,0};
> for (i = 0; i < this->Target->GetNumberOfPoints(); i++)
> {
> this->Target->GetPoint(i, p2);
> 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)
> {
> double outPoint[3];
> accumulate->InternalTransformPoint(this->Source->GetPoint(j),
> outPoint);
> points1->SetPoint(i, outPoint);
> }
> }
> else
> {
> for (i = 0, j = 0; i < nb_points; i++, j += step)
> {
> points1->SetPoint(i, this->Source->GetPoint(j));
> }
> }
>
> // Go
>
> vtkIdType cell_id;
> int sub_id;
> double dist2, totaldist = 0;
> double outPoint[3];
>
> vtkIdType outPointId = -1;
> int targetNumPoints = this->Target->GetNumberOfPoints();
> 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++)
> {
> outPointId = this->Locator->FindClosestPoint(a->GetPoint(i));
>
> if ( outPointId < 0 || outPointId >= targetNumPoints )
> {
> vtkErrorMacro(<<"Cannot find the closest point");
> return;
> }
>
> this->Target->GetPoint(outPointId,outPoint);
> closestp->SetPoint(i, outPoint);
> }
>
> // Build the landmark transform
>
> this->LandmarkTransform->SetSourceLandmarks(a);
> this->LandmarkTransform->SetTargetLandmarks(closestp);
> this->LandmarkTransform->Update();
>
> // Concatenate (can't use this->Concatenate directly)
>
> accumulate->Concatenate(this->LandmarkTransform->GetMatrix());
>
> this->NumberOfIterations++;
>
> vtkDebugMacro(<< "Iteration: " << this->NumberOfIterations);
> if (this->NumberOfIterations >= this->MaximumNumberOfIterations)
> {
> break;
> }
>
> // Move mesh and compute mean distance if needed
>
> if (this->CheckMeanDistance)
> {
> totaldist = 0.0;
> }
>
> for(i = 0; i < nb_points; i++)
> {
> a->GetPoint(i, p1);
> this->LandmarkTransform->InternalTransformPoint(p1, p2);
> b->SetPoint(i, p2);
> if (this->CheckMeanDistance)
> {
> if (this->MeanDistanceMode == VTK_ICP_MODE_RMS)
> {
> totaldist += vtkMath::Distance2BetweenPoints(p1, p2);
> } else {
> totaldist += sqrt(vtkMath::Distance2BetweenPoints(p1, p2));
> }
> }
> }
>
> if (this->CheckMeanDistance)
> {
> if (this->MeanDistanceMode == VTK_ICP_MODE_RMS)
> {
> this->MeanDistance = sqrt(totaldist / (double)nb_points);
> } else {
> this->MeanDistance = totaldist / (double)nb_points;
> }
>
> vtkDebugMacro("Mean distance: " << this->MeanDistance);
> if (this->MeanDistance <= this->MaximumMeanDistance)
> {
> break;
> }
> }
>
> temp = a;
> a = b;
> b = temp;
> }
> while (1);
>
> // Now recover accumulated result
>
> this->Matrix->DeepCopy(accumulate->GetMatrix());
>
> accumulate->Delete();
> points1->Delete();
> closestp->Delete();
> points2->Delete();
>}
>
>//----------------------------------------------------------------------------
>
>void vtkIterativeClosestPointTransform2::PrintSelf(ostream& os,
>vtkIndent indent)
>{
> this->Superclass::PrintSelf(os,indent);
>
> if ( this->Source )
> {
> os << indent << "Source: " << this->Source << "\n";
> }
> else
> {
> os << indent << "Source: (none)\n";
> }
>
> if ( this->Target )
> {
> os << indent << "Target: " << this->Target << "\n";
> }
> else
> {
> os << indent << "Target: (none)\n";
> }
>
> if ( this->Locator )
> {
> os << indent << "Locator: " << this->Locator << "\n";
> }
> else
> {
> os << indent << "Locator: (none)\n";
> }
>
> os << indent << "MaximumNumberOfIterations: " <<
> this->MaximumNumberOfIterations << "\n";
> os << indent << "CheckMeanDistance: " << this->CheckMeanDistance << "\n";
> os << indent << "MeanDistanceMode: " <<
> this->GetMeanDistanceModeAsString() << "\n";
> os << indent << "MaximumMeanDistance: " <<
> this->MaximumMeanDistance << "\n";
> os << indent << "MaximumNumberOfLandmarks: " <<
> this->MaximumNumberOfLandmarks << "\n";
> os << indent << "StartByMatchingCentroids: " <<
> this->StartByMatchingCentroids << "\n";
> os << indent << "NumberOfIterations: " << this->NumberOfIterations << "\n";
> os << indent << "MeanDistance: " << this->MeanDistance << "\n";
> if(this->LandmarkTransform)
> {
> os << indent << "LandmarkTransform:\n";
> this->LandmarkTransform->PrintSelf(os, indent.GetNextIndent());
> }
>}
>
>/*=========================================================================
>
> Program: Visualization Toolkit
> Module: $RCSfile: vtkIterativeClosestPointTransform2.h,v $
>
> Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
> All rights reserved.
> See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
>
> This software is distributed WITHOUT ANY WARRANTY; without even
> the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
> PURPOSE. See the above copyright notice for more information.
>
>=========================================================================*/
>
>// .NAME vtkIterativeClosestPointTransform2 - Implementation of the
>ICP algorithm.
>// .SECTION Description
>// 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.
>// .SECTION Note
>// Use vtkTransformPolyDataFilter to apply the resulting ICP transform to
>// your data. You might also set it to your actor's user transform.
>// .SECTION Note
>// 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.
>// .SECTION see also
>// vtkLandmarkTransform
>
>
>#ifndef __vtkIterativeClosestPointTransform2_h
>#define __vtkIterativeClosestPointTransform2_h
>
>#include "vtkLinearTransform.h"
>
>#define VTK_ICP_MODE_RMS 0
>#define VTK_ICP_MODE_AV 1
>
>class vtkPointLocator;
>class vtkLandmarkTransform;
>class vtkDataSet;
>
>//class VTK_HYBRID_EXPORT vtkIterativeClosestPointTransform2 :
>public vtkLinearTransform
>class vtkIterativeClosestPointTransform2 : public vtkLinearTransform
>{
>public:
> static vtkIterativeClosestPointTransform2 *New();
>
>vtkTypeRevisionMacro(vtkIterativeClosestPointTransform2,vtkLinearTransform);
> void PrintSelf(ostream& os, vtkIndent indent);
>
> // Description:
> // Specify the source and target data sets.
> void SetSource(vtkDataSet *source);
> void SetTarget(vtkDataSet *target);
> vtkGetObjectMacro(Source, vtkDataSet);
> vtkGetObjectMacro(Target, vtkDataSet);
>
> // Description:
> // Set/Get a spatial locator for speeding up the search process.
> // An instance of vtkPointLocator is used by default.
> void SetLocator(vtkPointLocator *locator);
> vtkGetObjectMacro(Locator,vtkPointLocator);
>
> // Description:
> // Set/Get the maximum number of iterations
> vtkSetMacro(MaximumNumberOfIterations, int);
> vtkGetMacro(MaximumNumberOfIterations, int);
>
> // Description:
> // Get the number of iterations since the last update
> vtkGetMacro(NumberOfIterations, int);
>
> // Description:
> // Force the algorithm to check the mean distance between two iteration.
> vtkSetMacro(CheckMeanDistance, int);
> vtkGetMacro(CheckMeanDistance, int);
> vtkBooleanMacro(CheckMeanDistance, int);
>
> // Description:
> // Specify the mean distance mode. This mode expresses how the mean
> // distance is computed. The RMS mode is the square root of the average
> // of the sum of squares of the closest point distances. The Absolute
> // Value mode is the mean of the sum of absolute values of the closest
> // point distances.
> vtkSetClampMacro(MeanDistanceMode,int,
> VTK_ICP_MODE_RMS,VTK_ICP_MODE_AV);
> vtkGetMacro(MeanDistanceMode,int);
> void SetMeanDistanceModeToRMS()
> {this->SetMeanDistanceMode(VTK_ICP_MODE_RMS);}
> void SetMeanDistanceModeToAbsoluteValue()
> {this->SetMeanDistanceMode(VTK_ICP_MODE_AV);}
> const char *GetMeanDistanceModeAsString();
>
> // Description:
> // Set/Get the maximum mean distance between two iteration. If the mean
> // distance is lower than this, the convergence stops.
> vtkSetMacro(MaximumMeanDistance, double);
> vtkGetMacro(MaximumMeanDistance, double);
>
> // Description:
> // Get the mean distance between the last two iterations.
> vtkGetMacro(MeanDistance, double);
>
> // 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);
> vtkGetMacro(MaximumNumberOfLandmarks, int);
>
> // Description:
> // Starts the process by translating source centroid to target centroid.
> vtkSetMacro(StartByMatchingCentroids, int);
> vtkGetMacro(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();
>
>protected:
>
> // Description:
> // Release source and target
> void ReleaseSource(void);
> void ReleaseTarget(void);
>
> // Description:
> // Release locator
> 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.
> unsigned long int GetMTime();
>
> vtkIterativeClosestPointTransform2();
> ~vtkIterativeClosestPointTransform2();
>
> void InternalUpdate();
>
> // Description:
> // This method does no type checking, use DeepCopy instead.
> void InternalDeepCopy(vtkAbstractTransform *transform);
>
> vtkDataSet* Source;
> vtkDataSet* Target;
> vtkPointLocator *Locator;
> int MaximumNumberOfIterations;
> int CheckMeanDistance;
> int MeanDistanceMode;
> double MaximumMeanDistance;
> int MaximumNumberOfLandmarks;
> int StartByMatchingCentroids;
>
> int NumberOfIterations;
> double MeanDistance;
> vtkLandmarkTransform *LandmarkTransform;
>private:
> vtkIterativeClosestPointTransform2(const
> vtkIterativeClosestPointTransform2&); // Not implemented.
> void operator=(const vtkIterativeClosestPointTransform2&); //
> Not implemented.
>};
>
>#endif
>
>_______________________________________________
>This is the private VTK discussion list.
>Please keep messages on-topic. Check the FAQ at:
>http://www.vtk.org/Wiki/VTK_FAQ
>Follow this link to subscribe/unsubscribe:
>http://www.vtk.org/mailman/listinfo/vtkusers
More information about the vtkusers
mailing list