OTB  9.0.0
Orfeo Toolbox
otbImageToSIFTKeyPointSetFilter.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2022 Centre National d'Etudes Spatiales (CNES)
3  *
4  * This file is part of Orfeo Toolbox
5  *
6  * https://www.orfeo-toolbox.org/
7  *
8  * Licensed under the Apache License, Version 2.0 (the "License");
9  * you may not use this file except in compliance with the License.
10  * You may obtain a copy of the License at
11  *
12  * http://www.apache.org/licenses/LICENSE-2.0
13  *
14  * Unless required by applicable law or agreed to in writing, software
15  * distributed under the License is distributed on an "AS IS" BASIS,
16  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
17  * See the License for the specific language governing permissions and
18  * limitations under the License.
19  */
20 
21 #ifndef otbImageToSIFTKeyPointSetFilter_hxx
22 #define otbImageToSIFTKeyPointSetFilter_hxx
23 
24 #include "otbImageToSIFTKeyPointSetFilter.h"
25 
26 #include "itkMatrix.h"
27 #include "itkProcessObject.h"
28 
29 namespace otb
30 {
31 
35 template <class TInputImage, class TOutputPointSet>
36 ImageToSIFTKeyPointSetFilter<TInputImage, TOutputPointSet>::ImageToSIFTKeyPointSetFilter()
37  :
38  m_OctavesNumber(1),
39  m_ScalesNumber(3),
40  m_ExpandFactors(2),
41  m_ShrinkFactors(2),
42  m_DoGThreshold(0.03),
43  m_EdgeThreshold(10),
44  m_RatioEdgeThreshold(0),
45  m_GradientMagnitudeThreshold(0.2),
46  m_Sigma0(1.6),
47  m_Sigmak(0),
48  m_SigmaFactorOrientation(3),
49  m_SigmaFactorDescriptor(1.5),
50 
51  m_ExpandFilter(ExpandFilterType::New()),
52 
53  m_ValidatedKeyPoints(0),
54  m_DifferentSamplePoints(0),
55  m_DiscardedKeyPoints(0),
56  m_ChangeSamplePointsMax(2),
57 
58  m_HistogramGaussianWeights({
59  2.3771112282795414e-07, 3.8860734758633732e-07, 6.2655544995978937e-07, 9.9631120821413786e-07, 1.5624909838697011e-06, 2.4167238265599128e-06,
60  3.6865788528530121e-06, 5.5463469229192623e-06, 8.2295774080263437e-06, 1.2043009749602365e-05, 1.738119136656513e-05, 2.4740646513897326e-05,
61  3.4731980398846277e-05, 4.808781565748272e-05, 6.5664032975164266e-05, 8.8431512984476723e-05, 0.00011745555408931643, 0.00015386047198026335,
62  0.00019877765486783745, 0.00025327659834301937, 0.00031828015928190065, 0.00039446735551235698, 0.00048216931692246382, 0.00058126620279441276,
63  0.00069109471776775144, 0.00081037694122312908, 0.00093718121775182789, 0.0010689246133776746, 0.0012024238404411182, 0.0013339976954896103,
64  0.0014596192424447215, 0.0015751106965100009, 0.0016763688464699555, 0.0017596045720966803, 0.0018215772013714365, 0.0018598035923515156,
65  0.0018727231637146351, 0.0018598035923515156, 0.0018215772013714365, 0.0017596045720966803, 0.0016763688464699555, 0.0015751106965100009,
66  0.0014596192424447215, 0.0013339976954896103, 0.0012024238404411182, 0.0010689246133776746, 0.00093718121775182789, 0.00081037694122312908,
67  0.00069109471776775144, 0.00058126620279441276, 0.00048216931692246382, 0.00039446735551235698, 0.00031828015928190065, 0.00025327659834301937,
68  0.00019877765486783745, 0.00015386047198026335, 0.00011745555408931643, 8.8431512984476723e-05, 6.5664032975164266e-05, 4.808781565748272e-05,
69  3.4731980398846277e-05, 2.4740646513897326e-05, 1.738119136656513e-05, 1.2043009749602365e-05, 8.2295774080263437e-06, 5.5463469229192623e-06,
70  3.6865788528530121e-06, 2.4167238265599128e-06, 1.5624909838697011e-06, 9.9631120821413786e-07, 6.2655544995978937e-07, 3.8860734758633732e-07,
71  2.3771112282795414e-07})
72 {
73  m_Offsets[0][0] = -1;
74  m_Offsets[0][1] = -1;
75  m_Offsets[1][0] = -1;
76  m_Offsets[1][1] = 0;
77  m_Offsets[2][0] = -1;
78  m_Offsets[2][1] = 1;
79  m_Offsets[3][0] = 0;
80  m_Offsets[3][1] = -1;
81  m_Offsets[4][0] = 0;
82  m_Offsets[4][1] = 1;
83  m_Offsets[5][0] = 1;
84  m_Offsets[5][1] = -1;
85  m_Offsets[6][0] = 1;
86  m_Offsets[6][1] = 0;
87  m_Offsets[7][0] = 1;
88  m_Offsets[7][1] = 1;
89 }
90 
91 template <class TInputImage, class TOutputPointSet>
92 void ImageToSIFTKeyPointSetFilter<TInputImage, TOutputPointSet>::GenerateData()
93 {
94  // First, subsample the input image
95  InitializeInputImage();
96 
97  InputImagePointerType input = m_ExpandFilter->GetOutput();
98  typename InputImageType::PointType point;
99  typename InputImageType::IndexType index;
100  index[0] = 0;
101  index[1] = 0;
102 
103  m_ExpandFilter->GetOutput()->TransformIndexToPhysicalPoint(index, point);
104 
105  // for each octave, compute the difference of gaussian
106  unsigned int lOctave = 0;
107  m_Sigmak = std::pow(2, static_cast<double>(1 / (double)(m_ScalesNumber + 1)));
108  m_RatioEdgeThreshold = (m_EdgeThreshold + 1) * (m_EdgeThreshold + 1) / m_EdgeThreshold;
109 
110  for (lOctave = 0; lOctave != m_OctavesNumber; lOctave++)
111  {
112  m_DifferentSamplePoints = 0;
113  m_DiscardedKeyPoints = 0;
114 
115  typename InputImageType::PointType origin0 = input->GetOrigin();
116 
117  ComputeDifferenceOfGaussian(input);
118  DetectKeyPoint(lOctave);
119 
120  // Get the last gaussian for subsample and
121  // repeat the process
122  m_ShrinkFilter = ShrinkFilterType::New();
123  m_ShrinkFilter->SetInput(m_LastGaussian);
124  m_ShrinkFilter->SetShrinkFactors(m_ShrinkFactors);
125  m_ShrinkFilter->Update();
126 
127  input = m_ShrinkFilter->GetOutput();
128 
129  typename InputImageType::PointType origin1;
130  typename InputImageType::SpacingType spacing = input->GetSignedSpacing();
131 
132  origin1[0] = origin0[0] + spacing[0] * 0.25;
133  origin1[1] = origin0[1] + spacing[1] * 0.25;
134 
135  input->SetOrigin(origin1);
136 
137  otbGenericMsgDebugMacro(<< "ImageToSIFTKeyPointSetFilter:: Number total key points : " << m_ValidatedKeyPoints);
138  otbGenericMsgDebugMacro(<< "ImageToSIFTKeyPointSetFilter:: Number different sample key points per octave : " << m_DifferentSamplePoints);
139  otbGenericMsgDebugMacro(<< "ImageToSIFTKeyPointSetFilter:: Number discarded key points per octave : " << m_DiscardedKeyPoints);
140  otbGenericMsgDebugMacro(<< "ImageToSIFTKeyPointSetFilter:: Resample image factor : " << m_ShrinkFactors);
141  }
142 
143  otbGenericMsgDebugMacro(<< "ImageToSIFTKeyPointSetFilter:: Total number key points : " << this->GetOutput()->GetNumberOfPoints());
144 }
145 
149 template <class TInputImage, class TOutputPointSet>
150 void ImageToSIFTKeyPointSetFilter<TInputImage, TOutputPointSet>::InitializeInputImage()
151 {
152  m_ExpandFilter->SetInput(this->GetInput());
153  m_ExpandFilter->SetExpandFactors(m_ExpandFactors);
154  m_ExpandFilter->Update();
156 
157  typename InputImageType::PointType origin0 = this->GetInput()->GetOrigin();
158  typename InputImageType::PointType origin1;
159  typename InputImageType::SpacingType spacing = m_ExpandFilter->GetOutput()->GetSignedSpacing();
160 
161  origin1[0] = origin0[0] - spacing[0] * 0.5;
162  origin1[1] = origin0[1] - spacing[1] * 0.5;
163 
164  m_ExpandFilter->GetOutput()->SetOrigin(origin1);
165 }
166 
170 template <class TInputImage, class TOutputPointSet>
171 void ImageToSIFTKeyPointSetFilter<TInputImage, TOutputPointSet>::ComputeDifferenceOfGaussian(InputImagePointerType input)
172 {
173  unsigned int lScale = 0;
174  InputImagePointerType previousGaussian;
175 
176  m_DoGList = ImageListType::New();
177 
178  m_MagnitudeList = ImageListType::New();
179  m_OrientationList = ImageListType::New();
180  // m_GaussianWeightOrientationList = ImageListType::New();
181  // m_GaussianWeightDescriptorList = ImageListType::New();
182 
183  // itkRecursiveGaussian use spacing to compute
184  // length sigma gaussian (in mm)
185  // sigma = sigma/spacing
186  //
187  // with multiply by spacing before filtering, length sigma gaussian
188  // is compute in pixel
189  double xsigman = std::abs(input->GetSignedSpacing()[0]) * m_Sigma0;
190  double ysigman = std::abs(input->GetSignedSpacing()[1]) * m_Sigma0;
191 
192  for (lScale = 0; lScale != m_ScalesNumber + 2; lScale++)
193  {
194  m_XGaussianFilter = GaussianFilterType::New();
195  m_YGaussianFilter = GaussianFilterType::New();
196 
197  m_XGaussianFilter->SetSigma(xsigman);
198  m_XGaussianFilter->SetDirection(0);
199  m_XGaussianFilter->SetInput(input);
200 
201  m_YGaussianFilter->SetSigma(ysigman);
202  m_YGaussianFilter->SetDirection(1);
203  m_YGaussianFilter->SetInput(m_XGaussianFilter->GetOutput());
204 
205  m_YGaussianFilter->Update();
206 
207  m_GradientFilter = GradientFilterType::New();
208  m_MagnitudeFilter = MagnitudeFilterType::New();
209  m_OrientationFilter = OrientationFilterType::New();
210 
211  m_GradientFilter->SetInput(m_YGaussianFilter->GetOutput());
212  m_MagnitudeFilter->SetInput(m_GradientFilter->GetOutput());
213  m_OrientationFilter->SetInput(m_GradientFilter->GetOutput());
214 
215  m_MagnitudeFilter->Update();
216  m_OrientationFilter->Update();
217 
218  m_MagnitudeList->PushBack(m_MagnitudeFilter->GetOutput());
219  m_OrientationList->PushBack(m_OrientationFilter->GetOutput());
220 
221  if (lScale > 0)
222  {
223  m_SubtractFilter = SubtractFilterType::New();
224  m_SubtractFilter->SetInput1(m_YGaussianFilter->GetOutput());
225  m_SubtractFilter->SetInput2(previousGaussian);
226  m_SubtractFilter->Update();
227  m_DoGList->PushBack(m_SubtractFilter->GetOutput());
228  }
229 
230  previousGaussian = m_YGaussianFilter->GetOutput();
231  xsigman = xsigman * m_Sigmak;
232  ysigman = ysigman * m_Sigmak;
233  }
234  m_LastGaussian = previousGaussian;
235  otbGenericMsgDebugMacro(<< "ImageToSIFTKeyPointSetFilter:: Number of DoG " << m_DoGList->Size());
236 }
237 
241 template <class TInputImage, class TOutputPointSet>
242 void ImageToSIFTKeyPointSetFilter<TInputImage, TOutputPointSet>::DetectKeyPoint(const unsigned int octave)
243 {
244  // need at least 3 DoG, ie 2 scales
245  if (m_ScalesNumber > 1)
246  {
247  typename ImageListType::Iterator lIterDoG = m_DoGList->Begin() + 1;
248  unsigned int lScale = 1;
249  OutputPointSetPointerType outputPointSet = this->GetOutput();
250  typename InputImageType::SpacingType spacing = lIterDoG.Get()->GetSignedSpacing();
252 
253  while ((lIterDoG + 1) != m_DoGList->End())
254  {
255  otbGenericMsgDebugMacro(<< "ImageToSIFTKeyPointSetFilter:: octave: " << octave << " scale: " << lScale);
256  otbUnusedMacro(octave);
257  // Compute max of DoG
258  MinimumMaximumCalculatorPointerType lMaximumCalculator = MinimumMaximumCalculatorType::New();
259  lMaximumCalculator->SetImage(lIterDoG.Get());
260  lMaximumCalculator->Compute();
261 
262  typename InputImageType::SizeType lRadius;
263  lRadius.Fill(1);
264  typename ImageListType::Iterator lIterNext = lIterDoG + 1;
265  typename ImageListType::Iterator lIterPrev = lIterDoG - 1;
266 
267  NeighborhoodIteratorType lIterCurrent(lRadius, lIterDoG.Get(), lIterDoG.Get()->GetLargestPossibleRegion());
268  NeighborhoodIteratorType lIterLowerAdjacent(lRadius, lIterPrev.Get(), lIterPrev.Get()->GetLargestPossibleRegion());
269  NeighborhoodIteratorType lIterUpperAdjacent(lRadius, lIterNext.Get(), lIterNext.Get()->GetLargestPossibleRegion());
270 
271  while (!lIterCurrent.IsAtEnd() && !lIterLowerAdjacent.IsAtEnd() && !lIterUpperAdjacent.IsAtEnd())
272  {
273  // check local min/max
274  if (IsLocalExtremum(lIterCurrent, lIterLowerAdjacent, lIterUpperAdjacent))
275  {
276  VectorPointType lTranslation(PixelType(0));
277  OffsetType lOffsetZero = {{0, 0}};
278 
279  unsigned int lChangeSamplePoints = 0;
280  NeighborhoodIteratorType neighborCurrentScale(lIterCurrent);
281  NeighborhoodIteratorType neighborPreviousScale(lIterLowerAdjacent);
282  NeighborhoodIteratorType neighborNextScale(lIterUpperAdjacent);
283 
284  bool accepted = false;
285  bool changed = true;
286  while (lChangeSamplePoints < m_ChangeSamplePointsMax && changed)
287  {
288  accepted = RefineLocationKeyPoint(neighborCurrentScale, neighborPreviousScale, neighborNextScale, lTranslation);
289 
290  OffsetType lTranslateOffset = {{0, 0}};
291 
292  lTranslateOffset[0] += static_cast<int>(lTranslation[0] > 0.5);
293  lTranslateOffset[0] += -static_cast<int>(lTranslation[0] < -0.5);
294 
295  lTranslateOffset[1] += static_cast<int>(lTranslation[1] > 0.5);
296  lTranslateOffset[1] += -static_cast<int>(lTranslation[1] < -0.5);
297 
298  NeighborhoodIteratorType moveIterator = neighborCurrentScale + lTranslateOffset;
299 
300  if (moveIterator.InBounds())
301  {
302  changed = lTranslateOffset != lOffsetZero;
303 
304  // move iterator
305  neighborCurrentScale += lTranslateOffset;
306  neighborPreviousScale += lTranslateOffset;
307  neighborNextScale += lTranslateOffset;
308  }
309  else
310  {
311  changed = false;
312  }
313  ++lChangeSamplePoints;
314  }
315  if (changed)
316  {
317  ++m_DifferentSamplePoints;
318  }
319 
320  // add key point
321  if (accepted)
322  {
323  std::vector<PixelType> lOrientations = ComputeKeyPointOrientations(neighborCurrentScale, lScale, lTranslation[2]);
324 
325  // for each main orientation
326  for (typename std::vector<PixelType>::iterator orientationIt = lOrientations.begin(); orientationIt != lOrientations.end(); ++orientationIt)
327  {
328 
329  std::vector<PixelType> lDescriptors = ComputeKeyPointDescriptor(neighborCurrentScale, lScale, *orientationIt);
330 
331  OutputPointType keyPoint;
332 
333  lIterDoG.Get()->TransformIndexToPhysicalPoint(neighborCurrentScale.GetIndex(), keyPoint);
334  keyPoint[0] += spacing[0] * lTranslation[0];
335  keyPoint[1] += spacing[1] * lTranslation[1];
336 
337  outputPointSet->SetPoint(m_ValidatedKeyPoints, keyPoint);
338 
339  OutputPixelType data;
340  data.SetSize(128);
341  // check this, compute scale
342  // real scale = octave*scale
343  typename std::vector<PixelType>::const_iterator lIterDescriptor = lDescriptors.begin();
344 
345  unsigned int lIndDesc = 0;
346  while (lIterDescriptor != lDescriptors.end())
347  {
348  data.SetElement(lIndDesc, *lIterDescriptor);
349  ++lIndDesc;
350  ++lIterDescriptor;
351  }
352  outputPointSet->SetPointData(m_ValidatedKeyPoints, data);
353 
354  ++m_ValidatedKeyPoints;
355  }
356  }
357  }
358 
359  ++lIterCurrent;
360  ++lIterLowerAdjacent;
361  ++lIterUpperAdjacent;
362  }
363 
364  ++lIterDoG;
365  ++lScale;
366  }
367  }
368 }
369 
373 template <class TInputImage, class TOutputPointSet>
374 bool ImageToSIFTKeyPointSetFilter<TInputImage, TOutputPointSet>::IsLocalExtremum(const NeighborhoodIteratorType& currentScale,
375  const NeighborhoodIteratorType& previousScale,
376  const NeighborhoodIteratorType& nextScale) const
377 {
378  bool isMin = currentScale.GetCenterPixel() < currentScale.GetPixel(m_Offsets[0]);
379  bool isMax = currentScale.GetCenterPixel() > currentScale.GetPixel(m_Offsets[0]);
380  bool isExtremum = isMin || isMax;
381  unsigned int lIterOffset = 0;
383 
384  while (isExtremum && lIterOffset != 8)
385  {
386  OffsetType off = m_Offsets[lIterOffset];
387  if (isMin)
388  {
389  isExtremum = currentScale.GetCenterPixel() < currentScale.GetPixel(off) && currentScale.GetCenterPixel() < previousScale.GetPixel(off) &&
390  currentScale.GetCenterPixel() < nextScale.GetPixel(off);
391  }
392  else if (isMax)
393  {
394  isExtremum = currentScale.GetCenterPixel() > currentScale.GetPixel(off) && currentScale.GetCenterPixel() > previousScale.GetPixel(off) &&
395  currentScale.GetCenterPixel() > nextScale.GetPixel(off);
396  }
397  lIterOffset++;
398  }
399  if (isExtremum && isMin)
400  {
401  isExtremum = currentScale.GetCenterPixel() < previousScale.GetCenterPixel() && currentScale.GetCenterPixel() < nextScale.GetCenterPixel();
402  }
403  else if (isExtremum && isMax)
404  {
405  isExtremum = currentScale.GetCenterPixel() > previousScale.GetCenterPixel() && currentScale.GetCenterPixel() > nextScale.GetCenterPixel();
406  }
407  return isExtremum;
408 }
409 
413 template <class TInputImage, class TOutputPointSet>
414 bool ImageToSIFTKeyPointSetFilter<TInputImage, TOutputPointSet>::RefineLocationKeyPoint(const NeighborhoodIteratorType& currentScale,
415  const NeighborhoodIteratorType& previousScale,
416  const NeighborhoodIteratorType& nextScale, VectorPointType& solution)
417 {
418  bool accepted = true;
419  solution = VectorPointType(PixelType(0));
421 
422  PixelType dx = 0.5 * (currentScale.GetPixel(m_Offsets[6]) - currentScale.GetPixel(m_Offsets[1]));
423 
424  PixelType dy = 0.5 * (currentScale.GetPixel(m_Offsets[4]) - currentScale.GetPixel(m_Offsets[3]));
425 
426  PixelType ds = 0.5 * (nextScale.GetCenterPixel() - previousScale.GetCenterPixel());
427 
428  PixelType dxx = currentScale.GetPixel(m_Offsets[6]) - 2 * currentScale.GetCenterPixel() + currentScale.GetPixel(m_Offsets[1]);
429 
430  PixelType dyy = currentScale.GetPixel(m_Offsets[3]) - 2 * currentScale.GetCenterPixel() + currentScale.GetPixel(m_Offsets[4]);
431 
432  PixelType dss = previousScale.GetCenterPixel() - 2 * currentScale.GetCenterPixel() + nextScale.GetCenterPixel();
433 
434  PixelType dxy = 0.25 * (currentScale.GetPixel(m_Offsets[7]) + currentScale.GetPixel(m_Offsets[0]) - currentScale.GetPixel(m_Offsets[2]) -
435  currentScale.GetPixel(m_Offsets[5]));
436 
437  PixelType dxs = 0.25 * (nextScale.GetPixel(m_Offsets[6]) + previousScale.GetPixel(m_Offsets[1]) - nextScale.GetPixel(m_Offsets[1]) -
438  previousScale.GetPixel(m_Offsets[6]));
439 
440  PixelType dys = 0.25 * (nextScale.GetPixel(m_Offsets[4]) + previousScale.GetPixel(m_Offsets[3]) - nextScale.GetPixel(m_Offsets[3]) -
441  previousScale.GetPixel(m_Offsets[4]));
442 
443  // Compute matrice determinant
444  double det = dxx * (dyy * dss - dys * dys) - dxy * (dxy * dss - dxs * dys) + dxs * (dxy * dys - dxs * dyy);
445 
446  // Solve system, compute key point offset
447  solution[0] = -dx * (dyy * dss - dys * dys) - dy * (dxs * dys - dxy * dss) - ds * (dxy * dys - dyy * dxs);
448  solution[1] = -dx * (dys * dxs - dss * dxy) - dy * (dxx * dss - dxs * dxs) - ds * (dxs * dxy - dxx * dys);
449  solution[2] = -dx * (dxy * dys - dxs * dyy) - dy * (dxy * dxs - dxx * dys) - ds * (dxx * dyy - dxy * dxy);
450 
451  // Compute interpolated value DoG for lSolution (determinant factor)
452  PixelType lDoGInterpolated = det * currentScale.GetCenterPixel() + 0.5 * (dx * solution[0] + dy * solution[1] + ds * solution[2]);
453 
454  PixelType lHessianTrace2 = (dxx + dyy) * (dxx + dyy);
455  PixelType lHessianDet = dxx * dyy - dxy * dxy;
456  // DoG threshold
457 
458  accepted = fabs(lDoGInterpolated) >= fabs(det * m_DoGThreshold);
459 
460  // Eliminating edge response
461 
462  accepted = accepted && fabs(lHessianTrace2) < fabs(m_RatioEdgeThreshold * lHessianDet);
463 
464  if (!accepted)
465  {
466  ++m_DiscardedKeyPoints;
467  }
468  if (det < 1e-10f)
469  {
470  solution.Fill(0);
471  }
472  else
473  {
474  // normalize offset with determinant of derivative matrix
475  solution /= det;
476  }
477  return accepted;
478 }
479 
483 template <class TInputImage, class TOutputPointSet>
484 std::vector<typename ImageToSIFTKeyPointSetFilter<TInputImage, TOutputPointSet>::PixelType>
485 ImageToSIFTKeyPointSetFilter<TInputImage, TOutputPointSet>::ComputeKeyPointOrientations(const NeighborhoodIteratorType& currentScale, const unsigned int scale,
486  const PixelType itkNotUsed(translation))
487 {
488  // radius of the neighborhood
489  unsigned int radius = 4;
490  double lSigma = scale * 3;
491  unsigned int nbBins = 36;
492  double binWidth = 360. / nbBins;
493 
494  // initialize the histogram
495  std::vector<double> lHistogram(nbBins, 0.), lSmoothedHistogram(nbBins, 0.);
496 
497  // Build the region to examine
498  typename InputImageType::RegionType region;
499  typename InputImageType::RegionType::SizeType regionSize;
500  typename InputImageType::RegionType::IndexType regionIndex;
501  regionSize.Fill(2 * radius + 2);
502  region.SetSize(regionSize);
503  regionIndex[0] = currentScale.GetIndex()[0] - regionSize[0] / 2;
504  regionIndex[1] = currentScale.GetIndex()[1] - regionSize[1] / 2;
505  region.SetIndex(regionIndex);
506 
507  if (!region.Crop(m_OrientationList->GetNthElement(scale)->GetLargestPossibleRegion()))
508  {
509  itkExceptionMacro(<< "Region " << region << " is strictly outside the largest possible region!");
510  }
511 
512  // iterators on the orientation and the magnitude
513  RegionIteratorType lIterOrientation(m_OrientationList->GetNthElement(scale), region);
514  RegionIteratorType lIterMagn(m_MagnitudeList->GetNthElement(scale), region);
515  lIterOrientation.GoToBegin();
516  lIterMagn.GoToBegin();
517 
518  // For each pixel
519  while (!lIterOrientation.IsAtEnd() && !lIterMagn.IsAtEnd())
520  {
521 
522  // check if pixel is inside the circle of radius
523  float dx = lIterMagn.GetIndex()[0] - currentScale.GetIndex()[0];
524  float dy = lIterMagn.GetIndex()[1] - currentScale.GetIndex()[1];
525  float dist = std::sqrt(dx * dx + dy * dy);
526 
527  // If we are in the circle
528  if (dist < radius)
529  {
530  // Get the values
531  PixelType lOrientation = lIterOrientation.Get();
532  PixelType lMagnitude = lIterMagn.Get();
533 
534  // Compute the gaussian weight
535  double lWeightMagnitude = std::exp(-dist * dist / (2 * lSigma * lSigma));
536 
537  // Compute the histogram bin index
538  unsigned int lHistoIndex = static_cast<unsigned int>(std::floor(nbBins * lOrientation / (CONST_2PI)));
539 
540  // Update the histogram value
541  lHistogram[lHistoIndex] += lMagnitude * lWeightMagnitude;
542  }
543  ++lIterOrientation;
544  ++lIterMagn;
545  }
546 
547  // Computing smoothed histogram and looking for the maximum and a second maximum within 80% of the first
548  double max = 0;
549  double secondMax = 0;
550  double sum = 0;
551  int maxIndex = 0;
552  // int secondMaxIndex = -1;
553  int j = 0;
554  int i = 0;
555 
556  // Smoothing histogram
557  for (i = 0; i < static_cast<int>(nbBins); ++i)
558  {
559  sum = 0;
560  for (j = i - nbBins; j < i; ++j)
561  {
562  sum += lHistogram[i - j - 1] * m_HistogramGaussianWeights[j + nbBins];
563  }
564  lSmoothedHistogram[i] = sum;
565  }
566 
567  // looking for maximums
568  for (i = 0; i < static_cast<int>(nbBins); ++i)
569  {
570  if (lSmoothedHistogram[i] > max)
571  {
572  secondMax = max;
573  // secondMaxIndex = maxIndex;
574  max = lSmoothedHistogram[i];
575  maxIndex = i;
576  }
577  else if (sum > secondMax)
578  {
579  secondMax = lSmoothedHistogram[i];
580  // secondMaxIndex = i;
581  }
582  }
583  // This structure will hold the located maximums
584  std::vector<PixelType> orientations;
585 
586  // interpolate orientation maximum
587  double x1, x2, x3, y1, y2, y3, a, b, num, denom, orientation;
588  x1 = (maxIndex - 1) * binWidth + binWidth / 2;
589  y1 = lSmoothedHistogram[(maxIndex - 1) < 0 ? maxIndex - 1 + nbBins : maxIndex - 1];
590  x2 = (maxIndex)*binWidth + binWidth / 2;
591  y2 = lSmoothedHistogram[maxIndex];
592  x3 = (maxIndex + 1) * binWidth + binWidth / 2;
593  y3 = lSmoothedHistogram[maxIndex + 1 > static_cast<int>(nbBins) - 1 ? maxIndex + 1 - nbBins : maxIndex + 1];
594 
595  denom = x1 * x1 * x2 + x2 * x2 * x3 + x3 * x3 * x1 - x1 * x1 * x3 - x2 * x2 * x1 - x3 * x3 * x2;
596  num = y1 * x2 + y2 * x3 + y3 * x1 - y1 * x3 - y2 * x1 - y3 * x2;
597 
598  if (denom == 0 || num == 0)
599  {
600  // no main orientation, return an empty orientation vector
601  return orientations;
602  }
603 
604  a = num / denom;
605  b = ((y1 - y2) - a * (x1 * x1 - x2 * x2)) / (x1 - x2);
606 
607  orientation = -b / (2 * a);
608  if (orientation < 0)
609  {
610  orientation += 360;
611  }
612  else if (orientation >= 360)
613  {
614  orientation -= 360;
615  }
616 
617  // orientations.push_back( static_cast<PixelType>(maxIndex*binWidth + binWidth/2));
618  orientations.push_back(static_cast<PixelType>(orientation));
619 
620  // Second peak is disabled, since it seems to confuse the matching procedure.
621 
622  // if(secondMaxIndex>=0 && secondMax > 0.8 * max)
623  // {
624  // x1 = (secondMaxIndex-1)*binWidth+binWidth/2;
625  // y1 = lSmoothedHistogram[(secondMaxIndex-1)<0 ? secondMaxIndex-1+nbBins : secondMaxIndex-1];
626  // x2 = (secondMaxIndex)*binWidth+binWidth/2;
627  // y2 = lSmoothedHistogram[secondMaxIndex];
628  // x3 = (secondMaxIndex+1)*binWidth+binWidth/2;
629  // y3 = lSmoothedHistogram[secondMaxIndex+1>static_cast<int>(nbBins)-1 ? secondMaxIndex+1-nbBins : secondMaxIndex+1];
630 
631  // denom = x1*x1*x2 + x2*x2*x3 + x3*x3*x1 - x1*x1*x3 - x2*x2*x1 - x3*x3*x2;
632  // num = y1*x2 + y2 * x3 + y3*x1 - y1*x3 - y2*x1 - y3*x2;
633 
634  // if(denom == 0 || num == 0)
635  // {
636  // // no main orientation, return an empty orientation vector
637  // return orientations;
638  // }
639 
640  // a = num/denom;
641  // b = ((y1-y2)-a*(x1*x1-x2*x2))/(x1-x2);
642 
643  // orientation = -b/(2*a);
644  // if(orientation<0)
645  // {
646  // orientation+=360;
647  // }
648  // else if(orientation>=360)
649  // {
650  // orientation-=360;
651  // }
652  // // orientations.push_back( static_cast<PixelType>(secondMaxIndex*binWidth+binWidth/2));
653  // orientations.push_back(static_cast<PixelType>(orientation));
654  // }
655 
656  return orientations;
657 }
658 
662 template <class TInputImage, class TOutputPointSet>
663 std::vector<typename ImageToSIFTKeyPointSetFilter<TInputImage, TOutputPointSet>::PixelType>
664 ImageToSIFTKeyPointSetFilter<TInputImage, TOutputPointSet>::ComputeKeyPointDescriptor(const NeighborhoodIteratorType& currentScale, const unsigned int scale,
665  const PixelType& orientation)
666 {
667  std::vector<PixelType> lHistogram(128, 0.);
668 
669  typename InputImageType::RegionType region;
670  typename InputImageType::RegionType::SizeType regionSize;
671  typename InputImageType::RegionType::IndexType regionIndex;
672 
673  unsigned int nbHistograms = 4;
674  unsigned int nbPixelsPerHistogram = 4;
675  unsigned int nbBinsPerHistogram = 8;
676 
677  float radius = static_cast<float>(nbHistograms / 2 * nbPixelsPerHistogram);
678 
679  // std::cout<<"Radius: "<<radius<<std::endl;
680 
681  // 4 region of 4 pixels plus 2 pixels of margin
682  regionSize[0] = nbHistograms * nbPixelsPerHistogram + 2;
683  regionSize[1] = nbHistograms * nbPixelsPerHistogram + 2;
684 
685  // sigma set to one half the width of descriptor window
686  // TODO check this
687  double lSigma = radius;
688 
689  // index - regionSize/2
690  regionIndex[0] = currentScale.GetIndex()[0] - regionSize[0] / 2;
691  regionIndex[1] = currentScale.GetIndex()[1] - regionSize[1] / 2;
692 
693  region.SetIndex(regionIndex);
694  region.SetSize(regionSize);
695 
696  // Crop with largest region
697  if (!region.Crop(m_OrientationList->GetNthElement(scale)->GetLargestPossibleRegion()))
698  {
699  itkExceptionMacro(<< "Region " << region << " is outside of the largest possible region!");
700  }
701  RegionIteratorType lIterMagnitude(m_MagnitudeList->GetNthElement(scale), region);
702  RegionIteratorType lIterOrientation(m_OrientationList->GetNthElement(scale), region);
703  lIterMagnitude.GoToBegin();
704  lIterOrientation.GoToBegin();
705 
706  // For each pixel in the region
707  while (!lIterMagnitude.IsAtEnd() && !lIterOrientation.IsAtEnd())
708  {
709  // check if pixel is inside the circle of radius
710  float dx = lIterMagnitude.GetIndex()[0] - currentScale.GetIndex()[0];
711  float dy = lIterMagnitude.GetIndex()[1] - currentScale.GetIndex()[1];
712  float dist = std::sqrt(dx * dx + dy * dy);
713 
714  // If we are in the circle
715  if (dist < radius)
716  {
717  // rotate the pixel location to compensate sift orientation
718  float angle = orientation * CONST_PI_180;
719  float cosangle = std::cos(-angle);
720  float sinangle = std::sin(-angle);
721  float rdx = dx * cosangle - dy * sinangle;
722  float rdy = dx * sinangle + dy * cosangle;
723  // decide to which histogram the pixel contributes
724  unsigned int xHistogramIndex = static_cast<unsigned int>(std::floor((rdx + radius) / static_cast<float>(nbPixelsPerHistogram)));
725  unsigned int yHistogramIndex = static_cast<unsigned int>(std::floor((rdy + radius) / static_cast<float>(nbPixelsPerHistogram)));
726 
727  // decide to which bin of the histogram the pixel contributes
728  float compensatedOrientation = lIterOrientation.Get() - angle;
729  if (compensatedOrientation < 0)
730  {
731  compensatedOrientation += CONST_2PI;
732  }
733  if (compensatedOrientation >= CONST_2PI)
734  {
735  compensatedOrientation -= CONST_2PI;
736  }
737  unsigned int histogramBin = static_cast<unsigned int>(std::floor(compensatedOrientation * nbBinsPerHistogram / (CONST_2PI)));
738 
739  // Compute the wheight of the pixel in the histogram
740  double lWeightMagnitude = std::exp(-(dist * dist) / (2 * lSigma * lSigma));
741 
742  // Compute the global descriptor index
743  unsigned int descriptorIndex = yHistogramIndex * nbBinsPerHistogram * nbHistograms + xHistogramIndex * nbBinsPerHistogram + histogramBin;
744  lHistogram[descriptorIndex] += lIterMagnitude.Get() * lWeightMagnitude;
745  }
746 
747  ++lIterOrientation;
748  ++lIterMagnitude;
749  }
750 
751  // normalize histogram to unit length
752  typename std::vector<PixelType>::iterator lIterHisto = lHistogram.begin();
753  float lNorm = 0.0;
754 
755  while (lIterHisto != lHistogram.end())
756  {
757  lNorm = lNorm + (*lIterHisto) * (*lIterHisto);
758  ++lIterHisto;
759  }
760  lNorm = std::sqrt(lNorm);
761 
762  lIterHisto = lHistogram.begin();
763  while (lIterHisto != lHistogram.end())
764  {
765  if (lNorm > 0)
766  {
767  *lIterHisto = (*lIterHisto) / lNorm;
768  }
769  else
770  {
771  *lIterHisto = m_GradientMagnitudeThreshold;
772  }
773 
774  // threshold gradient magnitude
775  if (*lIterHisto > m_GradientMagnitudeThreshold)
776  {
777  *lIterHisto = m_GradientMagnitudeThreshold;
778  }
779  ++lIterHisto;
780  }
781 
782  // renormalize histogram to unit length
783  lIterHisto = lHistogram.begin();
784  lNorm = 0.0;
785 
786  while (lIterHisto != lHistogram.end())
787  {
788  lNorm = lNorm + (*lIterHisto) * (*lIterHisto);
789  ++lIterHisto;
790  }
791  lNorm = std::sqrt(lNorm);
792 
793  lIterHisto = lHistogram.begin();
794  while (lIterHisto != lHistogram.end())
795  {
796  *lIterHisto = (*lIterHisto) / lNorm;
797  ++lIterHisto;
798  }
799 
800  return lHistogram;
801 }
802 
806 template <class TInputImage, class TOutputPointSet>
807 void ImageToSIFTKeyPointSetFilter<TInputImage, TOutputPointSet>::PrintSelf(std::ostream& os, itk::Indent indent) const
808 {
809  const OutputPointSetType* output = dynamic_cast<const OutputPointSetType*>(this->Superclass::ProcessObjectType::GetOutput(0));
810 
811  Superclass::PrintSelf(os, indent);
812  os << indent << "Number of octaves: " << m_OctavesNumber << std::endl;
813  os << indent << "Number of scales: " << m_ScalesNumber << std::endl;
814 
815  os << indent << "Number of SIFT key points: " << output->GetNumberOfPoints() << std::endl;
816 }
817 
818 } // End namespace otb
819 
820 #endif
otb::ImageToSIFTKeyPointSetFilter::MinimumMaximumCalculatorPointerType
MinimumMaximumCalculatorType::Pointer MinimumMaximumCalculatorPointerType
Definition: otbImageToSIFTKeyPointSetFilter.h:216
otbUnusedMacro
#define otbUnusedMacro(x)
Definition: otbMacro.h:188
otb::ImageToSIFTKeyPointSetFilter::NeighborhoodIteratorType
itk::ConstNeighborhoodIterator< InputImageType > NeighborhoodIteratorType
Definition: otbImageToSIFTKeyPointSetFilter.h:209
otb::ImageToSIFTKeyPointSetFilter::RegionIteratorType
itk::ImageRegionConstIterator< InputImageType > RegionIteratorType
Definition: otbImageToSIFTKeyPointSetFilter.h:213
otb::ImageToPointSetFilter::OutputPointSetType
Superclass::OutputPointSetType OutputPointSetType
Definition: otbImageToPointSetFilter.h:65
otb::ImageToSIFTKeyPointSetFilter::VectorPointType
itk::Vector< PixelType, 3 > VectorPointType
Definition: otbImageToSIFTKeyPointSetFilter.h:142
otb::ImageToSIFTKeyPointSetFilter::OutputPointType
TOutputPointSet::PointType OutputPointType
Definition: otbImageToSIFTKeyPointSetFilter.h:139
otb::ImageToSIFTKeyPointSetFilter::GenerateData
void GenerateData() override
Definition: otbImageToSIFTKeyPointSetFilter.hxx:92
otb::ImageToSIFTKeyPointSetFilter::IsLocalExtremum
bool IsLocalExtremum(const NeighborhoodIteratorType &currentScale, const NeighborhoodIteratorType &previousScale, const NeighborhoodIteratorType &nextScale) const
Definition: otbImageToSIFTKeyPointSetFilter.hxx:374
otb
The "otb" namespace contains all Orfeo Toolbox (OTB) classes.
Definition: otbJoinContainer.h:32
otb::ImageToSIFTKeyPointSetFilter::ExpandFilterType
itk::ExpandImageFilter< TInputImage, TInputImage > ExpandFilterType
Definition: otbImageToSIFTKeyPointSetFilter.h:190
otbGenericMsgDebugMacro
#define otbGenericMsgDebugMacro(x)
Definition: otbMacro.h:63
otb::ImageToSIFTKeyPointSetFilter::ComputeKeyPointOrientations
std::vector< PixelType > ComputeKeyPointOrientations(const NeighborhoodIteratorType &currentScale, const unsigned int scale, const PixelType translation)
Definition: otbImageToSIFTKeyPointSetFilter.hxx:485
otb::ImageToSIFTKeyPointSetFilter::PrintSelf
void PrintSelf(std::ostream &os, itk::Indent indent) const override
Definition: otbImageToSIFTKeyPointSetFilter.hxx:807
otb::ImageToSIFTKeyPointSetFilter::ComputeDifferenceOfGaussian
void ComputeDifferenceOfGaussian(InputImagePointerType input)
Definition: otbImageToSIFTKeyPointSetFilter.hxx:171
otb::ImageToSIFTKeyPointSetFilter::DetectKeyPoint
void DetectKeyPoint(const unsigned int octave)
Definition: otbImageToSIFTKeyPointSetFilter.hxx:242
otb::ImageToSIFTKeyPointSetFilter::OutputPointSetPointerType
TOutputPointSet::Pointer OutputPointSetPointerType
Definition: otbImageToSIFTKeyPointSetFilter.h:137
otb::CONST_PI_180
constexpr double CONST_PI_180
Definition: otbMath.h:56
otb::ImageToSIFTKeyPointSetFilter::InputImagePointerType
TInputImage::Pointer InputImagePointerType
Definition: otbImageToSIFTKeyPointSetFilter.h:133
otb::CONST_2PI
constexpr double CONST_2PI
Definition: otbMath.h:55
otb::ImageToSIFTKeyPointSetFilter::RefineLocationKeyPoint
bool RefineLocationKeyPoint(const NeighborhoodIteratorType &currentScale, const NeighborhoodIteratorType &previousScale, const NeighborhoodIteratorType &nextScale, VectorPointType &solution)
Definition: otbImageToSIFTKeyPointSetFilter.hxx:414
otbImageToSIFTKeyPointSetFilter.h
otb::ImageToSIFTKeyPointSetFilter::InitializeInputImage
void InitializeInputImage()
Definition: otbImageToSIFTKeyPointSetFilter.hxx:150
otb::ImageToSIFTKeyPointSetFilter::OffsetType
NeighborhoodType::OffsetType OffsetType
Definition: otbImageToSIFTKeyPointSetFilter.h:211
otb::ImageToSIFTKeyPointSetFilter::PixelType
TInputImage::PixelType PixelType
Definition: otbImageToSIFTKeyPointSetFilter.h:134
otb::ImageList::Iterator
Superclass::Iterator Iterator
Definition: otbImageList.h:57
otb::ImageToSIFTKeyPointSetFilter::ComputeKeyPointDescriptor
std::vector< PixelType > ComputeKeyPointDescriptor(const NeighborhoodIteratorType &currentScale, const unsigned int scale, const PixelType &orientation)
Definition: otbImageToSIFTKeyPointSetFilter.hxx:664
otb::ImageToSIFTKeyPointSetFilter::ImageToSIFTKeyPointSetFilter
ImageToSIFTKeyPointSetFilter()
Definition: otbImageToSIFTKeyPointSetFilter.hxx:36
otb::ImageToSIFTKeyPointSetFilter::OutputPixelType
TOutputPointSet::PixelType OutputPixelType
Definition: otbImageToSIFTKeyPointSetFilter.h:138
Generated at Wed Feb 28 2024 15:22:37 for OTB with   doxygen 1.8.17