OTB  9.0.0
Orfeo Toolbox
otbImageToSURFKeyPointSetFilter.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 otbImageToSURFKeyPointSetFilter_hxx
22 #define otbImageToSURFKeyPointSetFilter_hxx
23 
24 #include "otbImageToSURFKeyPointSetFilter.h"
25 #include "itkCenteredRigid2DTransform.h"
26 
27 namespace otb
28 {
32 template <class TInputImage, class TOutputPointSet>
33 ImageToSURFKeyPointSetFilter<TInputImage, TOutputPointSet>::ImageToSURFKeyPointSetFilter()
34 {
35  m_Offsets[0][0] = -1;
36  m_Offsets[0][1] = -1;
37  m_Offsets[1][0] = -1;
38  m_Offsets[1][1] = 0;
39  m_Offsets[2][0] = -1;
40  m_Offsets[2][1] = 1;
41  m_Offsets[3][0] = 0;
42  m_Offsets[3][1] = -1;
43  m_Offsets[4][0] = 0;
44  m_Offsets[4][1] = 1;
45  m_Offsets[5][0] = 1;
46  m_Offsets[5][1] = -1;
47  m_Offsets[6][0] = 1;
48  m_Offsets[6][1] = 0;
49  m_Offsets[7][0] = 1;
50  m_Offsets[7][1] = 1;
51 
52  m_OctavesNumber = 1;
53  m_ScalesNumber = 3;
54  m_NumberOfPoints = 0;
55  m_DifferentSamplePoints = 0;
56  m_DoHThreshold = 0.03;
57  m_DetHessianFilter = ImageToDetHessianImageType::New();
58 }
59 
60 /*---------------------------------------------------------
61  * Destructor.c
62  ----------------------------------------------------------*/
63 
64 template <class TInputImage, class TOutputPointSet>
65 ImageToSURFKeyPointSetFilter<TInputImage, TOutputPointSet>::~ImageToSURFKeyPointSetFilter()
66 {
67 }
68 
69 /*-------------------------------------------------------
70  * Generate Data
71  --------------------------------------------------------*/
72 template <class TInputImage, class TOutputPointSet>
73 void ImageToSURFKeyPointSetFilter<TInputImage, TOutputPointSet>::GenerateData(void)
74 {
75  double k;
76  double sigma_in = 2.;
77  SizeType radius;
78  SpacingType spacing;
79  /*Output*/
80  OutputPointSetPointerType outputPointSet = this->GetOutput();
81 
85  if (m_ScalesNumber > 1)
86  k = (double)std::pow(2.0, (double)(1 / (double)(m_ScalesNumber - 1)));
87  else
88  k = 3;
90 
91  /* Computation loop over octaves*/
92  for (int i = 0; i < m_OctavesNumber; ++i)
93  {
94 
95  sigma_in = 2.;
96  m_ImageList = ImageListType::New();
97  spacing = this->GetInput()->GetSignedSpacing();
98 
99  /*--------------------------------------------------------
100  Octave per octave
101  --------------------------------------------------------*/
102  if (i > 0)
103  {
104 
105  m_ResampleFilter = ResampleFilterType::New();
106  m_ResampleFilter->SetInput(this->GetInput());
107 
108  SizeType size = this->GetInput()->GetLargestPossibleRegion().GetSize();
109  for (int l = 0; l < 2; ++l)
110  size[l] = (unsigned int)floor(size[l] / std::pow(2.0, i));
111  m_ResampleFilter->SetSize(size);
112 
113 
114  for (int l = 0; l < 2; ++l)
115  spacing[l] = (spacing[l] * std::pow(2.0, i));
116  m_ResampleFilter->SetOutputSpacing(spacing);
117  /* necessary to handle images where the origin is not (0, 0) */
118  m_ResampleFilter->SetOutputOrigin(this->GetInput()->GetOrigin());
119 
120  m_ResampleFilter->SetDefaultPixelValue(0);
121  m_ResampleFilter->Update();
122  m_determinantImage = m_ResampleFilter->GetOutput();
123 
124  otbGenericMsgDebugMacro(<< "ImageToSURFKeyPointSetFilter:: Size of the image at the octave : " << i << " is "
125  << m_determinantImage->GetLargestPossibleRegion().GetSize());
126  }
127 
128  for (int j = 0; j < m_ScalesNumber; ++j)
129  {
135  if ((i != 0 && j != 0) || (i == 0 && (i + j != 0)) || (m_ScalesNumber == 1 && i != 0))
136  sigma_in *= k;
137 
143  m_DetHessianFilter = ImageToDetHessianImageType::New();
144 
145  if (i == 0)
146  m_DetHessianFilter->SetInput(this->GetInput());
147  else
148  m_DetHessianFilter->SetInput(m_determinantImage);
149 
150  m_DetHessianFilter->SetSigma(sigma_in);
151  m_DetHessianFilter->Update();
152  m_determinantImage = m_DetHessianFilter->GetOutput();
153 
154  if (i + j == 0)
155  {
156  otbGenericMsgDebugMacro(<< "ImageToSURFKeyPointSetFilter:: Size of the image at the octave : " << i << " is "
157  << m_determinantImage->GetLargestPossibleRegion().GetSize());
158  }
159 
161  m_ImageList->PushBack(m_determinantImage);
162  }
163 
164  /*----------------------------------------------------*/
165  /* extremum Search over octave's scales */
166  /*----------------------------------------------------*/
167 
168  for (int jj = 1; jj < (int)(m_ImageList->Size() - 1); jj++)
169  {
170  m_ImageCurrent = m_ImageList->GetNthElement(jj);
171  m_ImageMovedPrev = m_ImageList->GetNthElement(jj - 1);
172  m_ImageMovedNext = m_ImageList->GetNthElement(jj + 1);
173 
175  radius.Fill(1);
176  NeighborhoodIteratorType it(radius, m_ImageCurrent, m_ImageCurrent->GetLargestPossibleRegion());
177  it.GoToBegin();
179 
180  /* NeighborhoodIterator Adjacents parameters*/
181  NeighborhoodIteratorType itNeighPrev(radius, m_ImageMovedPrev, m_ImageMovedPrev->GetLargestPossibleRegion());
182  itNeighPrev.GoToBegin();
183 
184  NeighborhoodIteratorType itNeighNext(radius, m_ImageMovedNext, m_ImageMovedNext->GetLargestPossibleRegion());
185  itNeighNext.GoToBegin();
186 
187  while (!it.IsAtEnd())
188  {
193  if (IsLocalExtremum(it.GetNeighborhood()) && IsLocalExtremumAround(itNeighPrev.GetNeighborhood(), m_ImageCurrent->GetPixel(it.GetIndex())) &&
194  IsLocalExtremumAround(itNeighNext.GetNeighborhood(), m_ImageCurrent->GetPixel(it.GetIndex())))
195  {
196  OutputPointType keyPoint;
197  itNeighPrev.SetLocation(it.GetIndex());
198  itNeighNext.SetLocation(it.GetIndex());
199  VectorPointType lTranslation(PixelValue(0));
201 
202  /* Approximation de la position */
203  NeighborhoodIteratorType neighborCurrentScale(it);
204  NeighborhoodIteratorType neighborPreviousScale(itNeighPrev);
205  NeighborhoodIteratorType neighborNextScale(itNeighNext);
206  bool accepted = false;
207 
208  accepted = RefineLocationKeyPoint(neighborCurrentScale, neighborPreviousScale, neighborNextScale, lTranslation);
209 
210  OffsetType lTranslateOffset = {{0, 0}};
211 
212 
213  lTranslateOffset[0] += static_cast<int>(lTranslation[0] > 0.5);
214  lTranslateOffset[0] += -static_cast<int>(lTranslation[0] < -0.5);
215 
216  lTranslateOffset[1] += static_cast<int>(lTranslation[1] > 0.5);
217  lTranslateOffset[1] += -static_cast<int>(lTranslation[1] < -0.5);
218 
219  NeighborhoodIteratorType moveIterator = neighborCurrentScale + lTranslateOffset;
220 
221  if (moveIterator.InBounds())
222  {
223  // move iterator
224  neighborCurrentScale += lTranslateOffset;
225  neighborPreviousScale += lTranslateOffset;
226  neighborNextScale += lTranslateOffset;
227  lTranslation[0] -= lTranslateOffset[0];
228  lTranslation[1] -= lTranslateOffset[1];
229  }
230  else
231  {
232  lTranslation[0] = 0.0;
233  lTranslation[1] = 0.0;
234  }
235 
236  if (accepted == false)
237  {
238  ++it;
239  ++itNeighPrev;
240  ++itNeighNext;
241  continue;
242  }
243 
244 
245  typename InputImageType::IndexType indexKeyPoint;
246  indexKeyPoint[0] = neighborCurrentScale.GetIndex()[0];
247  indexKeyPoint[1] = neighborCurrentScale.GetIndex()[1];
248 
249  double sigmaDetected = sigma_in / pow(k, (double)jj) * pow(2., (double)i);
250 
251  radius.Fill(GetMin((int)(this->GetInput()->GetLargestPossibleRegion().GetSize()[0] - indexKeyPoint[0]),
252  (int)(this->GetInput()->GetLargestPossibleRegion().GetSize()[1] - indexKeyPoint[1]), (int)(6 * sigmaDetected)));
253 
254  /*
255  Computing the orientation of the key point detected
256  */
257  NeighborhoodIteratorType itNeighOrientation(radius, this->GetInput(), this->GetInput()->GetLargestPossibleRegion());
258 
259  itNeighOrientation.SetLocation(neighborCurrentScale.GetIndex());
260 
261  double orientationDetected = AssignOrientation(itNeighOrientation.GetNeighborhood(), sigmaDetected);
262 
263  /*Filling the Point pointSet Part*/
264  typename InputImageType::PointType physicalKeyPoint;
265  m_ImageCurrent->TransformIndexToPhysicalPoint(neighborCurrentScale.GetIndex(), keyPoint);
266 
267  physicalKeyPoint[0] = keyPoint[0] + spacing[0] * lTranslation[0];
268  physicalKeyPoint[1] = keyPoint[1] + spacing[1] * lTranslation[1];
269 
270  outputPointSet->SetPoint(m_NumberOfPoints, physicalKeyPoint);
271 
272  /*----------------------------------------*/
273  /* Descriptor Computation */
274  /*----------------------------------------*/
275 
276  radius.Fill(GetMin((int)(this->GetInput()->GetLargestPossibleRegion().GetSize()[0] - indexKeyPoint[0]),
277  (int)(this->GetInput()->GetLargestPossibleRegion().GetSize()[1] - indexKeyPoint[1]),
278  (int)(10 * sigmaDetected))); // TODO a changer sigmaDetected par Keypoint[2]
279 
280  NeighborhoodIteratorType itNeighDescriptor(radius, this->GetInput(), this->GetInput()->GetLargestPossibleRegion());
281  // itNeighDescriptor.SetLocation(it.GetIndex());
282  itNeighDescriptor.SetLocation(neighborCurrentScale.GetIndex());
283  VectorType descriptor;
284  descriptor.resize(64);
285  // descriptor = ComputeDescriptor(itNeighDescriptor.GetNeighborhood(), keyPoint[3], keyPoint[2]);
286  descriptor = ComputeDescriptor(itNeighDescriptor.GetNeighborhood(), orientationDetected, sigmaDetected);
287 
288  /*Updating the pointset data with values of the descriptor*/
289  OutputPixelType data;
290  data.SetSize(64);
291 
292  unsigned int IndDescriptor = 0;
293 
294  typename std::vector<double>::const_iterator itDescriptor = descriptor.begin();
295 
296  while (itDescriptor != descriptor.end())
297  {
298  data.SetElement(IndDescriptor, *itDescriptor);
299  IndDescriptor++;
300  itDescriptor++;
301  }
302  outputPointSet->SetPointData(m_NumberOfPoints, data);
303 
304  m_NumberOfPoints++;
305  }
306  ++it;
307  ++itNeighPrev;
308  ++itNeighNext;
309 
310  } /*End while for extremum search*/
311 
312  } /*End Iteration over scales */
313 
314  m_ImageList->Clear();
315 
316  } /* End Key Points search*/
317 
318  otbGenericMsgDebugMacro(<< "ImageToSURFKeyPointSetFilter:: Total Number of Key points " << m_NumberOfPoints);
319 
322 template <class TInputImage, class TOutputPointSet>
323 int ImageToSURFKeyPointSetFilter<TInputImage, TOutputPointSet>::GetMin(int a, int b, int c)
324 {
325  return std::min(a, std::min(b, c));
326 }
327 
328 /*-----------------------------------------------------------
329  * Find Local Extremum
330  -----------------------------------------------------------*/
331 template <class TInputImage, class TOutputPointSet>
332 bool ImageToSURFKeyPointSetFilter<TInputImage, TOutputPointSet>::IsLocalExtremum(const NeighborhoodType& neigh)
333 {
334  int centerIndex = neigh.GetCenterNeighborhoodIndex(), i = 0;
335  double centerValue = neigh[centerIndex];
336  bool max = false, min = false;
337  int flag_min = 0, flag_max = 0;
338 
339  while (i != (int)neigh.Size())
340  {
341  if (i != centerIndex)
342  {
343  if (centerValue > neigh[i] && flag_max == 0)
344  max = true;
345  else
346  {
347  max = false;
348  flag_max = 1;
349  }
350 
351  if (centerValue < neigh[i] && flag_min == 0 && centerValue < 0)
352  min = true;
353  else
354  {
355  min = false;
356  flag_min = 1;
357  }
358  }
359  ++i;
360  }
361 
362  return max || min;
363 }
364 
365 /*-----------------------------------------------------------
366  *Find Local Extremum Around
367  -----------------------------------------------------------*/
368 template <class TInputImage, class TOutputPointSet>
369 bool ImageToSURFKeyPointSetFilter<TInputImage, TOutputPointSet>::IsLocalExtremumAround(const NeighborhoodType& neigh, double CenterValue)
370 {
371 
372  int i = 0;
373  bool max = false, min = false;
374  int flag_min = 0, flag_max = 0;
375 
376  while (i != (int)neigh.Size())
377  {
378 
379  if (CenterValue > neigh[i] && flag_max == 0)
380  max = true;
381  else
382  {
383  max = false;
384  flag_max = 1;
385  }
386 
387  if (CenterValue < neigh[i] && flag_min == 0)
388  min = true;
389  else
390  {
391  min = false;
392  flag_min = 1;
393  }
394 
395  ++i;
396  }
397 
398  return max || min;
399 }
400 
405 template <class TInputImage, class TOutputPointSet>
406 bool ImageToSURFKeyPointSetFilter<TInputImage, TOutputPointSet>::RefineLocationKeyPoint(const NeighborhoodIteratorType& currentScale,
407  const NeighborhoodIteratorType& previousScale,
408  const NeighborhoodIteratorType& nextScale, VectorPointType& solution)
409 {
410  bool accepted = true;
411  solution = VectorPointType(PixelValue(0));
413 
414  PixelValue dx = 0.5 * (currentScale.GetPixel(m_Offsets[6]) - currentScale.GetPixel(m_Offsets[1]));
415 
416  PixelValue dy = 0.5 * (currentScale.GetPixel(m_Offsets[4]) - currentScale.GetPixel(m_Offsets[3]));
417 
418  PixelValue ds = 0.5 * (nextScale.GetCenterPixel() - previousScale.GetCenterPixel());
419 
420  PixelValue dxx = currentScale.GetPixel(m_Offsets[6]) - 2 * currentScale.GetCenterPixel() + currentScale.GetPixel(m_Offsets[1]);
421 
422  PixelValue dyy = currentScale.GetPixel(m_Offsets[3]) - 2 * currentScale.GetCenterPixel() + currentScale.GetPixel(m_Offsets[4]);
423 
424  PixelValue dss = previousScale.GetCenterPixel() - 2 * currentScale.GetCenterPixel() + nextScale.GetCenterPixel();
425 
426  PixelValue dxy = 0.25 * (currentScale.GetPixel(m_Offsets[7]) + currentScale.GetPixel(m_Offsets[0]) - currentScale.GetPixel(m_Offsets[2]) -
427  currentScale.GetPixel(m_Offsets[5]));
428 
429  PixelValue dxs = 0.25 * (nextScale.GetPixel(m_Offsets[6]) + previousScale.GetPixel(m_Offsets[1]) - nextScale.GetPixel(m_Offsets[1]) -
430  previousScale.GetPixel(m_Offsets[6]));
431 
432  PixelValue dys = 0.25 * (nextScale.GetPixel(m_Offsets[4]) + previousScale.GetPixel(m_Offsets[3]) - nextScale.GetPixel(m_Offsets[3]) -
433  previousScale.GetPixel(m_Offsets[4]));
434 
435  // We look for the "false" extrema : min-?-max and max-?-min
436  // Those are not compatible with the 3D quadric interpolation
437  PixelValue curr = currentScale.GetCenterPixel();
438  PixelValue prev = previousScale.GetCenterPixel();
439  PixelValue next = nextScale.GetCenterPixel();
440  if ((prev < curr && curr < next) || (next < curr && curr < prev))
441  {
442  // So we use only 2D interpolation in the current scale
443  dss = 1.0;
444  dxs = 0.0;
445  dys = 0.0;
446  ds = 0.0;
447  }
448 
449  // Compute matrice determinant
450  double det = dxx * (dyy * dss - dys * dys) - dxy * (dxy * dss - dxs * dys) + dxs * (dxy * dys - dxs * dyy);
451 
452  // Solve system, compute key point offset
453  solution[0] = -dx * (dyy * dss - dys * dys) - dy * (dxs * dys - dxy * dss) - ds * (dxy * dys - dyy * dxs);
454  solution[1] = -dx * (dys * dxs - dss * dxy) - dy * (dxx * dss - dxs * dxs) - ds * (dxs * dxy - dxx * dys);
455  solution[2] = -dx * (dxy * dys - dxs * dyy) - dy * (dxy * dxs - dxx * dys) - ds * (dxx * dyy - dxy * dxy);
456 
457  // Compute interpolated value Hessian Determinant for lSolution (determinant factor)
458  PixelValue lDoHInterpolated = det * currentScale.GetCenterPixel() + 0.5 * (dx * solution[0] + dy * solution[1] + ds * solution[2]);
459 
460  // DoG threshold : ignore detected extrema if is not sufficiently noticeable
461  accepted = std::abs(lDoHInterpolated) > std::abs(det * m_DoHThreshold);
462 
463  if (!accepted)
464  {
465  //++m_DiscardedKeyPoints; /* not used at the moment */
466  }
467  if (det < 1.0E-10) // this test cancels the shift for every "weak" extrema
468  {
469  solution.Fill(0);
470  }
471  else
472  {
473  // normalize offset with determinant of derivative matrix
474  solution /= det;
475  }
476  /* check that the translation found is inside a cubic pixel */
477  if (std::abs(static_cast<double>(solution[0])) > 1.0 || std::abs(static_cast<double>(solution[1])) > 1.0 || std::abs(static_cast<double>(solution[2])) > 1.0)
478  {
479  accepted = false;
480  // solution.Fill(0);
481  }
482 
483  return accepted;
484 }
485 
486 /*-----------------------------------------------------------
487  * Compute the orientation of The KeyPoint
488  -----------------------------------------------------------*/
489 
490 template <class TInputImage, class TOutputPointSet>
491 double ImageToSURFKeyPointSetFilter<TInputImage, TOutputPointSet>::AssignOrientation(const NeighborhoodType& neigh, double S)
492 {
493 
494  int i = 0;
495  int pas = ((i + S) - (int)(i + S) > 0.5) ? ((int)S + 1) : (int)S;
496  int Largeur = 2 * neigh.GetRadius()[0] + 1; // Width & length of a neighborhood
497  int rayon = neigh.GetRadius()[0]; // radius of the neigh
498  int col, raw;
499  double dist;
500  double w; // weight of the circular gaussian
501 
502  OutputPointType pt;
503 
504  // Gradient orientation histogram
505  double angle;
506  int bin = 0;
507  int Pi = 180;
508  int LengthBin = 60;
509  int NbBins = (2 * Pi / LengthBin);
510  std::vector<double> tab(NbBins * 2, 0.);
511 
512  while (i < (int)neigh.Size())
513  {
514  col = i % Largeur - rayon;
515  raw = i / Largeur - rayon;
516  dist = std::sqrt(static_cast<double>(col * col + raw * raw));
517  col += rayon;
518  raw += rayon; // Backup to the image coordinate axes
519 
520  if (dist < 6 * S)
521  {
522  // Haar Wavelets responses accumulated in an histogram with Pi/3 precision
523  if ((col > pas && col < Largeur - pas) && (raw > pas && raw < Largeur - pas))
524  {
525 
526  w = std::exp(-((col - rayon) * (col - rayon) + (raw - rayon) * (raw - rayon)) / (2 * 2.5 * 2.5 * S * S));
527  pt[0] = (neigh[(col + pas) + raw * Largeur] - neigh[(col - pas) + raw * Largeur]) * w;
528  pt[1] = (neigh[col + (raw + pas) * Largeur] - neigh[col + (raw - pas) * Largeur]) * w;
529 
530  if (pt[0] + pt[1] != 0)
531  {
532  angle = atan(pt[0] / pt[1]) * (Pi / otb::CONST_PI);
533  if (angle < 0)
534  angle += 2 * Pi;
535 
536  bin = (int)(angle / LengthBin);
537 
538  if (bin <= NbBins - 1 || bin >= 0)
539  {
540  tab[2 * bin] += pt[0];
541  tab[2 * bin + 1] += pt[1];
542  }
543  }
544  }
545  }
546  i += pas;
547  }
548 
549  // Find Orientation
550  double indice = 0;
551  double max = 0;
552  double length = 0;
553 
554  // Detecting current point orientation
555  for (int ii = 0; ii < NbBins * 2; ii = ii + 2)
556  {
557  length = std::sqrt(tab[ii] * tab[ii] + tab[ii + 1] * tab[ii + 1]);
558  if (length > max)
559  {
560  max = length;
561  indice = ii / 2;
562  }
563  }
564 
565  return (indice + 0.5) * LengthBin;
566 }
567 
568 /*-----------------------------------------------------------
569  * Compute the descriptor of The KeyPoint
570  -----------------------------------------------------------*/
571 
572 template <class TInputImage, class TOutputPointSet>
573 typename ImageToSURFKeyPointSetFilter<TInputImage, TOutputPointSet>::VectorType
574 ImageToSURFKeyPointSetFilter<TInputImage, TOutputPointSet>::ComputeDescriptor(const NeighborhoodType& neigh, double O, double S)
575 {
576 
577  typedef itk::CenteredRigid2DTransform<double> TransformType;
578  TransformType::Pointer eulerTransform = TransformType::New();
579  TransformType::ParametersType ParamVec(5);
580  PointImageType pSrc, pDst;
581  double angle = O * otb::CONST_PI / 180;
582 
583  int i = 0, col, raw, Nbin, pas = 1;
584  double xx = 0, yy = 0;
585  double dx, dy, w;
586  int Largeur = 2 * neigh.GetRadius()[0] + 1;
587  double rayon = static_cast<double>(Largeur) / 4.;
588  double r = neigh.GetRadius()[0];
589  double dist = 0;
590  double x0 = neigh.GetCenterNeighborhoodIndex() % Largeur;
591  double y0 = neigh.GetCenterNeighborhoodIndex() / Largeur;
592 
593  // std::cout << " x0 " << x0 << " y0 " << y0 << angle << std::endl;
594 
595  VectorType descriptorVector;
596  descriptorVector.resize(64);
597 
599  ParamVec[0] = angle;
600  ParamVec[1] = x0;
601  ParamVec[2] = y0;
602  ParamVec[3] = 0;
603  ParamVec[4] = 0;
604  eulerTransform->SetParameters(ParamVec);
605 
606  while (i < (int)neigh.Size())
607  {
608  col = i % Largeur;
609  raw = i / Largeur;
610 
611  if ((col > pas && col < Largeur - pas) && (raw > pas && raw < Largeur - pas))
612  {
613  double distanceX = (raw - r);
614  double distanceY = (col - r);
615  dist = std::sqrt(distanceX * distanceX + distanceY * distanceY);
616 
617  if (dist <= r)
618  {
619  /* Transform point to compensate the rotation the orientation */
620  pDst[0] = col;
621  pDst[1] = raw;
622  pSrc = eulerTransform->TransformPoint(pDst);
623 
625  col = static_cast<int>(std::floor(pSrc[0]));
626  raw = static_cast<int>(std::floor(pSrc[1]));
628 
629  if (raw == 0)
630  raw = +1;
631  if (col == 0)
632  col += 1;
633 
634  xx = static_cast<int>(pSrc[1] / rayon);
635  yy = static_cast<int>(pSrc[0] / rayon);
636  Nbin = static_cast<int>(xx + 4 * yy);
637 
638  if (Nbin < 16) // because 64 descriptor length
639  {
640  double distanceXcompensee_2 = (pSrc[0] - r) * (pSrc[0] - r);
641  double distanceYcompensee_2 = (pSrc[1] - r) * (pSrc[1] - r);
642 
643  w = std::exp(-(distanceXcompensee_2 + distanceYcompensee_2) / (2 * 3.3 * 3.3 * S * S));
644 
645  dx = 0.5 * (neigh[(col + pas) + raw * Largeur] - neigh[(col - pas) + raw * Largeur]) * w;
646  dy = 0.5 * (neigh[col + (raw + pas) * Largeur] - neigh[col + (raw - pas) * Largeur]) * w;
647 
648  descriptorVector[4 * Nbin] += dx;
649  descriptorVector[4 * Nbin + 1] += dy;
650  descriptorVector[4 * Nbin + 2] += std::abs(dx);
651  descriptorVector[4 * Nbin + 3] += std::abs(dy);
652  }
653  }
654  }
655  ++i;
656  }
657 
658  double accu = 0;
659  for (int ii = 0; ii < 64; ++ii)
660  accu += descriptorVector[ii] * descriptorVector[ii];
661 
662  for (int jj = 0; jj < 64; ++jj)
663  descriptorVector[jj] /= std::sqrt(accu);
664 
665  return descriptorVector;
666 }
667 
668 /*----------------------------------------------------------------
669  PrintSelf
670  -----------------------------------------------------------------*/
671 template <class TInputImage, class TOutputPointSet>
672 void ImageToSURFKeyPointSetFilter<TInputImage, TOutputPointSet>::PrintSelf(std::ostream& os, itk::Indent indent) const
673 {
674  Superclass::PrintSelf(os, indent);
675  os << indent << "Number of Key Points " << m_NumberOfPoints << std::endl;
676 }
677 }
678 #endif
otb::CONST_PI
constexpr double CONST_PI
Definition: otbMath.h:49
otb::ImageToSURFKeyPointSetFilter::AssignOrientation
virtual double AssignOrientation(const NeighborhoodType &neigh, double S)
otb::ImageToSURFKeyPointSetFilter::GetMin
virtual int GetMin(int a, int b, int c)
otbImageToSURFKeyPointSetFilter.h
otb
The "otb" namespace contains all Orfeo Toolbox (OTB) classes.
Definition: otbJoinContainer.h:32
otb::ImageToSURFKeyPointSetFilter::RefineLocationKeyPoint
bool RefineLocationKeyPoint(const NeighborhoodIteratorType &currentScale, const NeighborhoodIteratorType &previousScale, const NeighborhoodIteratorType &nextScale, VectorPointType &solution)
otb::ImageToSURFKeyPointSetFilter::ImageToSURFKeyPointSetFilter
ImageToSURFKeyPointSetFilter()
otb::MetaDataKey::VectorType
std::vector< double > VectorType
Definition: otbMetaDataKey.h:119
otbGenericMsgDebugMacro
#define otbGenericMsgDebugMacro(x)
Definition: otbMacro.h:63
otb::ImageToSURFKeyPointSetFilter::ComputeDescriptor
virtual VectorType ComputeDescriptor(const NeighborhoodType &neigh, double O, double S)
otb::ImageToSURFKeyPointSetFilter::IsLocalExtremum
virtual bool IsLocalExtremum(const NeighborhoodType &neigh)
otb::ImageToSURFKeyPointSetFilter::IsLocalExtremumAround
virtual bool IsLocalExtremumAround(const NeighborhoodType &neigh, double CenterValue)
otb::ImageToSURFKeyPointSetFilter::GenerateData
void GenerateData() override
otb::ImageToSURFKeyPointSetFilter::PrintSelf
void PrintSelf(std::ostream &os, itk::Indent indent) const override
otb::ImageToSURFKeyPointSetFilter::VectorType
std::vector< double > VectorType
Definition: otbImageToSURFKeyPointSetFilter.h:124
otb::ImageToSURFKeyPointSetFilter::~ImageToSURFKeyPointSetFilter
~ImageToSURFKeyPointSetFilter() override
Generated at Wed Feb 28 2024 15:22:37 for OTB with   doxygen 1.8.17