Orfeo Toolbox  4.0
otbKeyPointSetsMatchingFilter.txx
Go to the documentation of this file.
1 /*=========================================================================
2 
3  Program: ORFEO Toolbox
4  Language: C++
5  Date: $Date$
6  Version: $Revision$
7 
8 
9  Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
10  See OTBCopyright.txt for details.
11 
12 
13  This software is distributed WITHOUT ANY WARRANTY; without even
14  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
15  PURPOSE. See the above copyright notices for more information.
16 
17 =========================================================================*/
18 #ifndef __otbKeyPointSetsMatchingFilter_txx
19 #define __otbKeyPointSetsMatchingFilter_txx
20 
22 
23 namespace otb
24 {
25 
26 template <class TPointSet, class TDistance>
29 {
30  this->SetNumberOfRequiredInputs(2);
31  m_UseBackMatching = false;
32  m_DistanceThreshold = 0.6;
33  // Object used to measure distance
34  m_DistanceCalculator = DistanceType::New();
35 }
36 
37 template <class TPointSet, class TDistance>
39 ::PointSetType *
42 {
43  return static_cast<const PointSetType *>(this->itk::ProcessObject::GetInput(0));
44 }
45 
46 template <class TPointSet, class TDistance>
47 void
49 ::SetInput1(const PointSetType * pointset)
50 {
51  this->itk::ProcessObject::SetNthInput(0, const_cast<PointSetType *>(pointset));
52 }
53 
54 template <class TPointSet, class TDistance>
56 ::PointSetType *
59 {
60  return static_cast<const PointSetType *>(this->itk::ProcessObject::GetInput(1));
61 }
62 
63 template <class TPointSet, class TDistance>
64 void
66 ::SetInput2(const PointSetType * pointset)
67 {
68  this->itk::ProcessObject::SetNthInput(1, const_cast<PointSetType *>(pointset));
69 }
70 
71 template <class TPointSet, class TDistance>
72 void
75 {
76 // std::cout<<"GenerateData()"<<std::endl;
77 
78  // Get the input pointers
79  const PointSetType * ps1 = this->GetInput1();
80  const PointSetType * ps2 = this->GetInput2();
81 
82  // Check if one of the pointsets is empty
83  if (ps1->GetNumberOfPoints() == 0 || ps2->GetNumberOfPoints() == 0)
84  {
85  itkExceptionMacro(<< "Empty input pointset !");
86  }
87 
88  // Get the output pointer
89  LandmarkListPointerType landmarks = this->GetOutput();
90 
91  // Define iterators on points and point data.
92  PointsIteratorType pIt = ps1->GetPoints()->Begin();
93  PointDataIteratorType pdIt = ps1->GetPointData()->Begin();
94 
95  // iterate on pointset 1
96  while (pdIt != ps1->GetPointData()->End()
97  && pIt != ps1->GetPoints()->End())
98  {
99  // Get point and point data at current location
100  bool matchFound = false;
101  unsigned int currentIndex = pIt.Index();
102  PointDataType data = pdIt.Value();
103  PointType point = pIt.Value();
104 
105  // These variables will hold the matched point and point data
106  PointDataType dataMatch;
107  PointType pointMatch;
108 
109  // call to the matching routine
110  NeighborSearchResultType searchResult1 = NearestNeighbor(data, ps2);
111 
112  // Check if the neighbor distance is lower than the threshold
113  if (searchResult1.second < m_DistanceThreshold)
114  {
115  // Get the matched point and point data
116  dataMatch = ps2->GetPointData()->GetElement(searchResult1.first);
117  pointMatch = ps2->GetPoints()->GetElement(searchResult1.first);
118 
119  // If the back matching option is on
120  if (m_UseBackMatching)
121  {
122  // Peform the back search
123  NeighborSearchResultType searchResult2 = NearestNeighbor(dataMatch, ps1);
124 
125  // Test if back search finds the same match
126  if (currentIndex == searchResult2.first)
127  {
128  matchFound = true;
129  }
130  }
131  else // else back matching
132  {
133  matchFound = true;
134  }
135  }
136 
137  // If we found a match, add the proper landmark
138  if (matchFound)
139  {
140  LandmarkPointerType landmark = LandmarkType::New();
141  landmark->SetPoint1(point);
142  landmark->SetPointData1(data);
143  landmark->SetPoint2(pointMatch);
144  landmark->SetPointData2(dataMatch);
145  landmark->SetLandmarkData(searchResult1.second);
146 
147  // Add the new landmark to the landmark list
148  landmarks->PushBack(landmark);
149  }
150  ++pdIt;
151  ++pIt;
152  }
153 }
154 
155 template <class TPointSet, class TDistance>
158 ::NearestNeighbor(const PointDataType& data1, const PointSetType * pointset)
159 {
160 // std::cout<<"Call to NearestNeighbor()"<<std::endl;
161 // Declare the result
163 
164  // Define iterators on points and point data.
165  PointDataIteratorType pdIt = pointset->GetPointData()->Begin();
166 
167  // local variables
168  unsigned int nearestIndex = 0;
169  double d1 = m_DistanceCalculator->Evaluate(data1, pdIt.Value());
170  ++pdIt;
171  double d2 = m_DistanceCalculator->Evaluate(data1, pdIt.Value());
172  ++pdIt;
173 
174  if (d1 > d2)
175  {
176  nearestIndex = 1;
177  }
178  // Initialize distances
179  double nearestDistance = std::min(d1, d2);
180  double secondNearestDistance = std::max(d1, d2);
181  double distanceValue;
182 
183  // iterate on the pointset
184  while (pdIt != pointset->GetPointData()->End())
185  {
186  // Evaluate the distance
187  distanceValue = m_DistanceCalculator->Evaluate(data1, pdIt.Value());
188 
189 // std::cout<<nearestIndex<<" "<<nearestDistance<<" "<<secondNearestDistance<<std::endl;
190 
191  // Check if this point is the nearest neighbor
192  if (distanceValue < nearestDistance)
193  {
194  secondNearestDistance = nearestDistance;
195  nearestDistance = distanceValue;
196  nearestIndex = pdIt.Index();
197 
198  }
199  // Else check if it is the second nearest neighbor
200  else if (distanceValue < secondNearestDistance)
201  {
202  secondNearestDistance = distanceValue;
203  }
204  ++pdIt;
205  }
206 
207  // Fill results
208  result.first = nearestIndex;
209  if (secondNearestDistance == 0)
210  {
211  result.second = 1;
212  }
213  else
214  {
215  result.second = nearestDistance / secondNearestDistance;
216  }
217 
218  // return the result
219  return result;
220 
221 }
222 
223 template <class TPointSet, class TDistance>
224 void
226 ::PrintSelf(std::ostream& os, itk::Indent indent) const
227 {
228  Superclass::PrintSelf(os, indent);
229 }
230 
231 } // end namespace otb
232 
233 #endif

Generated at Sat Mar 8 2014 16:02:51 for Orfeo Toolbox with doxygen 1.8.3.1