KSquare Utilities
ContourFollower.cpp
Go to the documentation of this file.
1 /* ContourFollower.cpp -- Used to find the contour of image in a Raster object.
2  * Copyright (C) 1994-2011 Kurt Kramer
3  * For conditions of distribution and use, see copyright notice in KKB.h
4  */
5 #include "FirstIncludes.h"
6 #include <stdlib.h>
7 #include <memory>
8 #include <math.h>
9 #include <complex>
10 #include <fstream>
11 #include <iostream>
12 #include <string>
13 #include <vector>
14 #include "MemoryDebug.h"
15 using namespace std;
16 
17 #if defined(FFTW_AVAILABLE)
18 # include <fftw3.h>
19 //#else
20 //# include "kku_fftw.h"
21 #endif
22 #include "kku_fftw.h"
23 
24 
25 
26 #include "ContourFollower.h"
27 #include "Raster.h"
28 #include "KKStr.h"
29 using namespace KKB;
30 
31 
32 
33 const
34 MovDir movements[8] = {{-1, 0}, // Up
35  {-1, 1}, // Up-Right
36  { 0, 1}, // Right
37  { 1, 1}, // Down-Right
38  { 1, 0}, // Down
39  { 1, -1}, // Down-Left
40  { 0, -1}, // Left,
41  {-1, -1} // Up-Left
42  };
43 
44 
45 
46 
48  RunLog& _log
49  ):
50  backgroundPixelTH (_raster.BackgroundPixelTH ()),
51  curCol (-1),
52  curRow (-1),
53  fromDir (),
54  height (_raster.Height ()),
55  lastDir (0),
56  log (_log),
57  raster (_raster),
58  rows (_raster.Rows ()),
59  width (_raster.Width ())
60 {
61 }
62 
63 
64 
65 uchar ContourFollower::PixelValue (kkint32 row,
66  kkint32 col
67  )
68 {
69  if ((row < 0) || (row >= height)) return 0;
70  if ((col < 0) || (col >= width)) return 0;
71  return rows[row][col];
72 }
73 
74 
75 kkint32 ContourFollower::PixelCountIn9PixelNeighborhood (kkint32 row,
76  kkint32 col
77  )
78 {
79  kkint32 count = 0;
80  if (PixelValue (row - 1, col - 1) > backgroundPixelTH) ++count;
81  if (PixelValue (row - 1, col ) > backgroundPixelTH) ++count;
82  if (PixelValue (row - 1, col + 1) > backgroundPixelTH) ++count;
83  if (PixelValue (row , col - 1) > backgroundPixelTH) ++count;
84  if (PixelValue (row , col ) > backgroundPixelTH) ++count;
85  if (PixelValue (row , col + 1) > backgroundPixelTH) ++count;
86  if (PixelValue (row + 1, col - 1) > backgroundPixelTH) ++count;
87  if (PixelValue (row + 1, col ) > backgroundPixelTH) ++count;
88  if (PixelValue (row + 1, col + 1) > backgroundPixelTH) ++count;
89  return count;
90 }
91 
92 
93 
94 void ContourFollower::GetFirstPixel (kkint32& startRow,
95  kkint32& startCol
96  )
97 {
98  lastDir = 1;
99  fromDir = 0;
100 
101  startRow = height / 2;
102  startCol = 0;
103 
104  while (startCol < width)
105  {
106  if ((rows[startRow][startCol] > backgroundPixelTH) && (PixelCountIn9PixelNeighborhood (startRow, startCol) > 1))
107  break;
108  startCol++;
109  }
110 
111  if (startCol >= width)
112  {
113  // We did not find a Pixel in the Middle Row(height / 2) so now we will
114  // scan the image row, by row, Left to Right until we find an occupied
115  // pixel.
116 
117  bool found = false;
118  kkint32 row;
119  kkint32 col;
120 
121  for (row = 0; ((row < height) && (!found)); ++row)
122  {
123  for (col = 0; ((col < width) && (!found)); ++col)
124  {
125  if ((rows[row][col] > backgroundPixelTH) && (PixelCountIn9PixelNeighborhood (row, col) > 1))
126  {
127  found = true;
128  startRow = row;
129  startCol = col;
130  }
131  }
132  }
133 
134  if (!found)
135  {
136  startRow = startCol = -1;
137  return;
138  }
139  }
140 
141  curRow = startRow;
142  curCol = startCol;
143 } /* GetFirstPixel */
144 
145 
146 
147 void ContourFollower::GetNextPixel (kkint32& nextRow,
148  kkint32& nextCol
149  )
150 {
151  fromDir = lastDir + 4;
152 
153  if (fromDir > 7)
154  fromDir = fromDir - 8;
155 
156  bool nextPixelFound = false;
157 
158  kkint32 nextDir = fromDir + 2;
159 
160  while (!nextPixelFound)
161  {
162  if (nextDir > 7)
163  nextDir = nextDir - 8;
164 
165  nextRow = curRow + movements[nextDir].row;
166  nextCol = curCol + movements[nextDir].col;
167 
168  if ((nextRow < 0) ||
169  (nextRow >= height) ||
170  (nextCol < 0) ||
171  (nextCol >= width)
172  )
173  {
174  nextDir++;
175  }
176  else if (rows[nextRow][nextCol] > backgroundPixelTH)
177  {
178  nextPixelFound = true;
179  }
180  else
181  {
182  nextDir++;
183  }
184  }
185 
186  curRow = nextRow;
187  curCol = nextCol;
188 
189  lastDir = nextDir;
190 } /* GetNextPixel */
191 
192 
193 
194 #if defined(FFTW_AVAILABLE)
195  float CalcMagnitude (fftwf_complex* dest,
196  kkint32 index
197  )
198  {
199  float mag = 0.0f;
200  mag = (float)sqrt (dest[index][0] * dest[index][0] + dest[index][1] * dest[index][1]);
201  return mag;
202  } /* CalcMagnitude */
203 #else
204  float CalcMagnitude (KK_DFT1D_Float::DftComplexType* dest,
205  kkint32 index
206  )
207  {
208  double rp = (float)dest[index].real ();
209  double ip = (float)dest[index].imag ();
210  return (float)sqrt (rp * rp + ip * ip);
211  } /* CalcMagnitude */
212 
213 #endif
214 
215 
216 
217 
218 
219 
220 kkint32 ContourFollower::FollowContour (float* countourFreq,
221  float fourierDescriptors[15],
222  kkint32 totalPixels,
223  bool& successful
224  )
225 {
226  // startRow and startCol is assumed to come from the left (6)
227 
228  // I am having a problem;
229  // We sometimes locate a separate piece of the image and end of with a contour of only that
230  // one piece. This results in two problems.
231  // 1) The Fourier Descriptors do not even begin to represent the image
232  // 2) If the number of border pixels are small enough; then we could end
233  // up causing a memory access fault when we access a negative index.
234 
235 
236  kkint32 numOfBuckets = 5;
237 
238  kkint32 startRow;
239  kkint32 startCol;
240 
241  kkint32 scndRow;
242  kkint32 scndCol;
243 
244  kkint32 lastRow;
245  kkint32 lastCol;
246 
247  kkint32 nextRow;
248  kkint32 nextCol;
249 
250  kkint32 absoluteMaximumEdgePixels = totalPixels * 3;
251 
252  kkint32 maxNumOfBorderPoints = 3 * (height + width);
253  kkint32 numOfBorderPixels = 0;
254 
255  kkint32 x;
256 
257  successful = true;
258 
259  kkint32 totalRow = 0;
260  kkint32 totalCol = 0;
261 
262  #if defined(FFTW_AVAILABLE)
263  fftwf_complex* src = (fftwf_complex*)fftwf_malloc (sizeof (fftwf_complex) * maxNumOfBorderPoints);
264  #else
265  KK_DFT1D_Float::DftComplexType* src = new KK_DFT1D_Float::DftComplexType[maxNumOfBorderPoints];
266  #endif
267 
268 
269  GetFirstPixel (startRow, startCol);
270  if ((startRow < 0) || (startCol < 0) || (PixelCountIn9PixelNeighborhood (startRow, startCol) < 2))
271  {
272  cout << "Very Bad Starting Point" << std::endl;
273  successful = false;
274  return 0;
275  }
276 
277  GetNextPixel (scndRow, scndCol);
278 
279  nextRow = scndRow;
280  nextCol = scndCol;
281 
282  lastRow = startRow;
283  lastCol = startCol;
284 
285 
286  while (true)
287  {
288  if (numOfBorderPixels > absoluteMaximumEdgePixels)
289  {
290  // Something must be wrong; somehow we missed the starting point.
291 
292  log.Level (-1) << endl << endl << endl
293  << "ContourFollower::FollowContour ***ERROR***" << endl
294  << endl
295  << " We exceeded the absolute maximum number of possible edge pixels." << endl
296  << endl
297  << " FileName [" << raster.FileName () << "]" << endl
298  << " numOfBorderPixels [" << numOfBorderPixels << "]" << endl
299  << endl;
300 
301  ofstream r ("c:\\temp\\ContourFollowerFollowContour.log", std::ios_base::ate);
302  r << "FileName [" << raster.FileName () << "] ";
303  r << "totalPixels [" << totalPixels << "] ";
304  r << "numOfBorderPixels [" << numOfBorderPixels << "]";
305  r << std::endl;
306  r.close ();
307  break;
308  }
309 
310  lastRow = nextRow;
311  lastCol = nextCol;
312 
313  GetNextPixel (nextRow, nextCol);
314 
315  if ((nextRow == scndRow) && (nextCol == scndCol) &&
316  (lastRow == startRow) && (lastCol == startCol))
317  {
318  break;
319  }
320 
321  if (numOfBorderPixels >= maxNumOfBorderPoints)
322  {
323  kkint32 newMaxNumOfAngles = maxNumOfBorderPoints * 2;
324 
325  #if defined(FFTW_AVAILABLE)
326  fftwf_complex* newSrc = (fftwf_complex*)fftwf_malloc (sizeof (fftwf_complex) * newMaxNumOfAngles);
327  #else
328  KK_DFT1D_Float::DftComplexType* newSrc = new KK_DFT1D_Float::DftComplexType[newMaxNumOfAngles];
329  #endif
330 
331  if (newSrc == NULL)
332  {
333  log.Level (-1) << endl << endl
334  << "FollowContour ***ERROR***" << endl
335  << endl
336  << "Could not allocate memory needed for contour points." << endl
337  << endl;
338  successful = false;
339  return 0;
340  }
341 
342  kkint32 x;
343  for (x = 0; x < maxNumOfBorderPoints; x++)
344  {
345  #if defined(FFTW_AVAILABLE)
346  newSrc[x][0] = src[x][0];
347  newSrc[x][1] = src[x][1];
348  #else
349  newSrc[x].real (src[x].real ());
350  newSrc[x].imag (src[x].imag ());
351  #endif
352  }
353 
354 
355  #if defined(FFTW_AVAILABLE)
356  fftwf_free (src);
357  src = newSrc;
358  newSrc = NULL;
359  #else
360  delete src;
361  src = newSrc;
362  newSrc = NULL;
363  #endif
364 
365  maxNumOfBorderPoints = newMaxNumOfAngles;
366  }
367 
368  #if defined(FFTW_AVAILABLE)
369  src[numOfBorderPixels][0] = (float)nextRow;
370  src[numOfBorderPixels][1] = (float)nextCol;
371  #else
372  src[numOfBorderPixels].real ((float)nextRow);
373  src[numOfBorderPixels].imag ((float)nextCol);
374  #endif
375 
376  totalRow += nextRow;
377  totalCol += nextCol;
378 
379  numOfBorderPixels++;
380  }
381 
382  float centerRow = (float)totalRow / (float)numOfBorderPixels;
383  float centerCol = (float)totalCol / (float)numOfBorderPixels;
384 
385  float totalRe = 0.0;
386  float totalIm = 0.0;
387 
388  for (x = 0; x < numOfBorderPixels; x++)
389  {
390  #if defined(FFTW_AVAILABLE)
391  src[x][0] = (src[x][0] - centerRow);
392  src[x][1] = (src[x][1] - centerCol);
393  totalRe += src[x][0];
394  totalIm += src[x][1];
395  #else
396  src[x].real ((src[x].real () - centerRow));
397  src[x].imag ((src[x].imag () - centerCol));
398  totalRe += src[x].real ();
399  totalIm += src[x].imag ();
400  #endif
401  }
402 
403  #if defined(FFTW_AVAILABLE)
404  fftwf_complex* dest = (fftwf_complex*)fftwf_malloc (sizeof (fftwf_complex) * maxNumOfBorderPoints);
405  fftwf_plan plan = fftwCreateOneDPlan (numOfBorderPixels, src, dest, FFTW_FORWARD, FFTW_ESTIMATE);
406  fftwf_execute (plan);
407  fftwDestroyPlan (plan);
408  plan = NULL;
409  #else
410  KK_DFT1D_Float plan (numOfBorderPixels, true);
411  KK_DFT1D_Float::DftComplexType* dest = new KK_DFT1D_Float::DftComplexType[numOfBorderPixels];
412  plan.Transform (src, dest);
413  #endif
414 
415  kkint32 numOfedgePixels = numOfBorderPixels;
416 
417  kkint32* count = new kkint32 [numOfBuckets];
418 
419  for (x = 0; x < numOfBuckets; x++)
420  {
421  countourFreq[x] = 0.0f;
422  count[x] = 0;
423  }
424 
425  // Original Region Areas, as reflected in ICPR paper
426  float middle = numOfBorderPixels / (float)2.0;
427  float r1 = middle / (float)2.0;
428  float r2 = middle * ( (float)3.0 / (float)4.0);
429  float r3 = middle * ( (float)7.0 / (float)8.0);
430  float r4 = middle * ((float)15.0 / (float)16.0);
431 
432  float deltaX = 0.0f;
433  float mag = 0.0f;
434 
435  if (numOfBorderPixels < 8)
436  {
437  successful = false;
438  }
439  else
440  {
441  float normalizationFactor = CalcMagnitude (dest, 1);
442 
443  fourierDescriptors[ 0] = CalcMagnitude (dest, numOfBorderPixels - 1) / normalizationFactor;
444 
445  fourierDescriptors[ 1] = CalcMagnitude (dest, 2) / normalizationFactor;
446  fourierDescriptors[ 2] = CalcMagnitude (dest, numOfBorderPixels - 2) / normalizationFactor;
447 
448  fourierDescriptors[ 3] = CalcMagnitude (dest, 3) / normalizationFactor;
449  fourierDescriptors[ 4] = CalcMagnitude (dest, numOfBorderPixels - 3) / normalizationFactor;
450 
451  fourierDescriptors[ 5] = CalcMagnitude (dest, 4) / normalizationFactor;
452  fourierDescriptors[ 6] = CalcMagnitude (dest, numOfBorderPixels - 4) / normalizationFactor;
453 
454  fourierDescriptors[ 7] = CalcMagnitude (dest, 5) / normalizationFactor;
455  fourierDescriptors[ 8] = CalcMagnitude (dest, numOfBorderPixels - 5) / normalizationFactor;
456 
457  fourierDescriptors[ 9] = CalcMagnitude (dest, 6) / normalizationFactor;
458  fourierDescriptors[10] = CalcMagnitude (dest, numOfBorderPixels - 6) / normalizationFactor;
459 
460  fourierDescriptors[11] = CalcMagnitude (dest, 7) / normalizationFactor;
461  fourierDescriptors[12] = CalcMagnitude (dest, numOfBorderPixels - 7) / normalizationFactor;
462 
463  fourierDescriptors[13] = CalcMagnitude (dest, 8) / normalizationFactor;
464  fourierDescriptors[14] = CalcMagnitude (dest, numOfBorderPixels - 8) / normalizationFactor;
465  }
466 
467  for (x = 1; x < (numOfBorderPixels - 1); x++)
468  {
469  //mag = log (sqrt (dest[x].re * dest[x].re + dest[x].im * dest[x].im));
470  //mag = log (sqrt (dest[x].re * dest[x].re + dest[x].im * dest[x].im) + 1.0);
471 
472  mag = CalcMagnitude (dest, x);
473 
474  deltaX = (float)fabs ((float)x - middle);
475 
476  if (deltaX < r1)
477  {
478  countourFreq[0] = countourFreq[0] + mag;
479  count[0]++;
480  }
481 
482  else if (deltaX < r2)
483  {
484  countourFreq[1] = countourFreq[1] + mag;
485  count[1]++;
486  }
487 
488  else if (deltaX < r3)
489  {
490  countourFreq[2] = countourFreq[2] + mag;
491  count[2]++;
492  }
493 
494  else if (deltaX < r4)
495  {
496  countourFreq[3] = countourFreq[3] + mag;
497  count[3]++;
498  }
499 
500  else
501  {
502  countourFreq[4] = countourFreq[4] + mag;
503  count[4]++;
504  }
505  }
506 
507  for (x = 0; x < numOfBuckets; ++x)
508  {
509  if (count[x] <= 0)
510  {
511  countourFreq[x] = 0.0;
512  }
513  else
514  {
515  countourFreq[x] = countourFreq[x] / (float)count[x];
516  }
517  }
518 
519  #if defined(FFTW_AVAILABLE)
520  fftwf_free (src); src = NULL;
521  fftwf_free (dest); dest = NULL;
522  #else
523  delete src; src = NULL;
524  delete dest; dest = NULL;
525  #endif
526  delete[] count; count = NULL;
527 
528  return numOfedgePixels;
529 } /* FollowContour */
530 
531 
532 
533 kkint32 ContourFollower::FollowContour2 (float* countourFreq,
534  bool& successful
535  )
536 {
537  // A different Version of FollowContour
538 
539  // startRow and startCol is assumed to come from the left (6)
540 
541  kkint32 numOfBuckets = 16;
542 
543  kkint32 startRow;
544  kkint32 startCol;
545 
546  kkint32 scndRow;
547  kkint32 scndCol;
548 
549  kkint32 lastRow;
550  kkint32 lastCol;
551 
552  kkint32 nextRow;
553  kkint32 nextCol;
554 
555  kkint32 maxNumOfBorderPoints = 3 * (height + width);
556  kkint32 numOfBorderPixels = 0;
557 
558  kkint32 x;
559 
560  successful = true;
561 
562  kkint32 totalRow = 0;
563  kkint32 totalCol = 0;
564 
565 
566  #if defined(FFTW_AVAILABLE)
567  fftwf_complex* src = (fftwf_complex*)fftwf_malloc (sizeof (fftwf_complex) * maxNumOfBorderPoints);
568  #else
569  KK_DFT1D_Float::DftComplexType* src = new KK_DFT1D_Float::DftComplexType[maxNumOfBorderPoints];
570  #endif
571 
572  GetFirstPixel (startRow, startCol);
573  if ((startRow < 0) || (startCol < 0) || (PixelCountIn9PixelNeighborhood (startRow, startCol) < 2))
574  {
575  cout << "Very Bad Starting Point" << std::endl;
576  successful = false;
577  return 0;
578  }
579 
580  GetNextPixel (scndRow, scndCol);
581 
582  nextRow = scndRow;
583  nextCol = scndCol;
584 
585  lastRow = startRow;
586  lastCol = startCol;
587 
588  while (true)
589  {
590  lastRow = nextRow;
591  lastCol = nextCol;
592 
593  GetNextPixel (nextRow, nextCol);
594 
595  if ((nextRow == scndRow) && (nextCol == scndCol) &&
596  (lastRow == startRow) && (lastCol == startCol))
597  {
598  break;
599  }
600 
601  if (numOfBorderPixels >= maxNumOfBorderPoints)
602  {
603  kkint32 newMaxNumOfAngles = maxNumOfBorderPoints * 2;
604 
605  #if defined(FFTW_AVAILABLE)
606  fftwf_complex* newSrc = (fftwf_complex*)fftwf_malloc (sizeof (fftwf_complex) * newMaxNumOfAngles);
607  #else
608  KK_DFT1D_Float::DftComplexType* newSrc = new KK_DFT1D_Float::DftComplexType[newMaxNumOfAngles];
609  #endif
610 
611  kkint32 x;
612  for (x = 0; x < maxNumOfBorderPoints; x++)
613  {
614  #if defined(FFTW_AVAILABLE)
615  newSrc[x][0] = src[x][0];
616  newSrc[x][1] = src[x][1];
617  #else
618  newSrc[x].real (src[x].real ());
619  newSrc[x].imag (src[x].imag ());
620  #endif
621  }
622 
623  #if defined(FFTW_AVAILABLE)
624  fftwf_free (src);
625  src = newSrc;
626  newSrc = NULL;
627  #else
628  delete src;
629  src = newSrc;
630  newSrc = NULL;
631  #endif
632 
633  maxNumOfBorderPoints = newMaxNumOfAngles;
634  }
635 
636  #if defined(FFTW_AVAILABLE)
637  src[numOfBorderPixels][0] = (float)nextRow;
638  src[numOfBorderPixels][1] = (float)nextCol;
639  #else
640  src[numOfBorderPixels].real ((float)nextRow);
641  src[numOfBorderPixels].imag ((float)nextCol);
642  #endif
643 
644  totalRow += nextRow;
645  totalCol += nextCol;
646 
647  numOfBorderPixels++;
648  }
649 
650  float centerRow = (float)totalRow / (float)numOfBorderPixels;
651  float centerCol = (float)totalCol / (float)numOfBorderPixels;
652 
653  float totalRe = 0.0;
654  float totalIm = 0.0;
655 
656  for (x = 0; x < numOfBorderPixels; x++)
657  {
658  #if defined(FFTW_AVAILABLE)
659  src[x][0] = src[x][0] - centerRow;
660  src[x][1] = src[x][1] - centerCol;
661  totalRe+= (float)src[x][0];
662  totalIm+= (float)src[x][1];
663  #else
664  src[x].real (src[x].real () - centerRow);
665  src[x].imag (src[x].imag () - centerCol);
666  totalRe+= (float)src[x].real ();
667  totalIm+= (float)src[x].imag ();
668  #endif
669  }
670 
671  kkint32 numOfedgePixels = numOfBorderPixels;
672 
673  #if defined(FFTW_AVAILABLE)
674  fftwf_complex* dest = (fftwf_complex*)fftwf_malloc (sizeof (fftwf_complex) * maxNumOfBorderPoints);
675  fftwf_plan plan = fftwCreateOneDPlan (numOfBorderPixels, src, dest, FFTW_FORWARD, FFTW_ESTIMATE);
676  fftwf_execute (plan);
677  fftwDestroyPlan (plan);
678  plan = NULL;
679  #else
680  KK_DFT1D_Float plan (numOfBorderPixels, true);
681  KK_DFT1D_Float::DftComplexType* dest = new KK_DFT1D_Float::DftComplexType[numOfBorderPixels];
682  plan.Transform (src, dest);
683  #endif
684 
685  kkint32* count = new kkint32 [numOfBuckets];
686 
687  for (x = 0; x < numOfBuckets; x++)
688  {
689  countourFreq[x] = 0.0;
690  count[x] = 0;
691  }
692 
693  float middle = (float)numOfBorderPixels / (float)2.0;
694 
695  float r0 = middle * ( (float) 8.0 / (float)16.0);
696  float r1 = middle * ( (float)10.0 / (float)16.0);
697  float r2 = middle * ( (float)12.0 / (float)16.0);
698  float r3 = middle * ( (float)13.0 / (float)16.0);
699 
700 
701  float deltaX;
702  float mag;
703 
704  kkint32 region = 0;
705 
706  for (x = 1; x < numOfBorderPixels; x++)
707  {
708  //mag = log (sqrt (dest[x].re * dest[x].re + dest[x].im * dest[x].im));
709  //mag = log (sqrt (dest[x].re * dest[x].re + dest[x].im * dest[x].im) + 1.0);
710 
711 
712  #if defined(FFTW_AVAILABLE)
713  mag = (float)sqrt (dest[x][0] * dest[x][0] + dest[x][1] * dest[x][1]);
714  #else
715  mag = (float)sqrt (dest[x].real () * dest[x].real () + dest[x].imag () * dest[x].imag ());
716  #endif
717 
718  deltaX = (float)x - middle;
719 
720  if (fabs (deltaX) < r0)
721  continue;
722 
723  if (deltaX < 0)
724  {
725  // We are on the Left half.
726  deltaX = fabs (deltaX);
727 
728  if (x == 1) region = 0;
729  else if (x == 2) region = 1;
730  else if (x == 3) region = 2;
731  else if (x == 4) region = 3;
732  else if (deltaX >= r3) region = 4;
733  else if (deltaX >= r2) region = 5;
734  else if (deltaX >= r1) region = 6;
735  else region = 7;
736  }
737  else
738  {
739  // We are on Right Half of array
740  if (x == (numOfBorderPixels - 1)) region = 15;
741  else if (x == (numOfBorderPixels - 2)) region = 14;
742  else if (x == (numOfBorderPixels - 3)) region = 13;
743  else if (x == (numOfBorderPixels - 4)) region = 12;
744  else if (deltaX >= r3) region = 11;
745  else if (deltaX >= r2) region = 10;
746  else if (deltaX >= r1) region = 9;
747  else region = 8;
748  }
749 
750  countourFreq[region] = countourFreq[region] + mag;
751  count[region]++;
752  }
753 
754  for (x = 0; x < numOfBuckets; x++)
755  {
756  if (count[x] <= 0)
757  {
758  countourFreq[x] = 0.0;
759  }
760  else
761  {
762  countourFreq[x] = countourFreq[x] / (float)count[x];
763  }
764  }
765 
766  #if defined(FFTW_AVAILABLE)
767  fftwf_free (src); src = NULL;
768  fftwf_free (dest); dest = NULL;
769  #else
770  delete[] src; src = NULL;
771  delete[] dest; dest = NULL;
772  #endif
773 
774  delete[] count;
775  count = NULL;
776 
777 
778  return numOfedgePixels;
779 } /* FollowContour2 */
780 
781 
782 
784 {
785  // startRow and startCol is assumed to come from the left (6)
786 
787  kkint32 startRow = 0;
788  kkint32 startCol = 0;
789 
790  kkint32 scndRow = 0;
791  kkint32 scndCol = 0;
792 
793  kkint32 lastRow = 0;
794  kkint32 lastCol = 0;
795 
796  kkint32 nextRow = 0;
797  kkint32 nextCol = 0;
798 
799  PointListPtr points = new PointList (true);
800 
801  GetFirstPixel (startRow, startCol);
802  if ((startRow < 0) || (startCol < 0) || (PixelCountIn9PixelNeighborhood (startRow, startCol) < 2))
803  {
804  delete points;
805  points = NULL;
806  cout << "Very Bad Starting Point" << std::endl;
807  return NULL;
808  }
809 
810  GetNextPixel (scndRow, scndCol);
811 
812  nextRow = scndRow;
813  nextCol = scndCol;
814 
815  lastRow = startRow;
816  lastCol = startCol;
817 
818  while (true)
819  {
820  lastRow = nextRow;
821  lastCol = nextCol;
822 
823  GetNextPixel (nextRow, nextCol);
824 
825  if ((nextRow == scndRow) && (nextCol == scndCol) &&
826  (lastRow == startRow) && (lastCol == startCol))
827  {
828  break;
829  }
830 
831  points->PushOnBack (new Point (nextRow, nextCol));
832  }
833 
834  return points;
835 } /* GenerateContourList */
836 
837 
838 
839 
840 
841 
843 {
844  static
845  kkint32 curMaskSize = 0;
846 
847  static
848  ComplexDouble** fourierMask = NULL;
849 
850  if (size == curMaskSize)
851  return fourierMask;
852 
853  ComplexDouble N(size, 0);
854  ComplexDouble M(size, 0);
855 
856  kkint32 x;
857 
858  if (fourierMask)
859  {
860  for (x = 0; x < curMaskSize; x++)
861  {
862  delete fourierMask[x];
863  fourierMask[x] = NULL;
864  }
865 
866  delete[] fourierMask;
867  fourierMask = NULL;
868  }
869 
870  fourierMask = new ComplexDouble*[size];
871  for (x = 0; x < size; x++)
872  {
873  fourierMask[x] = new ComplexDouble[size];
874  }
875  curMaskSize = size;
876 
877 
878  ComplexDouble j;
879  ComplexDouble MinusOne (-1, 0);
880  j = sqrt (MinusOne);
881 
882  ComplexDouble Pi (3.14159265359, 0);
883  ComplexDouble One (1.0, 0);
884  ComplexDouble Two (2.0, 0);
885 
886  for (kkint32 m = 0; m < size; m++)
887  {
888  complex<double> mc (m, 0);
889 
890  for (kkint32 k = 0; k < size; k++)
891  {
892  complex<double> kc (k, 0);
893 
894  // double xxx = 2 * 3.14159265359 * (double)k * (double)m / (double)size;
895 
896  // fourierMask[m][k] = exp (MinusOne * j * Two * Pi * kc * mc / M);
897  fourierMask[m][k] = exp (MinusOne * j * Two * Pi * kc * mc / M);
898 
899  double exponetPart = 2.0 * 3.14159265359 * (double)k * (double)m / (double)size;
900  double realPart = cos (exponetPart);
901  double imgPart = -sin (exponetPart);
902 
903  if (realPart != fourierMask[m][k].real ())
904  {
905  continue;
906  }
907 
908  if (imgPart != fourierMask[m][k].imag ())
909  {
910  continue;
911  }
912  }
913  }
914 
915  return fourierMask;
916 } /* GetFourierOneDimMask */
917 
918 
919 
920 
921 ComplexDouble** GetRevFourierOneDimMask (kkint32 size) // For reverse Fourier
922 {
923  static
924  kkint32 curRevMaskSize = 0;
925 
926  static
927  ComplexDouble** revFourierMask = NULL;
928 
929 
930  if (size == curRevMaskSize)
931  return revFourierMask;
932 
933  ComplexDouble N(size, 0);
934  ComplexDouble M(size, 0);
935 
936  kkint32 x;
937 
938 
939  if (revFourierMask)
940  {
941  for (x = 0; x < curRevMaskSize; x++)
942  {
943  delete revFourierMask[x];
944  revFourierMask[x] = NULL;
945  }
946 
947  delete[] revFourierMask;
948  revFourierMask = NULL;
949  }
950 
951 
952  revFourierMask = new ComplexDouble*[size];
953  for (x = 0; x < size; x++)
954  {
955  revFourierMask[x] = new ComplexDouble[size];
956  }
957  curRevMaskSize = size;
958 
959  ComplexDouble j;
960  ComplexDouble MinusOne (-1, 0);
961  ComplexDouble PositiveOne (1, 0);
962  j = sqrt (MinusOne);
963 
964  ComplexDouble Pi (3.1415926, 0);
965  ComplexDouble One (1.0, 0);
966  ComplexDouble Two (2.0, 0);
967 
968  for (kkint32 m = 0; m < size; m++)
969  {
970  complex<double> mc (m, 0);
971 
972  for (kkint32 k = 0; k < size; k++)
973  {
974  complex<double> kc (k, 0);
975 
976  // double xxx = 2 * 3.14159265359 * (double)k * (double)m / (double)size;
977 
978  // fourierMask[m][k] = exp (MinusOne * j * Two * Pi * kc * mc / M);
979  revFourierMask[m][k] = exp (PositiveOne * j * Two * Pi * kc * mc / M);
980 
981 
982  double exponetPart = 2.0 * 3.14159265359 * (double)k * (double)m / (double)size;
983  double realPart = cos (exponetPart);
984  double imgPart = sin (exponetPart);
985 
986  if (realPart != revFourierMask[m][k].real ())
987  {
988  continue;
989  }
990 
991  if (imgPart != revFourierMask[m][k].imag ())
992  {
993  continue;
994  }
995  }
996  }
997 
998  return revFourierMask;
999 } /* GetRevFourierOneDimMask */
1000 
1001 
1002 
1003 
1005  float* countourFreq,
1006  bool& successful
1007  )
1008 {
1009  // startRow and startCol is assumed to come from the left (6)
1010  // kkint32 numOfBuckets = 8;
1011  kkint32 startRow;
1012  kkint32 startCol;
1013 
1014  kkint32 scndRow;
1015  kkint32 scndCol;
1016 
1017  kkint32 lastRow;
1018  kkint32 lastCol;
1019 
1020  kkint32 nextRow;
1021  kkint32 nextCol;
1022 
1023  kkint32 numOfBorderPixels = 0;
1024 
1025  kkint32 x;
1026 
1027  successful = true;
1028 
1029  PointListPtr points = new PointList (true);
1030 
1031  GetFirstPixel (startRow, startCol);
1032  if ((startRow < 0) || (startCol < 0) || (PixelCountIn9PixelNeighborhood (startRow, startCol) < 2))
1033  {
1034  cout << "Vary Bad Starting Point" << std::endl;
1035  successful = false;
1036  return 0;
1037  }
1038 
1039  GetNextPixel (scndRow, scndCol);
1040 
1041  nextRow = scndRow;
1042  nextCol = scndCol;
1043 
1044  lastRow = startRow;
1045  lastCol = startCol;
1046 
1047 
1048  while (true)
1049  {
1050  lastRow = nextRow;
1051  lastCol = nextCol;
1052 
1053  GetNextPixel (nextRow, nextCol);
1054 
1055  if ((nextRow == scndRow) && (nextCol == scndCol) &&
1056  (lastRow == startRow) && (lastCol == startCol))
1057  {
1058  break;
1059  }
1060 
1061  points->PushOnBack (new Point (nextRow, nextCol));
1062 
1063  numOfBorderPixels++;
1064  }
1065 
1066 
1067  #if defined(FFTW_AVAILABLE)
1068  fftwf_complex* src = (fftwf_complex*)fftwf_malloc (sizeof (fftwf_complex) * numOfBuckets);
1069  #else
1070  KK_DFT1D_Float::DftComplexType* src = new KK_DFT1D_Float::DftComplexType[numOfBuckets];
1071  #endif
1072 
1073  kkint32 totalRow = 0;
1074  kkint32 totalCol = 0;
1075 
1076  for (x = 0; x < numOfBuckets; x++)
1077  {
1078  kkint32 borderPixelIdx = (kkint32)(((double)x * (double)numOfBorderPixels) / (double)numOfBuckets);
1079 
1080  Point& point = (*points)[borderPixelIdx];
1081 
1082  #if defined(FFTW_AVAILABLE)
1083  src[x][0] = (float)point.Row ();
1084  src[x][1] = (float)point.Col ();
1085  #else
1086  src[x].real ((float)point.Row ());
1087  src[x].imag ((float)point.Col ());
1088  #endif
1089 
1090  totalRow += point.Row ();
1091  totalCol += point.Col ();
1092  }
1093 
1094  delete points;
1095  points = NULL;
1096 
1097  float centerRow = (float)totalRow / (float)numOfBuckets;
1098  float centerCol = (float)totalCol / (float)numOfBuckets;
1099 
1100  for (x = 0; x < numOfBuckets; x++)
1101  {
1102  #if defined(FFTW_AVAILABLE)
1103  src[x][0] = src[x][0] - centerRow;
1104  src[x][1] = src[x][1] - centerCol;
1105  #else
1106  src[x].real (src[x].real () - centerRow);
1107  src[x].imag (src[x].imag () - centerCol);
1108  #endif
1109  }
1110 
1111  #if defined(FFTW_AVAILABLE)
1112  fftwf_complex* dest = (fftwf_complex*)fftwf_malloc (sizeof (fftwf_complex) * numOfBuckets);
1113  fftwf_plan plan = fftwCreateOneDPlan (numOfBuckets, src, dest, FFTW_FORWARD, FFTW_ESTIMATE);
1114  fftwf_execute (plan);
1115  fftwDestroyPlan (plan);
1116  plan = NULL;
1117  #else
1118  KK_DFT1D_Float plan (numOfBuckets, true);
1119  KK_DFT1D_Float::DftComplexType* dest = new KK_DFT1D_Float::DftComplexType[numOfBuckets];
1120  plan.Transform (src, dest);
1121  #endif
1122 
1123  for (x = 0; x < numOfBuckets; x++)
1124  {
1125  #if defined(FFTW_AVAILABLE)
1126  float real = dest[x][0];
1127  float imag = dest[x][1];
1128  #else
1129  float real = dest[x].real ();
1130  float imag = dest[x].imag ();
1131  #endif
1132 
1133  countourFreq[x] = (float)(sqrt (real * real + imag * imag));
1134  }
1135 
1136  #if defined(FFTW_AVAILABLE)
1137  fftwf_free (src);
1138  fftwf_free (dest);
1139  #else
1140  delete[] src; src = NULL;
1141  delete[] dest; dest = NULL;
1142  #endif
1143 
1144  return numOfBorderPixels;
1145 } /* CreateFourierDescriptorBySampling */
1146 
1147 
1148 
1149 
1151  float pointCol,
1152  kkint32 numOfBuckets,
1153  kkint32* buckets,
1154  float& minDistance,
1155  float& maxDistance,
1156  kkint32& numOfEdgePixels
1157  )
1158 {
1159  PointListPtr points = GenerateContourList ();
1160 
1161  numOfEdgePixels = points->QueueSize ();
1162 
1163  kkint32 x;
1164 
1165  minDistance = FloatMax;
1166  maxDistance = FloatMin;
1167 
1168  for (x = 0; x < numOfBuckets; x++)
1169  buckets[x] = 0;
1170 
1171  float* distances = new float[points->QueueSize ()];
1172 
1173  for (x = 0; x < points->QueueSize (); x++)
1174  {
1175  Point& point = (*points)[x];
1176 
1177  float deltaCol = (float)point.Col () - pointCol;
1178  float deltaRow = (float)point.Row () - pointRow;
1179 
1180  float distance = (float)sqrt (deltaCol * deltaCol + deltaRow * deltaRow);
1181 
1182  distances[x] = distance;
1183 
1184  minDistance = Min (minDistance, distance);
1185  maxDistance = Max (maxDistance, distance);
1186  }
1187 
1188  float bucketSize = (maxDistance - minDistance) / numOfBuckets;
1189 
1190  if (bucketSize == 0.0)
1191  {
1192  buckets[numOfBuckets / 2] = points->QueueSize ();
1193  }
1194  else
1195  {
1196  for (x = 0; x < points->QueueSize (); x++)
1197  {
1198  kkint32 bucketIDX = (kkint32)((distances[x] - minDistance) / bucketSize);
1199  buckets[bucketIDX]++;
1200  }
1201  }
1202 
1203  delete points; points = NULL;
1204  delete[] distances; distances = NULL;
1205 
1206  return;
1207 } /* HistogramDistanceFromAPoint */
1208 
1209 
1210 
1212 {
1213  //ComplexDouble* src = new ComplexDouble[points.QueueSize ()];
1214 
1215  #if defined(FFTW_AVAILABLE)
1216  fftw_complex* src = (fftw_complex*)fftw_malloc (sizeof (fftw_complex) * points.QueueSize ());
1217  #else
1218  KK_DFT1D_Double::DftComplexType* src = new KK_DFT1D_Double::DftComplexType[points.QueueSize()];
1219  #endif
1220 
1221 
1222  kkint32 totalRow = 0;
1223  kkint32 totalCol = 0;
1224 
1225  kkint32 x = 0;
1226 
1227  for (x = 0; x < points.QueueSize (); x++)
1228  {
1229  Point& point (points[x]);
1230 
1231  //src[x] = ComplexDouble ((double)point.Row (), (double)point.Col ());
1232 
1233  #if defined(FFTW_AVAILABLE)
1234  src[x][0] = (double)point.Row ();
1235  src[x][1] = (double)point.Col ();
1236  #else
1237  src[x].real ((double)point.Row ());
1238  src[x].imag ((double)point.Col ());
1239  #endif
1240 
1241  totalRow += point.Row ();
1242  totalCol += point.Col ();
1243  }
1244 
1245  //double centerRow = (double)totalRow / (double)(points.QueueSize ());
1246  //double centerCol = (double)totalCol / (double)(points.QueueSize ());
1247 
1248 // for (x = 0; x < points.QueueSize (); x++)
1249 // {
1250 // src[x] = ComplexDouble ((double)(src[x].real () - centerRow), (double)(src[x].imag () - centerCol));
1251 // }
1252 //
1253 
1254 
1255  #if defined(FFTW_AVAILABLE)
1256  fftw_complex* destFFTW = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * points.QueueSize());
1257  fftw_plan plan;
1258  plan = fftw_plan_dft_1d (points.QueueSize (), src, destFFTW, FFTW_FORWARD, FFTW_ESTIMATE);
1259  fftw_execute (plan);
1260  fftw_destroy_plan (plan);
1261 #else
1262  KK_DFT1D_Double::DftComplexType* destFFTW = new KK_DFT1D_Double::DftComplexType[points.QueueSize()];
1263  KK_DFT1D_Double plan(points.QueueSize(), true);
1264  plan.Transform (src, destFFTW);
1265  #endif
1266 
1267  vector<ComplexDouble> dest;
1268 
1269  for (kkint32 l = 0; l < points.QueueSize (); l++)
1270  {
1271  #if defined(FFTW_AVAILABLE)
1272  dest.push_back (ComplexDouble (destFFTW[l][0] / (double)(points.QueueSize ()), destFFTW[l][1] / (double)(points.QueueSize ())));
1273  #else
1274  dest.push_back (ComplexDouble (destFFTW[l].real () / (double)(points.QueueSize ()), destFFTW[l].imag () / (double)(points.QueueSize ())));
1275  #endif
1276  }
1277  delete[] destFFTW;
1278  delete[] src;
1279 
1280  return dest;
1281 } /* CreateFourierFromPointList */
1282 
1283 
1284 
1285 
1286 
1287 PointListPtr ContourFollower::CreatePointListFromFourier (vector<ComplexDouble> fourier,
1288  PointList& origPointList
1289  )
1290 {
1291  kkint32 minRow, maxRow, minCol, maxCol;
1292  origPointList.BoxCoordinites (minRow, minCol, maxRow, maxCol);
1293 
1294  PointListPtr points = new PointList (true);
1295 
1296  size_t numOfEdgePixels = origPointList.size ();
1297 
1298  #if defined(FFTW_AVAILABLE)
1299  fftw_complex* src = (fftw_complex*)fftw_malloc (sizeof (fftw_complex) * numOfEdgePixels);
1300  #else
1301  KK_DFT1D_Double::DftComplexType* src = new KK_DFT1D_Double::DftComplexType[numOfEdgePixels];
1302  #endif
1303 
1304  for (kkint32 l = 0; l < (kkint32)fourier.size (); l++)
1305  {
1306  #if defined(FFTW_AVAILABLE)
1307  src[l][0] = fourier[l].real ();
1308  src[l][1] = fourier[l].imag ();
1309  #else
1310  src[l].real (fourier[l].real());
1311  src[l].imag (fourier[l].imag ());
1312  #endif
1313  }
1314 
1315  #if defined(FFTW_AVAILABLE)
1316  fftw_complex* destFFTW = (fftw_complex*)fftw_malloc (sizeof (fftw_complex) * numOfEdgePixels);
1317  fftw_plan plan;
1318  plan = fftw_plan_dft_1d (numOfEdgePixels, src, destFFTW, FFTW_FORWARD, FFTW_ESTIMATE);
1319  fftw_execute (plan);
1320  fftw_destroy_plan (plan);
1321  #else
1322  KK_DFT1D_Double::DftComplexType* destFFTW = new KK_DFT1D_Double::DftComplexType[numOfEdgePixels];
1323  KK_DFT1D_Double plan(numOfEdgePixels, false);
1324  plan.Transform(src, destFFTW);
1325  #endif
1326 
1327 
1328  kkint32 largestRow = -1;
1329  kkint32 largestCol = -1;
1330  kkint32 smallestRow = 999999;
1331  kkint32 smallestCol = 999999;
1332 
1333  for (kkint32 l = 0; l < (kkint32)fourier.size (); l++)
1334  {
1335 
1336  #if defined(FFTW_AVAILABLE)
1337  double realPart = (double)destFFTW[l][0];
1338  double imagPart = (double)destFFTW[l][1];
1339  ComplexDouble z (realPart, imagPart);
1340  #else
1341  double realPart = (double)destFFTW[l].real ();
1342  double imagPart = (double)destFFTW[l].imag ();
1343  ComplexDouble z (realPart, imagPart);
1344  #endif
1345 
1346 
1347  // z = z / (double)(fourier.size ());
1348 
1349  kkint32 row = (kkint32)(z.real () + 0.5);
1350  if (row > largestRow)
1351  largestRow = row;
1352  if (row < smallestRow)
1353  smallestRow = row;
1354 
1355  kkint32 col = (kkint32)(z.imag () + 0.5);
1356  if (col > largestCol)
1357  largestCol = col;
1358  if (col < smallestCol)
1359  smallestCol = col;
1360 
1361  PointPtr p = new Point (row, col);
1362  points->PushOnBack (p);
1363  }
1364 
1365  delete[] src;
1366  delete[] destFFTW;
1367 
1368  return points;
1369 } /* CreatePointListFromFourier */
kkint32 FollowContour2(float *countourFreq, bool &successful)
__int32 kkint32
Definition: KKBaseTypes.h:88
float FloatMin
Definition: KKBaseTypes.cpp:21
std::vector< ComplexDouble > CreateFourierFromPointList(const PointList &points)
kkint32 Col() const
Definition: Point.h:40
float FloatMax
Definition: KKBaseTypes.cpp:20
kkint32 Row() const
Definition: Point.h:39
A class that is used by to represent a single image in memory.
Definition: Raster.h:108
kkint32 CreateFourierDescriptorBySampling(kkint32 numOfBuckets, float *countourFreq, bool &successful)
PointList(bool _owner)
Definition: Point.cpp:105
kkint32 col
Definition: Raster.h:1381
float CalcMagnitude(KK_DFT1D_Float::DftComplexType *dest, kkint32 index)
Used by Raster class and MorphOp derived classes to denote a single pixel location in Raster image...
Definition: Point.h:20
kkint32 Height() const
Definition: Raster.h:319
KK_DFT1D< double > KK_DFT1D_Double
Definition: kku_fftw.h:173
std::complex< double > ComplexDouble
KKTHread * KKTHreadPtr
KK_DFT1D< float > KK_DFT1D_Float
Definition: kku_fftw.h:172
ComplexDouble ** GetRevFourierOneDimMask(kkint32 size)
unsigned char uchar
Unsigned character.
Definition: KKBaseTypes.h:77
ComplexDouble ** GetFourierOneDimMask(kkint32 size)
uchar BackgroundPixelTH() const
Definition: Raster.h:335
Point(kkint32 _row, kkint32 _col)
Definition: Point.cpp:43
void HistogramDistanceFromAPointOfEdge(float pointRow, float pointCol, kkint32 numOfBuckets, kkint32 *buckets, float &minDistance, float &maxDistance, kkint32 &numOfEdgePixels)
kkint32 Width() const
Definition: Raster.h:324
ContourFollower(Raster &_raster, RunLog &_log)
Used for logging messages.
Definition: RunLog.h:49
Container object used to maintaining a list of pixel locations.
Definition: Point.h:75
kkint32 FollowContour(float *countourFreq, float fourierDescriptors[15], kkint32 totalPixels, bool &successful)
void BoxCoordinites(kkint32 &minRow, kkint32 &minCol, kkint32 &maxRow, kkint32 &maxCol)
Definition: Point.cpp:111
PointListPtr GenerateContourList()
kkint32 row
Definition: Raster.h:1380
const MovDir movements[8]
uchar ** Rows() const
Definition: Raster.h:321
PointListPtr CreatePointListFromFourier(std::vector< ComplexDouble > fourier, PointList &origPointList)