Orfeo Toolbox  4.0
otbImageToPathListAlignFilter.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 __otbImageToPathListAlignFilter_txx
19 #define __otbImageToPathListAlignFilter_txx
20 
24 #include "itkPathIterator.h"
25 #include "itkImage.h"
26 #include "itkNumericTraits.h"
29 #include "otbMath.h"
30 
31 namespace otb
32 {
33 
35 {
36  short start; /* starting position (distance from border) */
37  short end; /* ending position (hence, length is end-start+1) */
38  double nfa; /* number of false alarms */
39  char ok;
40 };
41 
43 template <class TInputImage, class TOutputPath>
46 {
47  this->SetNumberOfRequiredInputs(1);
48  m_Size.Fill(0);
49  m_isMeaningfulSegment = false;
50  m_NbGradDirection = 16;
51  m_NbLineDirection = 96;
52  m_MinGradNorm = 2.0;
53  m_Eps = 0.0;
54 
55  for (unsigned int i = 0; i < InputImageDimension; ++i)
56  {
57  // Set an image spacing for the user
58  m_Spacing[i] = 1.0;
59  m_Origin[i] = 0;
60  }
61 
63  m_BackgroundValue = itk::NumericTraits<ValueType>::Zero;
64 }
65 
67 template <class TInputImage, class TOutputPath>
70 {
71 }
72 
73 //----------------------------------------------------------------------------
74 template <class TInputImage, class TOutputPath>
75 void
77 ::SetSpacing(const double* spacing)
78 {
79  unsigned int i;
80  for (i = 0; i < InputImageDimension; ++i)
81  {
82  if (spacing[i] != m_Spacing[i])
83  {
84  break;
85  }
86  }
87  if (i < InputImageDimension)
88  {
89  for (i = 0; i < InputImageDimension; ++i)
90  {
91  m_Spacing[i] = spacing[i];
92  }
93  this->Modified();
94  }
95 }
96 
97 template <class TInputImage, class TOutputPath>
98 void
100 ::SetSpacing(const float* spacing)
101 {
102  unsigned int i;
103  for (i = 0; i < InputImageDimension; ++i)
104  {
105  if ((double) spacing[i] != m_Spacing[i])
106  {
107  break;
108  }
109  }
110  if (i < InputImageDimension)
111  {
112  for (i = 0; i < InputImageDimension; ++i)
113  {
114  m_Spacing[i] = spacing[i];
115  }
116  this->Modified();
117  }
118 }
119 
120 template <class TInputImage, class TOutputPath>
121 const double *
124 {
125  return m_Spacing;
126 }
127 
128 //----------------------------------------------------------------------------
129 template <class TInputImage, class TOutputPath>
130 void
132 ::SetOrigin(const double* origin)
133 {
134  unsigned int i;
135  for (i = 0; i < InputImageDimension; ++i)
136  {
137  if (origin[i] != m_Origin[i])
138  {
139  break;
140  }
141  }
142  if (i < InputImageDimension)
143  {
144  for (i = 0; i < InputImageDimension; ++i)
145  {
146  m_Origin[i] = origin[i];
147  }
148  }
149 }
150 
151 template <class TInputImage, class TOutputPath>
152 void
154 ::SetOrigin(const float* origin)
155 {
156  unsigned int i;
157  for (i = 0; i < InputImageDimension; ++i)
158  {
159  if ((double) origin[i] != m_Origin[i])
160  {
161  break;
162  }
163  }
164  if (i < InputImageDimension)
165  {
166  for (i = 0; i < InputImageDimension; ++i)
167  {
168  m_Origin[i] = origin[i];
169  }
170  }
171 }
172 
173 template <class TInputImage, class TOutputPath>
174 const double *
176 ::GetOrigin() const
177 {
178  return m_Origin;
179 }
180 
181 //----------------------------------------------------------------------------
182 /* Algorithm */
183 template <class TInputImage, class TOutputPath>
184 std::vector<double>
186 ::tab(int n, double p, double m)
187 {
188  std::vector<double> out;
189  int adr1, adr2, x, y;
190 // double lambda;
191  double q;
192 
193  q = 1.0 - p;
194  out.resize((n + 1) * (n + 1));
195  adr1 = 0;
196 
197  /*** compute proba (=x among y) ***/
198  out[0] = 1.0;
199  for (y = 1, adr2 = 0; y <= n; ++y)
200  {
201  adr1 = adr2;
202  adr2 += n + 1;
203  out[adr2] = q * out[adr1];
204  for (x = 1; x <= y; ++x)
205  out[adr2 + x] = p * out[adr1 + x - 1] + q * out[adr1 + x];
206  }
207 
208  /*** sum to obtain proba (>=k among y) ***/
209  for (y = 1, adr1 = n + 1; y <= n; ++y, adr1 += n + 1)
210  for (x = y - 1; x >= 0; x--)
211  out[adr1 + x] += out[adr1 + x + 1];
212 
213  /*** multiply by m (number of segments) to obtain expectation***/
214  for (adr1 = (n + 1) * (n + 1); --adr1 >= 0; )
215  out[adr1] *= m;
216 
217  return out;
218 }
219 
220 template <class TInputImage, class TOutputPath>
221 void
224 {
225  double threshold;
226  int n, p, x, y;
227 
228  typename InputImageType::SizeType Taille;
229  typename RealImageType::IndexType IndexOut;
230 
231  Taille = InputImage->GetLargestPossibleRegion().GetSize();
232 
233  typename RealImageType::RegionType region;
234  region.SetSize(InputImage->GetLargestPossibleRegion().GetSize());
235  IndexOut[0] = 0;
236  IndexOut[1] = 0;
237 // region.SetIndex(InputImage->GetLargestPossibleRegion().GetIndex());
238  region.SetIndex(IndexOut);
239  m_AngleImage->SetRegions(region);
240  m_AngleImage->SetOrigin(InputImage->GetOrigin());
241  m_AngleImage->SetSpacing(InputImage->GetSpacing());
242  m_AngleImage->Allocate();
243 
244  n = Taille[0];
245  p = Taille[1];
246 
247  threshold = m_MinGradNorm;
248  threshold *= threshold;
249 
250  typename InputImageType::IndexType idx;
251 
252  for (x = 0; x < p; ++x)
253  {
254  idx[0] = (n - 1);
255  idx[1] = x;
256 // indice = (n-1)*p +x
257  m_AngleImage->SetPixel(idx, static_cast<RealType>(-1000.0));
258  }
259  for (y = 0; y < n; ++y)
260  {
261  idx[0] = y;
262  idx[1] = p - 1;
263 // indice = p*y+p-1
264  m_AngleImage->SetPixel(idx, static_cast<RealType>(-1000.0));
265  }
266 
267  typename InputImageType::IndexType adr;
268  RealType PixelA, PixelB, PixelC, PixelD;
269  RealType com1, com2, gx, gy, norm;
270 
271  for (x = 0; x < p - 1; ++x)
272  for (y = 0; y < n - 1; ++y)
273  {
274 // indice = y*p+x
275  adr[0] = y;
276  adr[1] = x;
277  idx[0] = adr[0] + 1;
278  idx[1] = adr[1] + 1;
279  PixelA = static_cast<RealType>(InputImage->GetPixel(idx));
280  idx[0] = adr[0];
281  idx[1] = adr[1];
282  assert(idx[0] < n);
283  assert(idx[1] < p);
284  assert(idx[0] >= 0);
285  assert(idx[1] >= 0);
286  PixelB = static_cast<RealType>(InputImage->GetPixel(idx));
287  idx[0] = adr[0] + 1;
288  idx[1] = adr[1];
289  assert(idx[0] < n);
290  assert(idx[1] < p);
291  assert(idx[0] >= 0);
292  assert(idx[1] >= 0);
293  PixelC = static_cast<RealType>(InputImage->GetPixel(idx));
294  idx[0] = adr[0];
295  idx[1] = adr[1] + 1;
296  assert(idx[0] < n);
297  assert(idx[1] < p);
298  assert(idx[0] >= 0);
299  assert(idx[1] >= 0);
300  PixelD = static_cast<RealType>(InputImage->GetPixel(idx));
301  com1 = PixelA - PixelB;
302  com2 = PixelC - PixelD;
303  gx = 0.5 * (com1 + com2);
304  gy = 0.5 * (com1 - com2);
305  norm = gx * gx + gy * gy;
306 
307  if (norm <= threshold) m_AngleImage->SetPixel(adr, static_cast<RealType>(-1000.0));
308  else m_AngleImage->SetPixel(adr, static_cast<RealType>(vcl_atan2(gx, -gy)));
309  }
310 }
311 
312 template <class TInputImage, class TOutputPath>
313 void
316 {
317 // SizeType size;
318 // double origin[InputImageDimension];
319  typename InputImageType::SizeType Taille;
320  RealImageTypeIndexType indexAngle;
321 // Flist result;
322  int iseglist, size_seglist; /* associated counter and dynamic size */
323  int iseg, size_seg;
324  double nfa, max_nfa;
325  std::vector<double> test;
326  std::vector<int> count, startbloc, endbloc;
327  std::vector<double> seglist; /* list of recorded segments */
328  std::vector<one_segment> seg;
329  int mx, my, ox, oy, nx, ny, n;
330  int xx, yy, pos, posmax, nblocs, inbloc, max_nblocs;
331  int cur, i, j, side, l, lphase;
332  // int tmp;
333  int itheta, ntheta;
334  double theta, theta0, dtheta, dx, dy, prec;
335  double error = 0.0;
336  itkDebugMacro(<< "ImageToPathListAlignFilter::GenerateData() called");
337 
338  // Get the input and output pointers
339  const InputImageType * InputImage = this->GetInput();
340  OutputPathListType * OutputPath = this->GetOutput();
341  // Generate the image
342 
343  /* Filter algorithm */
344 
345  Taille = InputImage->GetLargestPossibleRegion().GetSize();
346  nx = Taille[0];
347  ny = Taille[1];
348  max_nfa = vcl_pow(10.0, -(m_Eps));
349 
350 // typename InputImageType::IndexType adr;
351 
352  /*** maximal length for a line */
353  n = (int) vcl_ceil(hypot((double) nx, (double) ny)) + 1;
354 
355  /*** compute angle map of u ***/
356  RealImageTypePointer lAngleImagePointer = RealImageType::New();
357  m_AngleImage = static_cast<RealImageType*>(lAngleImagePointer.GetPointer());
358  this->AngleCalculate(InputImage);
359 
360  /*** compute P(k, l) ***/
361  test = tab(n, 1.0 / (double) (m_NbGradDirection), (double) (nx * ny) * (double) (nx * ny));
362 
363  /*** initialization ***/
364  prec = CONST_PI / (double) (m_NbGradDirection);
365  ntheta = m_NbLineDirection / 2; /* i.e. # directions of NON-ORIENTED lines */
366  dtheta = CONST_PI / (double) ntheta;
367 
368  /******************** memory allocation ********************/
369 
370  max_nblocs = n / 2 + 1; /* maximal number of blocs */
371  count.resize(max_nblocs);
372  startbloc.resize(max_nblocs);
373  endbloc.resize(max_nblocs);
374 
375  size_seg = 10000; /* initial allocation (may reallocate later) */
376  seg.resize(size_seg);
377 
378  size_seglist = 10000; /* initial allocation (may reallocate later) */
379  seglist.resize(5 * size_seglist);
380 
381  /* counter for recorded segments (seglist) */
382  iseglist = 0;
383 
384  /******************** first loop : the four sides ********************/
385 
386  for (side = 0; side < 4; side++)
387  {
388  printf("side %d/4 ", side + 1);
389 
390  theta0 = CONST_PI_2 * (double) side;
391  mx = ((side == 0 || side == 2) ? 1 : 0);
392  my = ((side == 1 || side == 3) ? 1 : 0);
393  ox = ((side == 1) ? nx - 1 : 0);
394  oy = ((side == 2) ? ny - 1 : 0);
395 
396  posmax = nx * mx + ny * my;
397 
398  /*** second loop : angles ***/
399  for (itheta = 0; itheta < ntheta; itheta++)
400  {
401  printf(".");
402  fflush(stdout);
403  theta = theta0 + (double) (itheta) * dtheta;
404  dx = (double) vcl_cos((double) theta);
405  dy = (double) vcl_sin((double) theta);
406 
407  /*** third loop : start positions ***/
408  for (pos = 0; pos < posmax; ++pos)
409  {
410 
411  /* clear segment array */
412  iseg = 0;
413 
414  /*** fourth loop : phase for two-spaced pixels ***/
415  for (lphase = 0; lphase < 2; lphase++)
416  {
417 
418  /*** detect aligned points by blocs ***/
419  inbloc = nblocs = cur = l = count[0] = 0;
420  xx = ox + pos * mx + (int) (dx * (double) (l * 2 + lphase));
421  yy = oy + pos * my + (int) (dy * (double) (l * 2 + lphase));
422 
423  for (; xx >= 0 && xx < nx && yy >= 0 && yy < ny; )
424  {
425  indexAngle[0] = xx;
426  indexAngle[1] = yy;
427  // indice = yy*nx+xx
428  assert(indexAngle[0] < nx);
429  assert(indexAngle[1] < ny);
430  assert(indexAngle[0] >= 0);
431  assert(indexAngle[1] >= 0);
432 
433  error = static_cast<double>(m_AngleImage->GetPixel(indexAngle));
434  if (error > -100.0)
435  {
436  error -= theta;
437  while (error <= -CONST_PI)
438  error += CONST_2PI;
439  while (error > CONST_PI)
440  error -= CONST_2PI;
441  if (error < 0.0) error = -error;
442  if (error < prec)
443  {
444  ++cur;
445  if (!inbloc)
446  {
447  startbloc[nblocs] = l;
448  inbloc = 1;
449  }
450  }
451  else
452  {
453  if (inbloc)
454  {
455  endbloc[nblocs] = l - 1;
456  ++nblocs;
457  count[nblocs] = cur;
458  }
459  inbloc = 0;
460  }
461  }
462  /* compute next point */
463  ++l;
464  xx = ox + pos * mx + (int) (dx * (double) (l * 2 + lphase));
465  yy = oy + pos * my + (int) (dy * (double) (l * 2 + lphase));
466  }
467 
468  /*** detect meaningful segments ***/
469  for (i = 0; i < nblocs; ++i)
470  for (j = i; j < nblocs; ++j)
471  if ((nfa = test[count[j + 1] - count[i]
472  + (n + 1) * (1 + endbloc[j] - startbloc[i])]) < max_nfa)
473  {
474  seg[iseg].start = startbloc[i] * 2 + lphase;
475  seg[iseg].end = endbloc[j] * 2 + lphase;
476  seg[iseg].nfa = nfa;
477  seg[iseg].ok = 1;
478  iseg++;
479  /* reallocate if necessary */
480  if (iseg == size_seg)
481  {
482  size_seg = (size_seg * 3) / 2;
483  seg.resize(size_seg);
484 // if (!seg)
485 // mwerror(FATAL, 1,"Not enough memory.");
486  }
487  }
488  }
489  /*** end of phase loop ***/
490 
491  /*** remove non-maximal segments ***/
492  if (!m_isMeaningfulSegment)
493  for (i = 0; i < iseg; ++i)
494  for (j = 0; j < iseg; ++j)
495  if (i != j)
496 
497  /* seg[i] is included in seg[j] ? */
498  if (seg[i].start >= seg[j].start && seg[i].end <= seg[j].end)
499  {
500 
501  /* remove the less meaningful of seg[i] and seg[j] */
502  if (seg[i].nfa < seg[j].nfa) seg[j].ok = 0;
503  else seg[i].ok = 0;
504 
505  }
506 
507  /*** store detected segments ***/
508  for (i = 0; i < iseg; ++i)
509  if (seg[i].ok)
510  {
511  seglist[iseglist * 5] = (double) (ox + pos * mx) + dx * (double) (seg[i].start);
512  seglist[iseglist * 5 + 1] = (double) (oy + pos * my) + dy * (double) (seg[i].start);
513  seglist[iseglist * 5 + 2] = (double) (ox + pos * mx) + dx * (double) (seg[i].end);
514  seglist[iseglist * 5 + 3] = (double) (oy + pos * my) + dy * (double) (seg[i].end);
515  seglist[iseglist * 5 + 4] = -(double) log10(seg[i].nfa);
516  iseglist++;
517  /* reallocate seglist if necessary */
518  if (iseglist == size_seglist)
519  {
520  size_seglist = (size_seglist * 3) / 2;
521  seglist.resize(size_seglist);
522 // if (!seglist)
523 // mwerror(FATAL, 1,"Not enough memory.");
524  }
525  }
526  }
527  }
528  /*** end of second loop ***/
529 
530  printf(" nb de segments: %d\n", iseglist);
531  }
532  /******************** end of first loop ********************/
533 
534  seg.clear();
535  endbloc.clear();
536  startbloc.clear();
537  count.clear();
538  test.clear();
539 
540  /* build segments list */
541  seglist.resize(5 * iseglist);
542 
543  /* build segments list */
544  OutputPath->Clear();
545 // OutputPath->Resize(iseglist);
546 
547  typedef typename OutputPathType::ContinuousIndexType ContinuousIndexType;
548  typename InputImageType::PointType point;
549 
550  ContinuousIndexType cindex;
551  for (i = 0; i < iseglist; ++i)
552  {
553 
554  OutputPathPointerType path = OutputPathType::New();
555 
556  path->Initialize();
557  point[0] = seglist[i * 5];
558  point[1] = seglist[i * 5 + 1];
559  InputImage->TransformPhysicalPointToContinuousIndex(point, cindex);
560  path->AddVertex(cindex);
561  cindex[0] = seglist[i * 5 + 2];
562  cindex[1] = seglist[i * 5 + 3];
563  path->AddVertex(cindex);
564 
565  OutputPath->PushBack(path);
566  }
567  itkDebugMacro(<< "ImageToPathListAlignFilter::GenerateData() finished");
568 
569 } // end update function
570 
571 template <class TInputImage, class TOutputPath>
572 void
574 ::PrintSelf(std::ostream& os, itk::Indent indent) const
575 {
576  Superclass::PrintSelf(os, indent);
577 // os << indent << "Size : " << m_Size << std::endl;
578 // os << indent << "Path Value : " << m_PathValue << std::endl;
579 // os << indent << "Background Value : " << m_BackgroundValue << std::endl;
580 }
581 
582 } // end namespace otb
583 
584 #endif

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