lib/creaVtk/vtkHausdorffDistancePointSetFilter.cxx

   1 /*=========================================================================
   2
   3   Program:   Visualization Toolkit
   4   Module:    vtkHausdorffDistancePointSetFilter.cxx
   5
   6   Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
   7   All rights reserved.
   8   See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
   9
  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.
  13
  14 =========================================================================*/
  15 // Copyright (c) 2011 LTSI INSERM U642
  16 // All rights reserved.
  17 //
  18 // Redistribution and use in source and binary forms, with or without modification,
  19 // are permitted provided that the following conditions are met:
  20 //
  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.
  29 //
  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.
  40
  41 #include "vtkHausdorffDistancePointSetFilter.h"
  42
  43 #include "vtkDoubleArray.h"
  44 #include "vtkInformation.h"
  45 #include "vtkInformationVector.h"
  46 #include "vtkObjectFactory.h"
  47 #include "vtkPointData.h"
  48
  49 #include "vtkCellLocator.h"
  50 #include "vtkGenericCell.h"
  51 #include "vtkKdTreePointLocator.h"
  52 #include "vtkPointSet.h"
  53 #include "vtkSmartPointer.h"
  54
  55 vtkStandardNewMacro(vtkHausdorffDistancePointSetFilter);
  56
  57 //------------------------------------------------------------------------------
  58 vtkHausdorffDistancePointSetFilter::vtkHausdorffDistancePointSetFilter()
  59 {
  60   this->RelativeDistance[0] = 0.0;
  61   this->RelativeDistance[1] = 0.0;
  62   this->HausdorffDistance = 0.0;
  63
  64   this->SetNumberOfInputPorts(2);
  65   this->SetNumberOfInputConnections(0, 1);
  66   this->SetNumberOfInputConnections(1, 1);
  67
  68   this->SetNumberOfOutputPorts(2);
  69
  70   this->TargetDistanceMethod = POINT_TO_POINT;
  71 }
  72
  73 //------------------------------------------------------------------------------
  74 vtkHausdorffDistancePointSetFilter::~vtkHausdorffDistancePointSetFilter() = default;
  75
  76 //------------------------------------------------------------------------------
  77 int vtkHausdorffDistancePointSetFilter::RequestData(vtkInformation* vtkNotUsed(request),
  78   vtkInformationVector** inputVector, vtkInformationVector* outputVector)
  79 {
  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);
  85
  86   if (inInfoA == nullptr || inInfoB == nullptr)
  87   {
  88     return 0;
  89   }
  90
  91   // Get the input
  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()));
  96
  97   if (inputA->GetNumberOfPoints() == 0 || inputB->GetNumberOfPoints() == 0)
  98   {
  99     return 0;
 100   }
 101
 102   // Re-initialize the distances
 103   this->RelativeDistance[0] = 0.0;
 104   this->RelativeDistance[1] = 0.0;
 105   this->HausdorffDistance = 0.0;
 106
 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();
 113
 114   vtkSmartPointer<vtkCellLocator> cellLocatorA = vtkSmartPointer<vtkCellLocator>::New();
 115   vtkSmartPointer<vtkCellLocator> cellLocatorB = vtkSmartPointer<vtkCellLocator>::New();
 116
 117   if (this->TargetDistanceMethod == POINT_TO_POINT)
 118   {
 119     pointLocatorA->SetDataSet(inputA);
 120     pointLocatorA->BuildLocator();
 121     pointLocatorB->SetDataSet(inputB);
 122     pointLocatorB->BuildLocator();
 123   }
 124   else
 125   {
 126     cellLocatorA->SetDataSet(inputA);
 127     cellLocatorA->BuildLocator();
 128     cellLocatorB->SetDataSet(inputB);
 129     cellLocatorB->BuildLocator();
 130   }
 131
 132   double dist;
 133   double currentPoint[3];
 134   double closestPoint[3];
 135   vtkIdType cellId;
 136   vtkSmartPointer<vtkGenericCell> cell = vtkSmartPointer<vtkGenericCell>::New();
 137   int subId;
 138
 139   vtkSmartPointer<vtkDoubleArray> distanceAToB = vtkSmartPointer<vtkDoubleArray>::New();
 140   distanceAToB->SetNumberOfComponents(1);
 141   distanceAToB->SetNumberOfTuples(inputA->GetNumberOfPoints());
 142   distanceAToB->SetName("Distance");
 143
 144   vtkSmartPointer<vtkDoubleArray> distanceBToA = vtkSmartPointer<vtkDoubleArray>::New();
 145   distanceBToA->SetNumberOfComponents(1);
 146   distanceBToA->SetNumberOfTuples(inputB->GetNumberOfPoints());
 147   distanceBToA->SetName("Distance");
 148
 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++)
 152   {
 153     inputA->GetPoint(i, currentPoint);
 154     if (this->TargetDistanceMethod == POINT_TO_POINT)
 155     {
 156       vtkIdType closestPointId = pointLocatorB->FindClosestPoint(currentPoint);
 157       inputB->GetPoint(closestPointId, closestPoint);
 158     }
 159     else
 160     {
 161       cellLocatorB->FindClosestPoint(currentPoint, closestPoint, cell, cellId, subId, dist);
 162     }
 163
 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);
 168
 169     if (dist > this->RelativeDistance[0])
 170     {
 171       this->RelativeDistance[0] = dist;
 172     }
 173   }
 174
 175   for (int i = 0; i < inputB->GetNumberOfPoints(); i++)
 176   {
 177     inputB->GetPoint(i, currentPoint);
 178     if (this->TargetDistanceMethod == POINT_TO_POINT)
 179     {
 180       vtkIdType closestPointId = pointLocatorA->FindClosestPoint(currentPoint);
 181       inputA->GetPoint(closestPointId, closestPoint);
 182     }
 183     else
 184     {
 185       cellLocatorA->FindClosestPoint(currentPoint, closestPoint, cell, cellId, subId, dist);
 186     }
 187
 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);
 192
 193     if (dist > this->RelativeDistance[1])
 194     {
 195       this->RelativeDistance[1] = dist;
 196     }
 197   }
 198
 199   if (this->RelativeDistance[0] >= RelativeDistance[1])
 200   {
 201     this->HausdorffDistance = this->RelativeDistance[0];
 202   }
 203   else
 204   {
 205     this->HausdorffDistance = this->RelativeDistance[1];
 206   }
 207
 208   vtkSmartPointer<vtkDoubleArray> relativeDistanceAtoB = vtkSmartPointer<vtkDoubleArray>::New();
 209   relativeDistanceAtoB->SetNumberOfComponents(1);
 210   relativeDistanceAtoB->SetName("RelativeDistanceAtoB");
 211   relativeDistanceAtoB->InsertNextValue(RelativeDistance[0]);
 212
 213   vtkSmartPointer<vtkDoubleArray> relativeDistanceBtoA = vtkSmartPointer<vtkDoubleArray>::New();
 214   relativeDistanceBtoA->SetNumberOfComponents(1);
 215   relativeDistanceBtoA->SetName("RelativeDistanceBtoA");
 216   relativeDistanceBtoA->InsertNextValue(RelativeDistance[1]);
 217
 218   vtkSmartPointer<vtkDoubleArray> hausdorffDistanceFieldDataA =
 219     vtkSmartPointer<vtkDoubleArray>::New();
 220   hausdorffDistanceFieldDataA->SetNumberOfComponents(1);
 221   hausdorffDistanceFieldDataA->SetName("HausdorffDistance");
 222   hausdorffDistanceFieldDataA->InsertNextValue(HausdorffDistance);
 223
 224   vtkSmartPointer<vtkDoubleArray> hausdorffDistanceFieldDataB =
 225     vtkSmartPointer<vtkDoubleArray>::New();
 226   hausdorffDistanceFieldDataB->SetNumberOfComponents(1);
 227   hausdorffDistanceFieldDataB->SetName("HausdorffDistance");
 228   hausdorffDistanceFieldDataB->InsertNextValue(HausdorffDistance);
 229
 230   outputA->DeepCopy(inputA);
 231   outputA->GetPointData()->AddArray(distanceAToB);
 232   outputA->GetFieldData()->AddArray(relativeDistanceAtoB);
 233   outputA->GetFieldData()->AddArray(hausdorffDistanceFieldDataA);
 234
 235   outputB->DeepCopy(inputB);
 236   outputB->GetPointData()->AddArray(distanceBToA);
 237   outputB->GetFieldData()->AddArray(relativeDistanceBtoA);
 238   outputB->GetFieldData()->AddArray(hausdorffDistanceFieldDataB);
 239
 240   return 1;
 241 }
 242
 243 //------------------------------------------------------------------------------
 244 int vtkHausdorffDistancePointSetFilter::FillInputPortInformation(int port, vtkInformation* info)
 245 {
 246   // The input should be two vtkPointsSets
 247   if (port == 0)
 248   {
 249     info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkPointSet");
 250     return 1;
 251   }
 252   if (port == 1)
 253   {
 254     info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkPointSet");
 255     return 1;
 256   }
 257   return 0;
 258 }
 259
 260 //------------------------------------------------------------------------------
 261 void vtkHausdorffDistancePointSetFilter::PrintSelf(ostream& os, vtkIndent indent)
 262 {
 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";
 268 }
 269