Basic Image AlgorithmS Library  2.8.0
 All Classes Namespaces Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
CornerDetectorGradient.cpp
1 /*
2 This file is part of the BIAS library (Basic ImageAlgorithmS).
3 
4 Copyright (C) 2003-2009 (see file CONTACT for details)
5  Multimediale Systeme der Informationsverarbeitung Institut fuer Informatik
6  Christian-Albrechts-Universitaet Kiel
7 
8 
9 BIAS is free software; you can redistribute it and/or modify
10 it under the terms of the GNU Lesser General Public License as published by
11 the Free Software Foundation; either version 2.1 of the License, or
12 (at your option) any later version.
13 
14 BIAS is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 GNU Lesser General Public License for more details.
18 
19 You should have received a copy of the GNU Lesser General Public License
20 along with BIAS; if not, write to the Free Software
21 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
22 */
23 
24 #ifdef WIN32
25 # include "Base/Common/W32Compat.hh"
26 #endif
27 
28 #include <algorithm>
29 #include "CornerDetectorGradient.hh"
30 #include <Base/Image/ImageConvert.hh>
31 
32 #ifdef BIAS_DEBUG
33 #include <Base/Image/ImageIO.hh>
34 #endif
35 
36  //#define TIME_MEASURE
37 #ifdef TIME_MEASURE
38 #include <Base/Debug/TimeMeasure.hh>
39 #endif
40 
41 #include <Base/Common/BIASpragma.hh>
42 
43 using namespace BIAS;
44 using namespace std;
45 
46 //////////////////////////////////////////////////////////////////////////
47 // implementation
48 //////////////////////////////////////////////////////////////////////////
49 
50 
51 template <class StorageType, class CalculationType>
52 CornerDetectorGradient<StorageType, CalculationType>::
53 CornerDetectorGradient()
54  : CornerDetectorBase<StorageType>(),
55  _st(),
56  _gauss()
57 {
58  _FeatureMap=NULL;
59  _MinDistance = 10;
60  _MinCornerness = 1.0;
61  _DoRefinement = false;
62 }
63 template <class StorageType, class CalculationType>
64 CornerDetectorGradient<StorageType, CalculationType>::~CornerDetectorGradient()
65 {
66  _DeleteInternalMem();
67 }
68 
69 template <class StorageType, class CalculationType>
70 int CornerDetectorGradient<StorageType, CalculationType>::
71 Detect(const Image<StorageType>& image, std::vector<HomgPoint2D>& p,
72  std::vector<QUAL>& quality, std::vector<Matrix2x2<double> > *cov)
73 {
74  Image<StorageType> fim;
75  _gauss.Filter(image, fim);
76  _gauss.SetBorderHandling(BIAS::FilterBase<StorageType, CalculationType>::TBH_same);
77 
78 
79 #ifdef BIAS_DEBUG
80  if (this->DebugLevelIsSet(D_CD_WRITE_DEBUG_IM)){
81  Image<StorageType> im2=image;
82  //ImageIO::Save("float-input.mip", im2);
83  //ImageIO::Save("gauss-filtered.mip", fim);
84  ImageIO::Save("input.mip", image);
85  ImageIO::Save("gauss-filtered.mip", fim);
86  }
87 #endif
88 
89  Image<CalculationType> *kim = NULL;
90  if (typeid(StorageType)==typeid(CalculationType)) {
91  kim = (Image<CalculationType> *)&fim;
92  } else {
93  kim = new Image<CalculationType>(fim.GetWidth(), fim.GetHeight(),
94  fim.GetChannelCount());
95  ImageConvert::ConvertST(fim, *kim, kim->GetStorageType());
96  }
97  int sres=this->_st.CalcStructureTensor(*kim, this->_sgxx, this->_sgxy,
98  this->_sgyy);
99  if (typeid(StorageType)!=typeid(CalculationType)) delete kim;
100 #ifdef BIAS_DEBUG
101  else if (this->DebugLevelIsSet(D_CD_WRITE_DEBUG_IM)) {
102  //ImageIO::Save("structure-tensor.mip", *kim);
103  ImageIO::Save("structure-tensor-input.mip", *kim);
104  }
105 #endif
106 
107  if (this->_sgxx.GetWidth()>this->_Cornerness.GetWidth() ||
108  this->_sgxx.GetHeight()>this->_Cornerness.GetHeight()){
109  _AllocInternalMem(this->_sgxx.GetWidth(), this->_sgxx.GetHeight());
110  }
111 
112  int mres = _ComputeCornerness(this->_Cornerness);
113 #ifdef BIAS_DEBUG
114  if (this->DebugLevelIsSet(D_CD_WRITE_DEBUG_IM)) {
115  //ImageIO::Save("cornerness.mip", this->_Cornerness);
116  ImageIO::Save("cornerness.mip", this->_Cornerness);
117  }
118 #endif
119  this->_ZeroFeatureMap();
120  // do we have old corners where no new ones are desired ?
121  if (p.size()>0) this->_FillFeatureMap(p);
122  int eres = this->_EnforceMinimumDistance(p, quality);
123 
124  // refine all selected corners
125  if (_DoRefinement) _RefineCorners(p, quality);
126 
127  // compute the covariances
128  if (cov) {
129  GetCov_(p, *cov);
130  }
131 
132  return (sres==0 && mres==0 && eres==0)?(0):-1;
133 }
134 
135 template <class StorageType, class CalculationType>
136 void CornerDetectorGradient<StorageType, CalculationType>::
137 GetCov_(const vector<HomgPoint2D>& p, vector<Matrix2x2<double> >& cov,
138  const double noise) const
139 {
140  const unsigned int num = p.size();
141  int x, y, ires;
142  const double square_noise = noise * noise;
143  Matrix2x2<double> c;
144  cov.resize(num);
145  for (unsigned int i = 0; i < num; i++){
146  x = (int)rint(p[i][0]);
147  y = (int)rint(p[i][1]);
148 #ifdef BIAS_DEBUG
149  if (!(x > 0 && y >= 0 &&
150  x < (int)_sgxx.GetWidth() && y < (int)_sgxx.GetHeight() &&
151  x < (int)_sgxy.GetWidth() && y < (int)_sgxy.GetHeight() &&
152  x < (int)_sgyy.GetWidth() && y < (int)_sgyy.GetHeight())) {
153  BIASERR("Invalid feature point position given (" << x
154  << ", " << y << "), setting covariance to max.!");
155  //BIASABORT;
156  cov[i][0][0] = cov[i][1][1] = DBL_MAX;
157  cov[i][1][0] = cov[i][0][1] = 0.0;
158  continue;
159  }
160 #endif
161  c[0][0] = _sgxx.GetImageDataArray()[y][x];
162  c[1][1] = _sgyy.GetImageDataArray()[y][x];
163  c[1][0] = c[0][1] = _sgxy.GetImageDataArray()[y][x];
164  ires = c.Invert(cov[i]);
165  if (ires != 0) {
166  BIASERR("Inversion of covariance matrix " << c << " failed (result of "
167  "Invert = " << ires << "), setting covariance matrix to max.!");
168  //BIASABORT;
169  cov[i][0][0] = cov[i][1][1] = DBL_MAX;
170  cov[i][1][0] = cov[i][0][1] = 0.0;
171  } else {
172  cov[i] *= square_noise;
173  }
174  }
175 }
176 
177 template <class StorageType, class CalculationType>
178 int CornerDetectorGradient<StorageType, CalculationType>::
179 Detect(const Image<CalculationType>& gradx,
180  const Image<CalculationType>& grady,
181  std::vector<HomgPoint2D>& p,std::vector<QUAL>& quality,
182  std::vector<Matrix2x2<double> > *cov)
183 {
184 #ifdef TIME_MEASURE
185  TimeMeasure timer2;
186  timer2.Start();
187 #endif
188 
189 #ifdef BIAS_DEBUG
190  if (this->DebugLevelIsSet(D_CD_WRITE_DEBUG_IM)){
191  Image<CalculationType> im2=gradx;
192  //ImageIO::Save("cd-gradx.mip", im2);
193  ImageIO::Save("cd-gradx.mip", im2);
194  im2=grady;
195  //ImageIO::Save("cd-grady.mip", im2);
196  ImageIO::Save("cd-grady.mip", im2);
197  }
198  BIASASSERT(gradx.SamePixelAndChannelCount(grady));
199 #endif
200 
201 
202 #ifdef TIME_MEASURE
203  timer2.Stop();
204  cerr << "allocating took "<<timer2.GetRealTime()<<" us\n";
205  timer2.Reset();
206  timer2.Start();
207 #endif
208 
209  int sres = 0;
210  // create instances of gradienttype images, so we can compare against
211  // storagetype
212  sres=this->_st.CalcStructureTensor(gradx, grady, this->_sgxx,
213  this->_sgxy, this->_sgyy);
214 
215 #ifdef BIAS_DEBUG
216  if (this->DebugLevelIsSet(D_CD_WRITE_DEBUG_IM)){
217  ImageIO::Save("cd-sgxx.mip", _sgxx);
218  ImageIO::Save("cd-sgxy.mip", _sgxy);
219  ImageIO::Save("cd-sgyy.mip", _sgyy);
220  }
221 #endif
222 
223 #ifdef TIME_MEASURE
224  timer2.Stop();
225  cerr << "CalcStructureTensor took "<<timer2.GetRealTime()<<" us\n";
226  timer2.Reset();
227  timer2.Start();
228 #endif
229 
230  if (this->_sgxx.GetWidth()>this->_Cornerness.GetWidth() ||
231  this->_sgxx.GetHeight()>this->_Cornerness.GetHeight()){
232  _AllocInternalMem(this->_sgxx.GetWidth(), this->_sgxx.GetHeight());
233  }
234  int mres = _ComputeCornerness(this->_Cornerness);
235 
236 #ifdef BIAS_DEBUG
237  if (this->DebugLevelIsSet(D_CD_WRITE_DEBUG_IM)){
238  //ImageIO::Save("cd-cornerness.mip", _Cornerness);
239  ImageIO::Save("cd-cornerness.mip", _Cornerness);
240  }
241 #endif
242 
243 #ifdef TIME_MEASURE
244  timer2.Stop();
245  cerr << "_ComputeCornerness took "<<timer2.GetRealTime()<<" us\n";
246  timer2.Reset();
247  timer2.Start();
248 #endif
249 
250  this->_ZeroFeatureMap();
251  if (p.size()>0) this->_FillFeatureMap(p);
252  int eres = this->_EnforceMinimumDistance(p, quality);
253 
254 #ifdef TIME_MEASURE
255  timer2.Stop();
256  cerr << "EnforceMinimumDistance took "<<timer2.GetRealTime()<<" us\n";
257  timer2.Reset();
258  timer2.Start();
259 #endif
260 
261  // refine all selected corners
262  if (_DoRefinement) _RefineCorners(p, quality);
263 
264  // compute the covariances
265  if (cov) {
266  GetCov_(p, *cov);
267  }
268 
269  return (sres==0 && mres==0 && eres==0)?(0):-1;
270 }
271 
272 template <class StorageType, class CalculationType>
273 int CornerDetectorGradient<StorageType, CalculationType>::
274 Cornerness(const Image<CalculationType>& gradx,
275  const Image<CalculationType>& grady,
276  Image<CalculationType>& cornerness)
277 {
278  int ires, meres;
279  if (!gradx.SamePixelAndChannelCount(this->_sgxx)){
280  _AllocInternalMem(gradx.GetWidth(), gradx.GetHeight());
281  }
282 
283  if ((ires=this->_st.CalcStructureTensor(gradx, grady, this->_sgxx,
284  this->_sgxy, this->_sgyy))!=0){
285  BIASERR("error calculating structure tensor");
286  }
287 
288  if (!cornerness.SamePixelAndChannelCount(gradx)){
289  if (!cornerness.IsEmpty()) cornerness.Release();
290  cornerness.Init(gradx.GetWidth(), gradx.GetHeight(), 1);
291  }
292 
293 
294  meres = _ComputeCornerness(cornerness);
295 
296  return (ires==0 && meres==0)?0:-1;
297 }
298 
299 template <class StorageType, class CalculationType>
300 int CornerDetectorGradient<StorageType, CalculationType>::
301 DetectFromCornerness(const Image<CalculationType>& cornerness,
302  std::vector<HomgPoint2D>& p, std::vector<QUAL>& quality)
303 {
304  // fill feature map from cornerness
305  _ExtractLocalMaxima(cornerness);
306 
307  // now do non maximum supression
308  this->_ZeroFeatureMap();
309  if (p.size()>0) this->_FillFeatureMap(p);
310  int res = this->_EnforceMinimumDistance(p, quality);
311 
312  // refine all selected corners
313  if (_DoRefinement) _RefineCorners(p, quality);
314 
315  return res;
316 }
317 
318 
319 template <class StorageType, class CalculationType>
320 int CornerDetectorGradient<StorageType, CalculationType>::
321 DetectFromCornerness(const Image<CalculationType> &cornerness,
322  const std::vector<Vector2<int> > &region,
323  std::vector<HomgPoint2D> &p,
324  std::vector<QUAL> &quality)
325 {
326 #ifdef BIAS_DEBUG
327  BIASASSERT( p.size() == quality.size() );
328 #endif
329 
330  const CalculationType **imagedatacorn = cornerness.GetImageDataArray();
331  const double mincorn = (double)this->_MinCornerness;
332 
333  // setup feature list
334  this->_FeatList.clear();
335  this->_FeatList.reserve(region.size());
336 
337  // fill feature list with features of min cornerness
338  Feat fp;
339  unsigned i;
340  for (i = 0; i < region.size(); i++) {
341  const int x = region[i][0];
342  const int y = region[i][1];
343  const double corn = imagedatacorn[y][x];
344 
345  if ( corn > mincorn ) {
346  fp.x = x;
347  fp.y = y;
348  fp.val = corn;
349 
350  this->_FeatList.push_back(fp);
351  }
352  }
353 
354  // init feature map, mask out given features by 'p'
355  this->_ZeroFeatureMap();
356  if ( p.size() > 0)
357  this->_FillFeatureMap(p);
358 
359  int res = this->_EnforceMinimumDistance(p, quality);
360 
361  // refine all selected corners
362  if (_DoRefinement) _RefineCorners(p, quality);
363 
364  return res;
365 }
366 
367 //////////////////////////////////////////////////////////////////////////
368 // protected functions
369 //////////////////////////////////////////////////////////////////////////
370 
371 template <class StorageType, class CalculationType>
372 int CornerDetectorGradient<StorageType, CalculationType>::
373 _EnforceMinimumDistance(vector<HomgPoint2D>& points, vector<QUAL>& qual)
374 {
375 #ifdef TIME_MEASURE
376  TimeMeasure timer;
377  timer.Start();
378 #endif
379 
380  //for (int i=0; i<50; i++) cout << i<< ": "<< _FeatList[i] << endl;
381  //int maxnum=(_InternalMemWidth*_InternalMemHeight)/8;
382  //partial_sort(_FeatList.begin(), _FeatList.begin()+maxnum, _FeatList.end());
383  sort(_FeatList.begin(), _FeatList.end());
384  // for (int i=0; i<50; i++) cout << i<< ": "<< _FeatList[i] << endl;
385 
386 #ifdef TIME_MEASURE
387  timer.Stop();
388  cerr << "Quicksort of "<<_FeatList.size()<<" took "<<timer.GetRealTime()<<" us\n";
389  timer.Reset();
390  timer.Start();
391 #endif
392 
393  // clear points and qual
394  points.clear();
395  qual.clear();
396 
397  int minx, maxx, miny, maxy, mx, my;
398  register int x, y;
399  register double mineig;
400  //register vector<class Feat>::iterator p=_FeatList.begin();
401  const int md=_MinDistance-1;
402  const int width=_Cornerness.GetWidth();
403  const int height=_Cornerness.GetHeight();
404  const int mwidth=width-md;
405  const int mheight=height-md;
406  HomgPoint2D hp;
407 
408  // 1 point in every md by md pixel patch
409  int num=2000;
410  if (md>0) num=width*height/(md*md);
411  points.reserve(num);
412  qual.reserve(num);
413 
414  //BIASASSERT(md>=0);
415 
416  // fill _FeatureMap and points
417  vector<class Feat>::iterator it;
418  for (it=_FeatList.begin(); it!=_FeatList.end(); it++)
419  {
420  // check if we have found enough corners:
421  if (_NumberFeaturesSoFar >= this->_MaxNum && this->_MaxNum > 0)
422  break;
423  x=it->x;
424  y=it->y;
425  mineig=it->val;
426  if (!_FeatureMap[y][x]) {
427  // fill feature map
428  miny=(y>md)?(y-md):(0);
429  maxy=(y>=mheight)?(height-1):(y+md);
430  minx=(x>md)?(x-md):(0);
431  maxx=(x>=mwidth)?(width-1):(x+md);
432  for (my=miny; my<=maxy; my++){
433  for (mx=minx; mx<=maxx; mx++){
434  _FeatureMap[my][mx]=true;
435  }
436  }
437  // add point to vector
438  hp.Set((double)(x), (double)(y), 1.0);
439 #ifdef BIAS_DEBUG
440  BIASASSERT(!BIAS_ISNAN(hp[0]) && !BIAS_ISINF(hp[0]));
441  BIASASSERT(!BIAS_ISNAN(hp[1]) && !BIAS_ISINF(hp[1]));
442  if (this->DebugLevelIsSet(D_CD_FEATURES) && points.size()<15)
443  cerr << points.size() << " : ("<<x<<", "<<y<<" ) qual "<<mineig<<endl;
444 #endif
445  points.push_back(hp);
446  qual.push_back((QUAL)mineig);
447  _NumberFeaturesSoFar++;
448  }
449  }
450 
451 #ifdef BIAS_DEBUG
452  for (unsigned i=0; i<points.size(); i++){
453  BIASASSERT(!BIAS_ISNAN(points[i][0]) && !BIAS_ISINF(points[i][0]));
454  BIASASSERT(!BIAS_ISNAN(points[i][1]) && !BIAS_ISINF(points[i][1]));
455  }
456 #endif
457 
458 #ifdef TIME_MEASURE
459  timer.Stop();
460  cerr << "EnforceMinimumDistance without Quicksort took "
461  <<timer.GetRealTime()<<" us\n";
462 #endif
463  return 0;
464 }
465 
466 
467 template <class StorageType, class CalculationType>
468 bool CornerDetectorGradient<StorageType, CalculationType>::
469 _RefineCornerPosition(int col, int row,
470  double &x, double &y, QUAL& NewCornerness) {
471  CalculationType **cornerness = _Cornerness.GetImageDataArray();
472  // cornerness values of 8-neighbourhood
473  const register CalculationType& DR_P11 = cornerness[row-1][col-1];
474  const register CalculationType& DR_P12 = cornerness[row-1][col ];
475  const register CalculationType& DR_P13 = cornerness[row-1][col+1];
476  const register CalculationType& DR_P21 = cornerness[row ][col-1];
477  const register CalculationType& DR_P22 = cornerness[row ][col ];
478  const register CalculationType& DR_P23 = cornerness[row ][col+1];
479  const register CalculationType& DR_P31 = cornerness[row+1][col-1];
480  const register CalculationType& DR_P32 = cornerness[row+1][col ];
481  const register CalculationType& DR_P33 = cornerness[row+1][col+1];
482  NewCornerness = (float)DR_P22;
483 
484  // Create weighting coefficients for the quadratic fit. They are obtained by
485  // applying these masks to the 3x3 window centered on pixel (row,col) :
486  float b,c,d,e,f;
487 
488  // 1/9 * [ -1 +2 -1 ]
489  // [ +2 +5 +2 ]
490  // [ -1 +2 -1 ]
491  // dont need "a" here
492 
493  // 1/6 * [ -1 0 +1 ]
494  // [ -1 0 +1 ]
495  // [ -1 0 +1 ]
496  b = (-DR_P11 + DR_P13 - DR_P21 +
497  DR_P23 - DR_P31 + DR_P33)/6.0f;
498 
499  // 1/6 * [ -1 -1 -1 ]
500  // [ 0 0 0 ]
501  // [ +1 +1 +1 ]
502  c = (-DR_P11 - DR_P12 - DR_P13 +
503  DR_P31 + DR_P32 + DR_P33)/6.0f;
504 
505  // 1/3 * [ +1 -2 +1 ]
506  // [ +1 -2 +1 ]
507  // [ +1 -2 +1 ]
508  d = (DR_P11 + DR_P13 +
509  DR_P21 + DR_P23 +
510  DR_P31 + DR_P33
511  -2.0f*(DR_P12+DR_P22+DR_P32))/3.0f;
512 
513  // 1/4 * [ +1 0 -1 ]
514  // [ 0 0 0 ]
515  // [ -1 0 +1 ]
516  e = (DR_P11 - DR_P31 -
517  DR_P13 + DR_P33)/4.0f;
518 
519  // 1/3 * [ +1 +1 +1 ]
520  // [ -2 -2 -2 ]
521  // [ +1 +1 +1 ]
522  f = ( DR_P11 + DR_P12 + DR_P13
523  -2.0f*(DR_P21+DR_P22+DR_P23)
524  +DR_P31 + DR_P32 + DR_P33)/3.0f;
525 
526  // Next find the extremum of the fitted surface by setting its
527  // partial derivatives to zero. We need to solve the following linear system:
528  //
529  // [ d e ] [ off_x ] + [ b ] = [ 0 ] (dS/dx = 0)
530  // [ e f ] [ off_y ] [ c ] [ 0 ] (dS/dy = 0)
531  //
532  // This implies that the fitted surface is
533  // S(x,y) = a + b x + c y + 1/2 d x^2 + e x y + 1/2 f y^2
534  // where we don't need to know the value of a.
535 
536  float det = d*f - e*e;
537  if (det > 0.0) {
538  float off_x = (c*e - b*f) / det;
539  float off_y = (b*e - c*d) / det;
540 
541  // more than one pixel away
542  if (fabs (off_x) > 1.0 || fabs (off_y) > 1.0)
543  return false;
544  else {
545  //cout<<"Succesfully refined point. old was "<<x<<" "<<y<<"with q="
546  // <<NewCornerness;
547  x = col+off_x;
548  y = row+off_y;
549  // now compute a because we need it for the cornerness value:
550  // 1/9 * [ -1 +2 -1 ]
551  // [ +2 +5 +2 ]
552  // [ -1 +2 -1 ]
553  float a = (- DR_P11 +2.0f *DR_P12 - DR_P13
554  +2.0f*DR_P21 +5.0f *DR_P22 +2.0f*DR_P23
555  - DR_P31 +2.0f *DR_P32 - DR_P33)/9.0f;
556 
557  NewCornerness = float( a + b*off_x + c*off_y + 0.5*d*off_x*off_x
558  + e*off_x*off_y + 0.5*f*off_y*off_y );
559  //cout<<" New is "<<x<<" "<<y<<"with q=" <<NewCornerness <<endl;
560  return true;
561  }
562  }
563  // it's a saddle surface, but the corner may be good anyway :
564  else if (det < 0.0) {
565  x = col;
566  y = row;
567  return true;
568  }
569 
570  // det is zero !
571  return false;
572 }
573 
574 
575 template <class StorageType, class CalculationType>
576 int CornerDetectorGradient<StorageType, CalculationType>::
577 _ComputeCornerness(Image<CalculationType>& i) {
578  BIASERR("Make instances of derived class, dont call this !!!");
579  BIASABORT;
580  return 0;
581 }
582 
583 template <class StorageType, class CalculationType>
584 void CornerDetectorGradient<StorageType, CalculationType>::
585 _FillFeatureMap(vector<HomgPoint2D>& p)
586 {
587  vector<HomgPoint2D>::iterator pit;
588  const int md=_MinDistance-1;
589  const int width=_Cornerness.GetWidth();
590  const int height=_Cornerness.GetHeight();
591  const int mwidth=width-md;
592  const int mheight=height-md;
593  int minx, maxx, miny, maxy, mx, my, x, y;
594  for (pit=p.begin(); pit!=p.end(); pit++){
595  if (!pit->IsAtInfinity()){
596  _NumberFeaturesSoFar++;
597  x=(int)rint((*pit)[0]);
598  y=(int)rint((*pit)[1]);
599  // fill feature map
600  miny=(y>md)?(y-md):(0);
601  maxy=(y>=mheight)?(height-1):(y+md);
602  minx=(x>md)?(x-md):(0);
603  maxx=(x>=mwidth)?(width-1):(x+md);
604  for (my=miny; my<=maxy; my++){
605  for (mx=minx; mx<=maxx; mx++){
606  _FeatureMap[my][mx]=true;
607  }
608  }
609  }
610  }
611 }
612 
613 template <class StorageType, class CalculationType>
614 void CornerDetectorGradient<StorageType, CalculationType>::
615 _DeleteInternalMem()
616 {
617  if (!_Cornerness.IsEmpty()) _Cornerness.Release();
618 
619  if (_FeatureMap!=NULL) {
620  delete[] _FeatureMap[0]; _FeatureMap[0]=NULL;
621  delete[] _FeatureMap; _FeatureMap=NULL;
622  }
623 }
624 
625 template <class StorageType, class CalculationType>
626 void CornerDetectorGradient<StorageType, CalculationType>::
627 _AllocInternalMem(const int width, const int height)
628 {
629  _DeleteInternalMem();
630  int size=width*height;
631  _Cornerness.Init(width, height, 1);
632 // _sgxx.Init(width, height, 1);
633 // _sgxy.Init(width, height, 1);
634 // _sgyy.Init(width, height, 1);
635  _FeatureMap = new bool*[height];
636  _FeatureMap[0] = new bool[size];
637  for (int i=1; i<height; i++) {
638  _FeatureMap[i]=_FeatureMap[i-1]+width;
639  }
640  _FeatList.reserve(size);
641 }
642 
643 template <class StorageType, class CalculationType>
644 int CornerDetectorGradient<StorageType, CalculationType>::
645 _CalcCornerness(const Image<StorageType>& image,
646  Image<CalculationType>& cornerness)
647 {
648  Image<StorageType> fim;
649  _gauss.Filter(image, fim);
650 #ifdef BIAS_DEBUG
651  if (this->DebugLevelIsSet(D_CD_WRITE_DEBUG_IM)){
652  Image<StorageType> im2=image;
653  //ImageIO::Save("float-input.mip", im2);
654  //ImageIO::Save("gauss-filtered.mip", fim);
655  ImageIO::Save("float-input.mip", im2);
656  ImageIO::Save("gauss-filtered.mip", fim);
657  }
658 #endif
659 
660 
661  Image<CalculationType> *kim = NULL;
662  if (typeid(StorageType)==typeid(CalculationType)) {
663  kim = (Image<CalculationType> *) &fim;
664  } else {
665  kim = new Image<CalculationType>(fim.GetWidth(), fim.GetHeight(),
666  fim.GetChannelCount());
667  ImageConvert::ConvertST(fim, *kim, kim->GetStorageType());
668  }
669  int sres=this->_st.CalcStructureTensor(*kim, this->_sgxx, this->_sgxy,
670  this->_sgyy);
671  if (typeid(StorageType)!=typeid(CalculationType)) delete kim;
672 
673  if (this->_sgxx.GetWidth()>this->_Cornerness.GetWidth() ||
674  this->_sgxx.GetHeight()>this->_Cornerness.GetHeight()){
675  _AllocInternalMem(this->_sgxx.GetWidth(), this->_sgxx.GetHeight());
676  }
677 
678 
679  int mres = _ComputeCornerness(cornerness);
680  return (sres==0 && mres==0)?0:-1;
681 }
682 
683 
684 
685 template <class StorageType, class CalculationType>
686 int CornerDetectorGradient<StorageType, CalculationType>::
687 _CalcCornerness(const Image<CalculationType>& gradx,
688  const Image<CalculationType>& grady,
689  Image<CalculationType>& cornerness) {
690  int ires, meres;
691  if ((ires=this->_st.CalcStructureTensorValid(gradx, grady,
692  this->_sgxx, this->_sgxy,
693  this->_sgyy))!=0){
694  BIASERR("error calculatimng structure tensor");
695  }
696  if (!cornerness.SamePixelAndChannelCount(gradx)){
697  if (!cornerness.IsEmpty()) cornerness.Release();
698  cornerness.Init(gradx.GetWidth(), gradx.GetHeight(), 1);
699  }
700  meres = _ComputeCornerness(cornerness);
701 
702  return (ires==0 && meres==0)?0:-1;
703 }
704 
705 
706 template <class StorageType, class CalculationType>
707 void CornerDetectorGradient<StorageType, CalculationType>::
708 _ExtractLocalMaxima(const Image<CalculationType>& mineig)
709 {
710  const ROI* roi = mineig.GetROI();
711  int tlx = 0;
712  int tly = 0;
713  int brx = (int)roi->GetWidth() - 1;
714  int bry = (int)roi->GetHeight() - 1;
715  if (roi->GetROIType() == ROI_Corners)
716  mineig.GetROI()->GetCorners(tlx, tly, brx, bry);
717  const int minx = tlx;
718  const int miny = tly;
719  const int maxx = brx;
720  const int maxy = bry;
721 
722  register int x, y, lx, ly;
723  register double ev;
724  register const double me = (double)this->_MinCornerness;
725  Feat fp;
726  this->_FeatList.clear();
727  this->_FeatList.reserve((brx-tlx)*(bry-tly));
728  const CalculationType **ime = mineig.GetImageDataArray();
729 
730  // test image due to different ROI types
731  // @todo Use more efficient strategies for ROI_Points and ROI_Rows!
732  switch (roi->GetROIType()) {
733  case ROI_Corners:
734  for (y=miny; y<maxy; y++){
735  for (x=minx; x<maxx; x++){
736  ev = ime[y][x];
737  if (ev>me){
738  const int hws_local_max = 1;
739  bool is_local_maximum = true;
740  for (ly=y-hws_local_max; ly<=y+hws_local_max && is_local_maximum; ly++){
741  if (ly<miny || ly>=maxy) continue;
742  for (lx=x-hws_local_max; lx<=x+hws_local_max && is_local_maximum; lx++){
743  if (lx<minx || lx>=maxx) continue;
744  if (lx==x && ly==y) continue;
745  if (ev <= ime[ly][lx])
746  is_local_maximum = false;
747  }
748  }
749  if (is_local_maximum){
750  fp.x=x;
751  fp.y=y;
752  fp.val=(ev>(double)INT_MAX)?(INT_MAX):((int)ev);
753  this->_FeatList.push_back(fp);
754  }
755  }
756  }
757  }
758  break;
759  case ROI_Mask:
760  case ROI_Points:
761  case ROI_Rows:
762  for (y=miny; y<maxy; y++){
763  for (x=minx; x<maxx; x++){
764  ev = ime[y][x];
765  if (ev>me && roi->IsInROI(x,y)){
766  const int hws_local_max = 1;
767  bool is_local_maximum = true;
768  for (ly=y-hws_local_max; ly<=y+hws_local_max && is_local_maximum; ly++){
769  if (ly<miny || ly>=maxy) continue;
770  for (lx=x-hws_local_max; lx<=x+hws_local_max && is_local_maximum; lx++){
771  if (lx<minx || lx>=maxx) continue;
772  if (lx==x && ly==y) continue;
773  if (!roi->IsInROI(lx,ly)) continue;
774  if (ev <= ime[ly][lx])
775  is_local_maximum = false;
776  }
777  }
778  if (is_local_maximum){
779  fp.x=x;
780  fp.y=y;
781  fp.val=(ev>(double)INT_MAX)?(INT_MAX):((int)ev);
782  this->_FeatList.push_back(fp);
783  }
784  }
785  }
786  }
787  break;
788  default:
789  BIASERR("Unknown ROI type " << roi->GetROIType() << " found!");
790  }
791 }
792 
793 
794 //////////////////////////////////////////////////////////////////////////
795 // instantiation
796 //////////////////////////////////////////////////////////////////////////
797 
798 namespace BIAS{
799 template class CornerDetectorGradient<unsigned char, float>;
800 template class CornerDetectorGradient<float, float>;
801 
802 // fill in instances as required
803 #ifdef BUILD_IMAGE_INT
804 template class CornerDetectorGradient<unsigned char, int>;
805 #endif
806 
807 #ifdef BUILD_IMAGE_CHAR
808 template class CornerDetectorGradient<unsigned char, char>;
809 #endif
810 
811 #ifdef BUILD_IMAGE_SHORT
812 #endif
813 
814 #ifdef BUILD_IMAGE_USHORT
815 #ifdef BUILD_IMAGE_INT
816 template class CornerDetectorGradient<unsigned short, int>;
817 #endif
818 template class CornerDetectorGradient<unsigned short, float>;
819 #endif
820 
821 #ifdef BUILD_IMAGE_UINT
822 #endif
823 
824 #ifdef BUILD_IMAGE_DOUBLE
825 #endif
826 }
BIAS::Image::Release
void Release()
reimplemented from ImageBase
Definition: Image.cpp:1579
BIAS::CornerDetectorGradient::_RefineCornerPosition
bool _RefineCornerPosition(int col, int row, double &x, double &y, QUAL &SubPixelCornerness)
fit quadratic surface around (row;col) to find the maximum cornerness subpixel position and save it i...
Definition: CornerDetectorGradient.cpp:469
BIAS::CornerDetectorGradient::_MinDistance
int _MinDistance
minimal distance between points
Definition: CornerDetectorGradient.hh:309
BIAS::Matrix2x2::Invert
int Invert(Matrix2x2< T > &result) const
analyticaly inverts matrix
Definition: Matrix2x2.cpp:115
BIAS::ROI
class for handling different region of interest (ROI) representations...
Definition: ROI.hh:118
BIAS::ROI::GetROIType
enum ERoiType GetROIType() const
is the mask image valid?
Definition: ROI.hh:303
BIAS::HomgPoint2D
class HomgPoint2D describes a point with 2 degrees of freedom in projective coordinates.
Definition: HomgPoint2D.hh:67
BIAS::ROI::GetCorners
void GetCorners(unsigned &UpperLeftX, unsigned &UpperLeftY, unsigned &LowerRightX, unsigned &LowerRightY) const
Return the region of interest, by saving the coordinates within the variables defined by the paramete...
Definition: ROI.hh:443
BIAS::CornerDetectorGradient::_DoRefinement
bool _DoRefinement
tells whether subpixel refinement should be done
Definition: CornerDetectorGradient.hh:311
BIAS::ImageBase::IsEmpty
bool IsEmpty() const
check if ImageData_ points to allocated image buffer or not
Definition: ImageBase.hh:245
BIAS::CornerDetectorGradient::_ComputeCornerness
virtual int _ComputeCornerness(Image< CalculationType > &Cornerness)
the main function which should be overwritten in derived classes, no implementation in this class ...
Definition: CornerDetectorGradient.cpp:577
BIAS::CornerDetectorGradient::_ExtractLocalMaxima
void _ExtractLocalMaxima(const Image< CalculationType > &Cornerness)
extracts the local maxima checking a 3x3 region into _FeatList
Definition: CornerDetectorGradient.cpp:708
BIAS::CornerDetectorGradient::DetectFromCornerness
virtual int DetectFromCornerness(const Image< CalculationType > &cornerness, std::vector< HomgPoint2D > &p, std::vector< QUAL > &quality)
detect corners avoiding double calculation of cornerness image
Definition: CornerDetectorGradient.cpp:301
BIAS::ROI::GetWidth
unsigned GetWidth() const
width capacity of roi (image width)
Definition: ROI.hh:250
BIAS::FilterBase
virtual parent class for API definition of all (future) filters
Definition: FilterBase.hh:77
BIAS::CornerDetectorGradient::_MinCornerness
double _MinCornerness
threshold to accept a feature
Definition: CornerDetectorGradient.hh:310
BIAS::ImageBase::GetWidth
unsigned int GetWidth() const
Definition: ImageBase.hh:312
BIAS::CornerDetectorGradient::GetCov_
void GetCov_(const std::vector< HomgPoint2D > &p, std::vector< Matrix2x2< double > > &cov, const double noise=5) const
returns the inverse structure tensor, multiplied by noise^2
Definition: CornerDetectorGradient.cpp:137
BIAS::CornerDetectorGradient::_FeatureMap
bool ** _FeatureMap
array used to enforce a min distance / neighbor suppression
Definition: CornerDetectorGradient.hh:296
BIAS::CornerDetectorGradient::_AllocInternalMem
void _AllocInternalMem(const int width, const int height)
allocates _Cornerness , _sgxx, _sgxy and _sgyy
Definition: CornerDetectorGradient.cpp:627
BIAS::ImageBase::GetROI
ROI * GetROI()
Returns a pointer to the roi object.
Definition: ImageBase.hh:615
BIAS::ImageConvert::ConvertST
static int ConvertST(const BIAS::ImageBase &source, BIAS::ImageBase &dest, ImageBase::EStorageType targetST)
Function to convert the storage type of images e.g.
Definition: ImageConvert.cpp:2841
BIAS::ROI::IsInROI
bool IsInROI(const double &x, const double &y) const
ROI check if pixel position is inside the ROI.
Definition: ROI.hh:463
BIAS::ImageBase::GetChannelCount
unsigned int GetChannelCount() const
returns the number of Color channels, e.g.
Definition: ImageBase.hh:382
BIAS::ImageBase::GetHeight
unsigned int GetHeight() const
Definition: ImageBase.hh:319
BIAS::Image
The image template class for specific storage types.
Definition: Image.hh:78
BIAS::CornerDetectorGradient::_FillFeatureMap
void _FillFeatureMap(std::vector< HomgPoint2D > &p)
fills feature map with every point from p which isnt at infinity
Definition: CornerDetectorGradient.cpp:585
BIAS::ImageBase::SamePixelAndChannelCount
bool SamePixelAndChannelCount(const ImageBase &Image) const
checks if data area has same &quot;size&quot; as Image of other type
Definition: ImageBase.hh:73
BIAS::ImageIO::Save
static int Save(const std::string &filename, const ImageBase &img, const enum TFileFormat FileFormat=FF_auto, const bool sync=BIAS_DEFAULT_SYNC, const int c_jpeg_quality=BIAS_DEFAULT_IMAGE_QUALITY, const bool forceNewID=BIAS_DEFAULT_FORCENEWID, const bool &writeMetaData=true)
Export image as file using extrnal libs.
Definition: ImageIO.cpp:725
BIAS::TimeMeasure::GetRealTime
double GetRealTime() const
return real time (=wall time clock) in usec JW For Win32: real-time is measured differently from user...
Definition: TimeMeasure.cpp:370
BIAS::Image::Init
void Init(unsigned int Width, unsigned int Height, unsigned int channels=1, enum EStorageType storageType=ST_unsignedchar, const bool interleaved=true)
calls Init from ImageBase storageType is ignored, just dummy argument
Definition: Image.cpp:421
BIAS::CornerDetectorBase
purly virtual interface defining class for corner detectors
Definition: CornerDetectorBase.hh:53
BIAS::ImageBase::GetStorageType
enum EStorageType GetStorageType() const
Definition: ImageBase.hh:414
BIAS::CornerDetectorGradient
base class for all gradient based corner detectors
Definition: CornerDetectorGradient.hh:87
BIAS::ROI::GetHeight
unsigned GetHeight() const
height capacity of roi (image height)
Definition: ROI.hh:254
BIAS::CornerDetectorGradient::_EnforceMinimumDistance
int _EnforceMinimumDistance(std::vector< HomgPoint2D > &p, std::vector< QUAL > &quality)
selects features sorted by quality from _FeatList such that they have a min distance as given in the ...
Definition: CornerDetectorGradient.cpp:373
BIAS::HomgPoint2D::Set
void Set(const HOMGPOINT2D_TYPE &x, const HOMGPOINT2D_TYPE &y)
set elementwise with given 2 euclidean scalar values.
Definition: HomgPoint2D.hh:174
BIAS::CornerDetectorGradient::_CalcCornerness
int _CalcCornerness(const Image< StorageType > &src, Image< CalculationType > &cornerness)
calculates the cornerness using the src image
Definition: CornerDetectorGradient.cpp:645
BIAS::TimeMeasure
class TimeMeasure contains functions for timing real time and cpu time.
Definition: TimeMeasure.hh:111
BIAS::CornerDetectorGradient::Cornerness
virtual int Cornerness(const Image< CalculationType > &gradx, const Image< CalculationType > &grady, Image< CalculationType > &cornerness)
calculates the cornerness using the gradients gx and gy
Definition: CornerDetectorGradient.cpp:274
BIAS::CornerDetectorGradient::_DeleteInternalMem
void _DeleteInternalMem()
frees _Cornerness, _sgxx, _sgxy and _sgyy
Definition: CornerDetectorGradient.cpp:615
BIAS::CornerDetectorGradient::Detect
virtual int Detect(const Image< StorageType > &image, std::vector< HomgPoint2D > &p, std::vector< QUAL > &quality, std::vector< Matrix2x2< double > > *cov=NULL)
detect corners in a grey image
Definition: CornerDetectorGradient.cpp:71
BIAS::Image::GetImageDataArray
const StorageType ** GetImageDataArray() const
overloaded GetImageDataArray() from ImageBase
Definition: Image.hh:153
Generated on Tue Feb 6 2018 02:01:34 for Basic Image AlgorithmS Library by   doxygen 1.8.5