1 /*=========================================================================
3 Program: Visualization Toolkit
4 Module: vtkHausdorffDistancePointSetFilter.cxx
6 Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
8 See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
10 This software is distributed WITHOUT ANY WARRANTY; without even
11 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12 PURPOSE. See the above copyright notice for more information.
14 =========================================================================*/
15 // Copyright (c) 2011 LTSI INSERM U642
16 // All rights reserved.
18 // Redistribution and use in source and binary forms, with or without modification,
19 // are permitted provided that the following conditions are met:
21 // * Redistributions of source code must retain the above copyright notice,
22 // this list of conditions and the following disclaimer.
23 // * Redistributions in binary form must reproduce the above copyright notice,
24 // this list of conditions and the following disclaimer in the documentation and/or
25 // other materials provided with the distribution.
26 // * Neither name of LTSI, INSERM nor the names
27 // of any contributors may be used to endorse or promote products derived from this
28 // software without specific prior written permission.
30 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ``AS IS''
31 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
32 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
33 // DISCLAIMED. IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE FOR ANY
34 // DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
35 // (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
36 // LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
37 // ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
38 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
39 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
41 #include "vtkHausdorffDistancePointSetFilter.h"
43 #include "vtkDoubleArray.h"
44 #include "vtkInformation.h"
45 #include "vtkInformationVector.h"
46 #include "vtkObjectFactory.h"
47 #include "vtkPointData.h"
49 #include "vtkCellLocator.h"
50 #include "vtkGenericCell.h"
51 #include "vtkKdTreePointLocator.h"
52 #include "vtkPointSet.h"
53 #include "vtkSmartPointer.h"
55 vtkStandardNewMacro(vtkHausdorffDistancePointSetFilter);
57 //------------------------------------------------------------------------------
58 vtkHausdorffDistancePointSetFilter::vtkHausdorffDistancePointSetFilter()
60 this->RelativeDistance[0] = 0.0;
61 this->RelativeDistance[1] = 0.0;
62 this->HausdorffDistance = 0.0;
64 this->SetNumberOfInputPorts(2);
65 this->SetNumberOfInputConnections(0, 1);
66 this->SetNumberOfInputConnections(1, 1);
68 this->SetNumberOfOutputPorts(2);
70 this->TargetDistanceMethod = POINT_TO_POINT;
73 //------------------------------------------------------------------------------
74 vtkHausdorffDistancePointSetFilter::~vtkHausdorffDistancePointSetFilter() = default;
76 //------------------------------------------------------------------------------
77 int vtkHausdorffDistancePointSetFilter::RequestData(vtkInformation* vtkNotUsed(request),
78 vtkInformationVector** inputVector, vtkInformationVector* outputVector)
80 // Get the info objects
81 vtkInformation* inInfoA = inputVector[0]->GetInformationObject(0);
82 vtkInformation* inInfoB = inputVector[1]->GetInformationObject(0);
83 vtkInformation* outInfoA = outputVector->GetInformationObject(0);
84 vtkInformation* outInfoB = outputVector->GetInformationObject(1);
86 if (inInfoA == nullptr || inInfoB == nullptr)
92 vtkPointSet* inputA = vtkPointSet::SafeDownCast(inInfoA->Get(vtkDataObject::DATA_OBJECT()));
93 vtkPointSet* inputB = vtkPointSet::SafeDownCast(inInfoB->Get(vtkDataObject::DATA_OBJECT()));
94 vtkPointSet* outputA = vtkPointSet::SafeDownCast(outInfoA->Get(vtkDataObject::DATA_OBJECT()));
95 vtkPointSet* outputB = vtkPointSet::SafeDownCast(outInfoB->Get(vtkDataObject::DATA_OBJECT()));
97 if (inputA->GetNumberOfPoints() == 0 || inputB->GetNumberOfPoints() == 0)
102 // Re-initialize the distances
103 this->RelativeDistance[0] = 0.0;
104 this->RelativeDistance[1] = 0.0;
105 this->HausdorffDistance = 0.0;
107 // TODO: using vtkStaticCellLocator, vtkStaticPointLocator is going to be much faster.
108 // Need to investigate and replace if appropriate.
109 vtkSmartPointer<vtkKdTreePointLocator> pointLocatorA =
110 vtkSmartPointer<vtkKdTreePointLocator>::New();
111 vtkSmartPointer<vtkKdTreePointLocator> pointLocatorB =
112 vtkSmartPointer<vtkKdTreePointLocator>::New();
114 vtkSmartPointer<vtkCellLocator> cellLocatorA = vtkSmartPointer<vtkCellLocator>::New();
115 vtkSmartPointer<vtkCellLocator> cellLocatorB = vtkSmartPointer<vtkCellLocator>::New();
117 if (this->TargetDistanceMethod == POINT_TO_POINT)
119 pointLocatorA->SetDataSet(inputA);
120 pointLocatorA->BuildLocator();
121 pointLocatorB->SetDataSet(inputB);
122 pointLocatorB->BuildLocator();
126 cellLocatorA->SetDataSet(inputA);
127 cellLocatorA->BuildLocator();
128 cellLocatorB->SetDataSet(inputB);
129 cellLocatorB->BuildLocator();
133 double currentPoint[3];
134 double closestPoint[3];
136 vtkSmartPointer<vtkGenericCell> cell = vtkSmartPointer<vtkGenericCell>::New();
139 vtkSmartPointer<vtkDoubleArray> distanceAToB = vtkSmartPointer<vtkDoubleArray>::New();
140 distanceAToB->SetNumberOfComponents(1);
141 distanceAToB->SetNumberOfTuples(inputA->GetNumberOfPoints());
142 distanceAToB->SetName("Distance");
144 vtkSmartPointer<vtkDoubleArray> distanceBToA = vtkSmartPointer<vtkDoubleArray>::New();
145 distanceBToA->SetNumberOfComponents(1);
146 distanceBToA->SetNumberOfTuples(inputB->GetNumberOfPoints());
147 distanceBToA->SetName("Distance");
149 // Find the nearest neighbors to each point and add edges between them,
150 // if they do not already exist and they are not self loops
151 for (int i = 0; i < inputA->GetNumberOfPoints(); i++)
153 inputA->GetPoint(i, currentPoint);
154 if (this->TargetDistanceMethod == POINT_TO_POINT)
156 vtkIdType closestPointId = pointLocatorB->FindClosestPoint(currentPoint);
157 inputB->GetPoint(closestPointId, closestPoint);
161 cellLocatorB->FindClosestPoint(currentPoint, closestPoint, cell, cellId, subId, dist);
164 dist = std::sqrt(std::pow(currentPoint[0] - closestPoint[0], 2) +
165 std::pow(currentPoint[1] - closestPoint[1], 2) +
166 std::pow(currentPoint[2] - closestPoint[2], 2));
167 distanceAToB->SetValue(i, dist);
169 if (dist > this->RelativeDistance[0])
171 this->RelativeDistance[0] = dist;
175 for (int i = 0; i < inputB->GetNumberOfPoints(); i++)
177 inputB->GetPoint(i, currentPoint);
178 if (this->TargetDistanceMethod == POINT_TO_POINT)
180 vtkIdType closestPointId = pointLocatorA->FindClosestPoint(currentPoint);
181 inputA->GetPoint(closestPointId, closestPoint);
185 cellLocatorA->FindClosestPoint(currentPoint, closestPoint, cell, cellId, subId, dist);
188 dist = std::sqrt(std::pow(currentPoint[0] - closestPoint[0], 2) +
189 std::pow(currentPoint[1] - closestPoint[1], 2) +
190 std::pow(currentPoint[2] - closestPoint[2], 2));
191 distanceBToA->SetValue(i, dist);
193 if (dist > this->RelativeDistance[1])
195 this->RelativeDistance[1] = dist;
199 if (this->RelativeDistance[0] >= RelativeDistance[1])
201 this->HausdorffDistance = this->RelativeDistance[0];
205 this->HausdorffDistance = this->RelativeDistance[1];
208 vtkSmartPointer<vtkDoubleArray> relativeDistanceAtoB = vtkSmartPointer<vtkDoubleArray>::New();
209 relativeDistanceAtoB->SetNumberOfComponents(1);
210 relativeDistanceAtoB->SetName("RelativeDistanceAtoB");
211 relativeDistanceAtoB->InsertNextValue(RelativeDistance[0]);
213 vtkSmartPointer<vtkDoubleArray> relativeDistanceBtoA = vtkSmartPointer<vtkDoubleArray>::New();
214 relativeDistanceBtoA->SetNumberOfComponents(1);
215 relativeDistanceBtoA->SetName("RelativeDistanceBtoA");
216 relativeDistanceBtoA->InsertNextValue(RelativeDistance[1]);
218 vtkSmartPointer<vtkDoubleArray> hausdorffDistanceFieldDataA =
219 vtkSmartPointer<vtkDoubleArray>::New();
220 hausdorffDistanceFieldDataA->SetNumberOfComponents(1);
221 hausdorffDistanceFieldDataA->SetName("HausdorffDistance");
222 hausdorffDistanceFieldDataA->InsertNextValue(HausdorffDistance);
224 vtkSmartPointer<vtkDoubleArray> hausdorffDistanceFieldDataB =
225 vtkSmartPointer<vtkDoubleArray>::New();
226 hausdorffDistanceFieldDataB->SetNumberOfComponents(1);
227 hausdorffDistanceFieldDataB->SetName("HausdorffDistance");
228 hausdorffDistanceFieldDataB->InsertNextValue(HausdorffDistance);
230 outputA->DeepCopy(inputA);
231 outputA->GetPointData()->AddArray(distanceAToB);
232 outputA->GetFieldData()->AddArray(relativeDistanceAtoB);
233 outputA->GetFieldData()->AddArray(hausdorffDistanceFieldDataA);
235 outputB->DeepCopy(inputB);
236 outputB->GetPointData()->AddArray(distanceBToA);
237 outputB->GetFieldData()->AddArray(relativeDistanceBtoA);
238 outputB->GetFieldData()->AddArray(hausdorffDistanceFieldDataB);
243 //------------------------------------------------------------------------------
244 int vtkHausdorffDistancePointSetFilter::FillInputPortInformation(int port, vtkInformation* info)
246 // The input should be two vtkPointsSets
249 info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkPointSet");
254 info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkPointSet");
260 //------------------------------------------------------------------------------
261 void vtkHausdorffDistancePointSetFilter::PrintSelf(ostream& os, vtkIndent indent)
263 this->Superclass::PrintSelf(os, indent);
264 os << indent << "HausdorffDistance: " << this->GetHausdorffDistance() << "\n";
265 os << indent << "RelativeDistance: " << this->GetRelativeDistance()[0] << ", "
266 << this->GetRelativeDistance()[1] << "\n";
267 os << indent << "TargetDistanceMethod: " << this->GetTargetDistanceMethodAsString() << "\n";