Basic Image AlgorithmS Library  2.8.0
 All Classes Namespaces Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StructureTensor.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
6  Institut fuer Informatik
7  Christian-Albrechts-Universitaet Kiel
8 
9 
10 BIAS is free software; you can redistribute it and/or modify
11 it under the terms of the GNU Lesser General Public License as published by
12 the Free Software Foundation; either version 2.1 of the License, or
13 (at your option) any later version.
14 
15 BIAS is distributed in the hope that it will be useful,
16 but WITHOUT ANY WARRANTY; without even the implied warranty of
17 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 GNU Lesser General Public License for more details.
19 
20 You should have received a copy of the GNU Lesser General Public License
21 along with BIAS; if not, write to the Free Software
22 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
23 */
24 
25 #include "StructureTensor.hh"
26 
27 #ifdef BIAS_DEBUG
28 #include <Base/Image/ImageIO.hh>
29 #endif
30 
31 #include <Base/Common/BIASpragma.hh>
32 
33 
34  //#define TIME_MEASURE
35 
36 #ifdef TIME_MEASURE
37 #include <Base/Debug/TimeMeasure.hh>
38 #endif
39 
40 using namespace BIAS;
41 using namespace std;
42 
43 //////////////////////////////////////////////////////////////////////////
44 // implementation
45 //////////////////////////////////////////////////////////////////////////
46 
47 template <class InputStorageType, class OutputStorageType>
50  FilterNTo3N<InputStorageType, OutputStorageType>()
51 {
53  _HalfWinSize=3;
55 }
56 
57 template <class InputStorageType, class OutputStorageType>
60  : _HalfWinSize(other._HalfWinSize),
61  _InternalMemWidth(other._InternalMemWidth),
62  _InternalMemHeight(other._InternalMemHeight)
63 {
64  if (other._grad)
65  _grad = other._grad->Clone();
66  else
67  _grad = NULL;
68 }
69 
70 template <class InputStorageType, class OutputStorageType>
72 {
73  if (_grad) delete _grad;
74 }
75 
76 template <class InputStorageType, class OutputStorageType>
79 {
80  BIASCDOUT(D_FILTERBASE_CALLSTACK, "StructureTensor::Filter(src, dst)\n");
81  BIASERR("unfinished");
82  BIASABORT;
83  return -10;
84 }
85 
86 template <class InputStorageType, class OutputStorageType>
89 {
90  BIASCDOUT(D_FILTERBASE_CALLSTACK, "StructureTensor::FilterInt(src, dst)\n");
91  BIASERR("unfinished");
92  BIASABORT;
93  return -10;
94 }
95 
96 template <class InputStorageType, class OutputStorageType>
99 {
100  BIASCDOUT(D_FILTERBASE_CALLSTACK, "StructureTensor::FilterFloat(src, dst)\n");
101  BIASERR("unfinished");
102  BIASABORT;
103  return -10;
104 }
105 
106 template <class InputStorageType, class OutputStorageType>
112 {
113  BIASCDOUT(D_FILTERBASE_CALLSTACK, "StructureTensor::Filter(src, d1, d2, d3)\n");
114  BIASERR("unfinished");
115  BIASABORT;
116  return -10;
117 }
118 
119 template <class InputStorageType, class OutputStorageType>
125 {
126  BIASCDOUT(D_FILTERBASE_CALLSTACK, "StructureTensor::FilterInt(src, d1, d2, d3)\n");
127  BIASERR("unfinished");
128  BIASABORT;
129  return -10;
130 }
131 
132 template <class InputStorageType, class OutputStorageType>
138 {
139  BIASCDOUT(D_FILTERBASE_CALLSTACK, "StructureTensor::FilterFloat(src, d1, d2, d3)\n");
140 
141  BIASERR("unfinished");
142  BIASABORT;
143  return -10;
144 }
145 
146 
147 
148 
149 
150 
151 
152 
153 
154 
155 
156 
157 
158 
159 
160 
161 
162 
163 
164 
165 
166 
167 template <class InputStorageType, class OutputStorageType>
173 {
174  int gres=_grad->Filter(image, _gx, _gy);
175 #ifdef BIAS_DEBUG
176  if (DebugLevelIsSet(D_ST_WRITE_DEBUG_IM)){
177  //if (ImageIO::Save("gx.mip", _gx)!=0){
178  if (ImageIO::Save("gx.mip", _gx)!=0){
179  BIASERR("error writing gx");
180  }
181  //if (ImageIO::Save("gy.mip", _gy)!=0){
182  if (ImageIO::Save("gy.mip", _gy)!=0){
183  BIASERR("error writing gy");
184  }
185  }
186 #endif
187 
188  int sres = CalcStructureTensor(_gx, _gy, sgxx, sgxy, sgyy);
189 
190 #ifdef BIAS_DEBUG
191  if (DebugLevelIsSet(D_ST_WRITE_DEBUG_IM)){
192  _WriteDebugImage(sgxx, "sgxx.mip");
193  _WriteDebugImage(sgxy, "sgxy.mip");
194  _WriteDebugImage(sgyy, "sgyy.mip");
195  }
196 #endif
197  return (sres==0 && gres==0)?0:-1;
198 }
199 
200 template <class InputStorageType, class OutputStorageType>
203  const Image<OutputStorageType>& gy,
207 {
208  int sres;
209  // initialize the destination images
210  if (!gx.SamePixelAndChannelCount(sgxx)){
211  if (!sgxx.IsEmpty()) sgxx.Release();
212  sgxx.Init(gx.GetWidth(), gx.GetHeight(), 1);
213  }
214  if (!gx.SamePixelAndChannelCount(sgxy)){
215  if (!sgxy.IsEmpty()) sgxy.Release();
216  sgxy.Init(gx.GetWidth(), gx.GetHeight(), 1);
217  }
218  if (!gx.SamePixelAndChannelCount(sgyy)){
219  if (!sgyy.IsEmpty()) sgyy.Release();
220  sgyy.Init(gx.GetWidth(), gx.GetHeight(), 1);
221  }
222  // calculate it
223  // if (*_FilterBorderHandling==FilterBase<InputStorageType>::TBH_valid){
224  switch (_HalfWinSize){
225  case 1:
226  sres=CalcStructureTensor3x3(gx, gy, sgxx, sgxy, sgyy);
227  break;
228  case 2:
229  sres=CalcStructureTensor5x5(gx, gy, sgxx, sgxy, sgyy);
230  break;
231  case 3:
232  sres=CalcStructureTensor7x7(gx, gy, sgxx, sgxy, sgyy);
233  break;
234  default:
235  sres=CalcStructureTensorValid(gx, gy, sgxx, sgxy, sgyy);
236  break;
237  }
238 
239 
240  return sres;
241 }
242 
243 template <class InputStorageType, class OutputStorageType>
246  const Image<OutputStorageType>& gy,
250 {
251 
252 #ifdef BIAS_DEBUG
253  BIASASSERT(gx.SamePixelAndChannelCount(gy));
254  BIASASSERT(gx.SamePixelAndChannelCount(sgxx));
255  BIASASSERT(gx.SamePixelAndChannelCount(sgxy));
256  BIASASSERT(gx.SamePixelAndChannelCount(sgyy));
257 #endif
258 
259 
260  const int width=gx.GetWidth();
261  const int height=gx.GetHeight();
262  if (width>_InternalMemWidth || height>_InternalMemHeight){
263  //_DeleteInternalMem();
264  _AllocInternalMem(width, height);
265  }
266 
267  _GradientProducts(gx, gy);
268 
269  // propagate ROI
270  int tlx, tly, brx, bry;
271  _gxx.GetROI()->GetCorners(tlx, tly, brx, bry);
272  const int hws=_HalfWinSize;
273  const int minx=tlx+hws;
274  const int maxx=brx-hws;
275  const int miny=tly+hws;
276  const int maxy=bry-hws;
277  sgxx.GetROI()->SetCorners(minx, miny, maxx, maxy);
278  sgxy.GetROI()->SetCorners(minx, miny, maxx, maxy);
279  sgyy.GetROI()->SetCorners(minx, miny, maxx, maxy);
280 
281  register int x, y, ix, iy, mx, my;
282  register OutputStorageType **isgxx=sgxx.GetImageDataArray();
283  register OutputStorageType **isgxy=sgxy.GetImageDataArray();
284  register OutputStorageType **isgyy=sgyy.GetImageDataArray();
285  register OutputStorageType **igxx=_gxx.GetImageDataArray();
286  register OutputStorageType **igxy=_gxy.GetImageDataArray();
287  register OutputStorageType **igyy=_gyy.GetImageDataArray();
288 
289  for (y=miny; y<maxy; y++){
290  for (x=minx; x<maxx; x++){
291  isgxx[y][x]=isgxy[y][x]=isgyy[y][x]=(OutputStorageType)0.0;
292  for (iy=-hws; iy<=hws; iy++){
293  my=y+iy;
294  for (ix=-hws; ix<=hws; ix++){
295  mx=x+ix;
296  isgxx[y][x]+=igxx[my][mx];
297  isgxy[y][x]+=igxy[my][mx];
298  isgyy[y][x]+=igyy[my][mx];
299  }
300  }
301 #ifdef BIAS_DEBUG
302  if (DebugLevelIsSet(D_ST_FEATURES) && y<miny+2 && x<minx+2){
303  cerr << "( "<<x<<", "<<y<<") sgxx "<<isgxx[y][x]<<" sgxy "
304  <<isgxy[y][x]<<" sgyy "<<isgyy[y][x]<<endl;
305  }
306 #endif
307  }
308  }
309  return 0;
310 }
311 
312 
313 template <class InputStorageType, class OutputStorageType>
316  const Image<OutputStorageType>& gy,
320 {
321 
322 #ifdef BIAS_DEBUG
323  BIASASSERT(gx.SamePixelAndChannelCount(gy));
324  BIASASSERT(gx.SamePixelAndChannelCount(sgxx));
325  BIASASSERT(gx.SamePixelAndChannelCount(sgxy));
326  BIASASSERT(gx.SamePixelAndChannelCount(sgyy));
327 #endif
328 
329  const int width=gx.GetWidth();
330  const int height=gx.GetHeight();
331  if (width>_InternalMemWidth || height>_InternalMemHeight){
332  //_DeleteInternalMem();
333  _AllocInternalMem(width, height);
334  }
335 
336  _GradientProducts(gx, gy);
337 
338  // propagate ROI
339  int tlx, tly, brx, bry;
340  _gxx.GetROI()->GetCorners(tlx, tly, brx, bry);
341  const int minx=tlx+1;
342  const int maxx=brx-1;
343  const int miny=tly+1;
344  const int maxy=bry-1;
345  sgxx.GetROI()->SetCorners(minx, miny, maxx, maxy);
346  sgxy.GetROI()->SetCorners(minx, miny, maxx, maxy);
347  sgyy.GetROI()->SetCorners(minx, miny, maxx, maxy);
348 
349  register int x, y, xm, xp, ym, yp;
350  register OutputStorageType **isgxx=sgxx.GetImageDataArray();
351  register OutputStorageType **isgxy=sgxy.GetImageDataArray();
352  register OutputStorageType **isgyy=sgyy.GetImageDataArray();
353  register OutputStorageType **igxx=_gxx.GetImageDataArray();
354  register OutputStorageType **igxy=_gxy.GetImageDataArray();
355  register OutputStorageType **igyy=_gyy.GetImageDataArray();
356  for (y=miny; y<maxy; y++){
357  for (x=minx; x<maxx; x++){
358  xm=x-1; xp=x+1;
359  ym=y-1; yp=y+1;
360  isgxx[y][x]=
361  igxx[ym][xm]+igxx[ym][x]+igxx[ym][xp]+
362  igxx[y][xm]+igxx[y][x]+igxx[y][xp]+
363  igxx[yp][xm]+igxx[yp][x]+igxx[yp][xp];
364  isgxy[y][x]=
365  igxy[ym][xm]+igxy[ym][x]+igxy[ym][xp]+
366  igxy[y][xm]+igxy[y][x]+igxy[y][xp]+
367  igxy[yp][xm]+igxy[yp][x]+igxy[yp][xp];
368  isgyy[y][x]=
369  igyy[ym][xm]+igyy[ym][x]+igyy[ym][xp]+
370  igyy[y][xm]+igyy[y][x]+igyy[y][xp]+
371  igyy[yp][xm]+igyy[yp][x]+igyy[yp][xp];
372  }
373  }
374 #ifdef BIAS_DEBUG
375  if (DebugLevelIsSet(D_ST_WRITE_DEBUG_IM)){
376  _WriteDebugImage(sgxx, "sgxx-3x3.mip");
377  _WriteDebugImage(sgxy, "sgxy-3x3.mip");
378  _WriteDebugImage(sgyy, "sgyy-3x3.mip");
379 
380  CalcStructureTensorValid(gx, gy, sgxx, sgxy, sgyy);
381  _WriteDebugImage(sgxx, "sgxx-st.mip");
382  _WriteDebugImage(sgxy, "sgxy-st.mip");
383  _WriteDebugImage(sgyy, "sgyy-st.mip");
384  }
385 #endif
386 
387  return 0;
388 }
389 
390 template <class InputStorageType, class OutputStorageType>
393  const Image<OutputStorageType>& gy,
397 {
398 
399 #ifdef TIME_MEASURE
400  TimeMeasure timer;
401  timer.Start();
402 #endif
403 
404 #ifdef BIAS_DEBUG
405  BIASASSERT(gx.SamePixelAndChannelCount(gy));
406  BIASASSERT(gx.SamePixelAndChannelCount(sgxx));
407  BIASASSERT(gx.SamePixelAndChannelCount(sgxy));
408  BIASASSERT(gx.SamePixelAndChannelCount(sgyy));
409 #endif
410 
411  const int width=gx.GetWidth();
412  const int height=gx.GetHeight();
413  if (width>_InternalMemWidth || height>_InternalMemHeight){
414  //_DeleteInternalMem();
415  _AllocInternalMem(width, height);
416  }
417 
418  _GradientProducts(gx, gy);
419 
420 #ifdef TIME_MEASURE
421  timer.Stop();
422  cerr << "_GradientProducts took "<<timer.GetRealTime()<<" us\n";
423  timer.Reset();
424  timer.Start();
425 #endif
426 
427  // propagate ROI
428  int tlx, tly, brx, bry;
429  _gxx.GetROI()->GetCorners(tlx, tly, brx, bry);
430  const int miny=tly;
431  const int maxy=bry;
432  const int minxp=tlx+2;
433  const int maxxm=brx-2;
434  const int minyp=tly+2;
435  const int maxym=bry-2;
436  sgxx.GetROI()->SetCorners(minxp, minyp, maxxm, maxym);
437  sgxy.GetROI()->SetCorners(minxp, minyp, maxxm, maxym);
438  sgyy.GetROI()->SetCorners(minxp, minyp, maxxm, maxym);
439 
440  register int x, y, xm, xmm, xp, xpp, ym, ymm, yp, ypp;
441  register OutputStorageType **isgxx=sgxx.GetImageDataArray();
442  register OutputStorageType **isgxy=sgxy.GetImageDataArray();
443  register OutputStorageType **isgyy=sgyy.GetImageDataArray();
444  register OutputStorageType **itsgxx=_tsgxx.GetImageDataArray();
445  register OutputStorageType **itsgxy=_tsgxy.GetImageDataArray();
446  register OutputStorageType **itsgyy=_tsgyy.GetImageDataArray();
447  register OutputStorageType **igxx=_gxx.GetImageDataArray();
448  register OutputStorageType **igxy=_gxy.GetImageDataArray();
449  register OutputStorageType **igyy=_gyy.GetImageDataArray();
450  // separate loops are faster
451  for (y=miny; y<maxy; y++){
452  for (x=minxp; x<maxxm; x++){
453  xmm=x-2; xm=x-1; xp=x+1; xpp=x+2;
454  itsgxx[y][x]=
455  igxx[y][xmm]+igxx[y][xm]+igxx[y][x]+igxx[y][xp]+igxx[y][xpp];
456  }
457  }
458  for (y=minyp; y<maxym; y++){
459  for (x=minxp; x<maxxm; x++){
460  ymm=y-2; ym=y-1; yp=y+1; ypp=y+2;
461  isgxx[y][x]=
462  itsgxx[ymm][x]+itsgxx[ym][x]+itsgxx[y][x]+itsgxx[yp][x]+itsgxx[ypp][x];
463  }
464  }
465  for (y=miny; y<maxy; y++){
466  for (x=minxp; x<maxxm; x++){
467  xmm=x-2; xm=x-1; xp=x+1; xpp=x+2;
468  itsgxy[y][x]=
469  igxy[y][xmm]+igxy[y][xm]+igxy[y][x]+igxy[y][xp]+igxy[y][xpp];
470  }
471  }
472  for (y=minyp; y<maxym; y++){
473  for (x=minxp; x<maxxm; x++){
474  ymm=y-2; ym=y-1; yp=y+1; ypp=y+2;
475  isgxy[y][x]=
476  itsgxy[ymm][x]+itsgxy[ym][x]+itsgxy[y][x]+itsgxy[yp][x]+itsgxy[ypp][x];
477  }
478  }
479  for (y=miny; y<maxy; y++){
480  for (x=minxp; x<maxxm; x++){
481  xmm=x-2; xm=x-1; xp=x+1; xpp=x+2;
482  itsgyy[y][x]=
483  igyy[y][xmm]+igyy[y][xm]+igyy[y][x]+igyy[y][xp]+igyy[y][xpp];
484  }
485  }
486  for (y=minyp; y<maxym; y++){
487  for (x=minxp; x<maxxm; x++){
488  ymm=y-2; ym=y-1; yp=y+1; ypp=y+2;
489  isgyy[y][x]=
490  itsgyy[ymm][x]+itsgyy[ym][x]+itsgyy[y][x]+itsgyy[yp][x]+itsgyy[ypp][x];
491  }
492  }
493 
494 #ifdef BIAS_DEBUG
495  if (DebugLevelIsSet(D_ST_WRITE_DEBUG_IM)){
496  _WriteDebugImage(sgxx, "sgxx-5x5.mip");
497  _WriteDebugImage(sgxy, "sgxy-5x5.mip");
498  _WriteDebugImage(sgyy, "sgyy-5x5.mip");
499 
500  CalcStructureTensorValid(gx, gy, sgxx, sgxy, sgyy);
501  _WriteDebugImage(sgxx, "sgxx-st.mip");
502  _WriteDebugImage(sgxy, "sgxy-st.mip");
503  _WriteDebugImage(sgyy, "sgyy-st.mip");
504  }
505 #endif
506 
507 #ifdef TIME_MEASURE
508  timer.Stop();
509  cerr << "CalcStructureTensor5x5 without _GradientProducts took "<<timer.GetRealTime()<<" us\n";
510 #endif
511  return 0;
512 }
513 
514 template <class InputStorageType, class OutputStorageType>
517  const Image<OutputStorageType>& gy,
521 {
522 #ifdef TIME_MEASURE
523  TimeMeasure timer;
524  timer.Start();
525 #endif
526 
527 #ifdef BIAS_DEBUG
528  BIASASSERT(gx.SamePixelAndChannelCount(gy));
529  BIASASSERT(gx.SamePixelAndChannelCount(sgxx));
530  BIASASSERT(gx.SamePixelAndChannelCount(sgxy));
531  BIASASSERT(gx.SamePixelAndChannelCount(sgyy));
532 #endif
533 
534  register int width=gx.GetWidth();
535  const int height=gx.GetHeight();
536  if (width>_InternalMemWidth || height>_InternalMemHeight){
537  _AllocInternalMem(width, height);
538  }
539 
540  _GradientProducts(gx, gy);
541 
542 #ifdef TIME_MEASURE
543  timer.Stop();
544  cerr << "_GradientProducts took "<<timer.GetRealTime()<<" us\n";
545  timer.Reset();
546  timer.Start();
547 #endif
548 
549  // propagate ROI
550  int tlx, tly, brx, bry;
551  _gxx.GetROI()->GetCorners(tlx, tly, brx, bry);
552  const int miny=tly;
553  const int maxy=bry;
554  const int minxp=tlx+3;
555  const int maxxm=brx-3;
556  const int minyp=tly+3;
557  const int maxym=bry-3;
558  sgxx.GetROI()->SetCorners(minxp, minyp, maxxm, maxym);
559  sgxy.GetROI()->SetCorners(minxp, minyp, maxxm, maxym);
560  sgyy.GetROI()->SetCorners(minxp, minyp, maxxm, maxym);
561 
562  {
563  register OutputStorageType *ptsgxx, *ptmsgxx, *pmgxx, *ppgxx, *lend, *end;
564  register OutputStorageType *ptsgxy, *ptmsgxy, *pmgxy, *ppgxy;
565  register OutputStorageType *ptsgyy, *ptmsgyy, *pmgyy, *ppgyy;
566  register const int tstep=width-(maxxm-minxp);
567  const int offset=(width*miny)+minxp;
568  ptmsgxx = _tsgxx.GetImageData() + offset;
569  ptmsgxy = _tsgxy.GetImageData() + offset;
570  ptmsgyy = _tsgyy.GetImageData() + offset;
571  pmgxx = _gxx.GetImageData() + offset-3;
572  pmgxy = _gxy.GetImageData() + offset-3;
573  pmgyy = _gyy.GetImageData() + offset-3;
574  lend = _tsgxx.GetImageData() + (width*miny)+maxxm;
575  end = _tsgxx.GetImageData() + ((maxy-1)*width)+maxxm;
576  while (ptmsgxx<end){
577  ppgxx = pmgxx;
578  ptsgxx = ptmsgxx;
579  *ptsgxx=(OutputStorageType)0.0;
580  *ptsgxx+=*ppgxx++;
581  *ptsgxx+=*ppgxx++;
582  *ptsgxx+=*ppgxx++;
583  *ptsgxx+=*ppgxx++;
584  *ptsgxx+=*ppgxx++;
585  *ptsgxx+=*ppgxx++;
586  *ptsgxx+=*ppgxx++;
587  ptsgxx++;
588 
589  ppgxy = pmgxy;
590  ptsgxy = ptmsgxy;
591  *ptsgxy=(OutputStorageType)0.0;
592  *ptsgxy+=*ppgxy++;
593  *ptsgxy+=*ppgxy++;
594  *ptsgxy+=*ppgxy++;
595  *ptsgxy+=*ppgxy++;
596  *ptsgxy+=*ppgxy++;
597  *ptsgxy+=*ppgxy++;
598  *ptsgxy+=*ppgxy++;
599  ptsgxy++;
600 
601  ppgyy = pmgyy;
602  ptsgyy = ptmsgyy;
603  *ptsgyy=(OutputStorageType)0.0;
604  *ptsgyy+=*ppgyy++;
605  *ptsgyy+=*ppgyy++;
606  *ptsgyy+=*ppgyy++;
607  *ptsgyy+=*ppgyy++;
608  *ptsgyy+=*ppgyy++;
609  *ptsgyy+=*ppgyy++;
610  *ptsgyy+=*ppgyy++;
611  ptsgyy++;
612  while (ptsgxx<lend){
613  *ptsgxx++ = *ptmsgxx++ - *pmgxx++ + *ppgxx++;
614  *ptsgxy++ = *ptmsgxy++ - *pmgxy++ + *ppgxy++;
615  *ptsgyy++ = *ptmsgyy++ - *pmgyy++ + *ppgyy++;
616  }
617  lend+=width;
618  pmgxx+=tstep;
619  ptmsgxx+=tstep;
620  pmgxy+=tstep;
621  ptmsgxy+=tstep;
622  pmgyy+=tstep;
623  ptmsgyy+=tstep;
624  }
625  }
626 
627 #ifdef TIME_MEASURE
628  timer.Stop();
629  cerr << "first convolution took "<<timer.GetRealTime()<<" us\n";
630  timer.Reset();
631  timer.Start();
632 #endif
633  /*
634  //const int max_cache_size=128000;
635  //const int dy=max_cache_size/width;
636  //int mminy=miny, mminyp=minyp, mmaxy, mmaxym;
637  //for (mmaxy=dy, mmaxym=dy-3; mmaxy<maxy;
638  // mmaxy+=dy, mmaxym+=dy, mminy+=dy, mmindyp+=dy)
639  { // this is very slow, maybe because of caches misses if done on the hole
640  // image
641  register OutputStorageType *ppsgxx, *pmsgxx, *ptmsgxx, *ptpsgxx, *lend, *end;
642  register OutputStorageType *ppsgxy, *pmsgxy, *ptmsgxy, *ptpsgxy;
643  register OutputStorageType *ppsgyy, *pmsgyy, *ptmsgyy, *ptpsgyy;
644  register const int tstep=(width*(maxym-minyp-1))-1;
645  const int offset=(width*miny)+minxp;
646 
647  ptmsgxx = _tsgxx.GetImageData() + offset;
648  ppsgxx = pmsgxx = sgxx.GetImageData() + (width*minyp)+minxp;
649  ptmsgxy = _tsgxy.GetImageData() + offset;
650  ppsgxy = pmsgxy = sgxy.GetImageData() + (width*minyp)+minxp;
651  ptmsgyy = _tsgyy.GetImageData() + offset;
652  ppsgyy = pmsgyy = sgyy.GetImageData() + (width*minyp)+minxp;
653 
654  lend = sgxx.GetImageData() + (width*(maxym))+minxp;
655  end = sgxx.GetImageData() + (width*(maxym))+maxxm-1;
656 
657  while (ppsgxx<end){
658  ppsgxx = pmsgxx + width;
659  ptpsgxx = ptmsgxx;
660  *pmsgxx=(OutputStorageType)0.0;
661  *pmsgxx+=*ptpsgxx; ptpsgxx+=width;
662  *pmsgxx+=*ptpsgxx; ptpsgxx+=width;
663  *pmsgxx+=*ptpsgxx; ptpsgxx+=width;
664  *pmsgxx+=*ptpsgxx; ptpsgxx+=width;
665  *pmsgxx+=*ptpsgxx; ptpsgxx+=width;
666  *pmsgxx+=*ptpsgxx; ptpsgxx+=width;
667  *pmsgxx+=*ptpsgxx; ptpsgxx+=width;
668 
669  ppsgxy = pmsgxy + width;
670  ptpsgxy = ptmsgxy;
671  *pmsgxy=(OutputStorageType)0.0;
672  *pmsgxy+=*ptpsgxy; ptpsgxy+=width;
673  *pmsgxy+=*ptpsgxy; ptpsgxy+=width;
674  *pmsgxy+=*ptpsgxy; ptpsgxy+=width;
675  *pmsgxy+=*ptpsgxy; ptpsgxy+=width;
676  *pmsgxy+=*ptpsgxy; ptpsgxy+=width;
677  *pmsgxy+=*ptpsgxy; ptpsgxy+=width;
678  *pmsgxy+=*ptpsgxy; ptpsgxy+=width;
679 
680  ppsgyy = pmsgyy + width;
681  ptpsgyy = ptmsgyy;
682  *pmsgyy=(OutputStorageType)0.0;
683  *pmsgyy+=*ptpsgyy; ptpsgyy+=width;
684  *pmsgyy+=*ptpsgyy; ptpsgyy+=width;
685  *pmsgyy+=*ptpsgyy; ptpsgyy+=width;
686  *pmsgyy+=*ptpsgyy; ptpsgyy+=width;
687  *pmsgyy+=*ptpsgyy; ptpsgyy+=width;
688  *pmsgyy+=*ptpsgyy; ptpsgyy+=width;
689  *pmsgyy+=*ptpsgyy; ptpsgyy+=width;
690 
691  while (ppsgxx<lend){
692  *ppsgxx = *pmsgxx - *ptmsgxx + *ptpsgxx;
693  ppsgxx+=width;
694  pmsgxx+=width;
695  ptmsgxx+=width;
696  ptpsgxx+=width;
697 
698  *ppsgxy = *pmsgxy - *ptmsgxy + *ptpsgxy;
699  ppsgxy+=width;
700  pmsgxy+=width;
701  ptmsgxy+=width;
702  ptpsgxy+=width;
703 
704  *ppsgyy = *pmsgyy - *ptmsgyy + *ptpsgyy;
705  ppsgyy+=width;
706  pmsgyy+=width;
707  ptmsgyy+=width;
708  ptpsgyy+=width;
709 
710  }
711  lend++;
712  ptmsgxx-=tstep;
713  pmsgxx-=tstep;
714  ptmsgxy-=tstep;
715  pmsgxy-=tstep;
716  ptmsgyy-=tstep;
717  pmsgyy-=tstep;
718  }
719  }
720 
721  */
722 
723 
724  // todo this really fast by pointer arithmetic
725  register int x, y, myp, mym, my1;
726  register const int ymmm=minyp-3, ymm=minyp-2, ym=minyp-1, yp=minyp+1,
727  ypp=minyp+2, yppp=minyp+3;
728  register OutputStorageType **isgxx=sgxx.GetImageDataArray();
729  register OutputStorageType **isgxy=sgxy.GetImageDataArray();
730  register OutputStorageType **isgyy=sgyy.GetImageDataArray();
731  register OutputStorageType **itsgxx=_tsgxx.GetImageDataArray();
732  register OutputStorageType **itsgxy=_tsgxy.GetImageDataArray();
733  register OutputStorageType **itsgyy=_tsgyy.GetImageDataArray();
734 
735  if (width*height>153600){ // the caches are running empty
736  // do concolution in two parts
737  const int medy=((maxym-minyp)>>1)+minyp;
738  const int eymmm=medy-3, eymm=medy-2, eym=medy-1, eyp=medy+1,
739  eypp=medy+2, eyppp=medy+3;
740  for (x=minxp; x<maxxm; x++){
741  isgxy[minyp][x]=itsgxy[ymmm][x]+itsgxy[ymm][x]+itsgxy[ym][x]+
742  itsgxy[minyp][x]+itsgxy[yp][x]+itsgxy[ypp][x]+itsgxy[yppp][x];
743  for (y=minyp+1; y<medy; y++){
744  myp=y+3; mym=y-4; my1=y-1;
745  isgxy[y][x]=isgxy[my1][x]-itsgxy[mym][x]+itsgxy[myp][x];
746  }
747  }
748  for (x=minxp; x<maxxm; x++){
749  isgxy[medy][x]=itsgxy[eymmm][x]+itsgxy[eymm][x]+itsgxy[eym][x]+
750  itsgxy[medy][x]+itsgxy[eyp][x]+itsgxy[eypp][x]+itsgxy[eyppp][x];
751  for (y=medy+1; y<maxym; y++){
752  myp=y+3; mym=y-4; my1=y-1;
753  isgxy[y][x]=isgxy[my1][x]-itsgxy[mym][x]+itsgxy[myp][x];
754  }
755  }
756 
757  for (x=minxp; x<maxxm; x++){
758  isgxx[minyp][x]=itsgxx[ymmm][x]+itsgxx[ymm][x]+itsgxx[ym][x]+
759  itsgxx[minyp][x]+itsgxx[yp][x]+itsgxx[ypp][x]+itsgxx[yppp][x];
760  for (y=minyp+1; y<medy; y++){
761  myp=y+3; mym=y-4; my1=y-1;
762  isgxx[y][x]=isgxx[my1][x]-itsgxx[mym][x]+itsgxx[myp][x];
763  }
764  }
765  for (x=minxp; x<maxxm; x++){
766  isgxx[medy][x]=itsgxx[eymmm][x]+itsgxx[eymm][x]+itsgxx[eym][x]+
767  itsgxx[medy][x]+itsgxx[eyp][x]+itsgxx[eypp][x]+itsgxx[eyppp][x];
768  for (y=medy+1; y<maxym; y++){
769  myp=y+3; mym=y-4; my1=y-1;
770  isgxx[y][x]=isgxx[my1][x]-itsgxx[mym][x]+itsgxx[myp][x];
771  }
772  }
773 
774  for (x=minxp; x<maxxm; x++){
775  isgyy[minyp][x]=itsgyy[ymmm][x]+itsgyy[ymm][x]+itsgyy[ym][x]+
776  itsgyy[minyp][x]+itsgyy[yp][x]+itsgyy[ypp][x]+itsgyy[yppp][x];
777  for (y=minyp+1; y<medy; y++){
778  myp=y+3; mym=y-4; my1=y-1;
779  isgyy[y][x]=isgyy[my1][x]-itsgyy[mym][x]+itsgyy[myp][x];
780  }
781  }
782  for (x=minxp; x<maxxm; x++){
783  isgyy[medy][x]=itsgyy[eymmm][x]+itsgyy[eymm][x]+itsgyy[eym][x]+
784  itsgyy[medy][x]+itsgyy[eyp][x]+itsgyy[eypp][x]+itsgyy[eyppp][x];
785  for (y=medy+1; y<maxym; y++){
786  myp=y+3; mym=y-4; my1=y-1;
787  isgyy[y][x]=isgyy[my1][x]-itsgyy[mym][x]+itsgyy[myp][x];
788  }
789  }
790  } else {
791  for (x=minxp; x<maxxm; x++){
792  isgxy[minyp][x]=itsgxy[ymmm][x]+itsgxy[ymm][x]+itsgxy[ym][x]+
793  itsgxy[minyp][x]+itsgxy[yp][x]+itsgxy[ypp][x]+itsgxy[yppp][x];
794  for (y=minyp+1; y<maxym; y++){
795  myp=y+3; mym=y-4; my1=y-1;
796  isgxy[y][x]=isgxy[my1][x]-itsgxy[mym][x]+itsgxy[myp][x];
797  }
798  }
799 
800  for (x=minxp; x<maxxm; x++){
801  isgxx[minyp][x]=itsgxx[ymmm][x]+itsgxx[ymm][x]+itsgxx[ym][x]+
802  itsgxx[minyp][x]+itsgxx[yp][x]+itsgxx[ypp][x]+itsgxx[yppp][x];
803  for (y=minyp+1; y<maxym; y++){
804  myp=y+3; mym=y-4; my1=y-1;
805  isgxx[y][x]=isgxx[my1][x]-itsgxx[mym][x]+itsgxx[myp][x];
806  }
807  }
808 
809  for (x=minxp; x<maxxm; x++){
810  isgyy[minyp][x]=itsgyy[ymmm][x]+itsgyy[ymm][x]+itsgyy[ym][x]+
811  itsgyy[minyp][x]+itsgyy[yp][x]+itsgyy[ypp][x]+itsgyy[yppp][x];
812  for (y=minyp+1; y<maxym; y++){
813  myp=y+3; mym=y-4; my1=y-1;
814  isgyy[y][x]=isgyy[my1][x]-itsgyy[mym][x]+itsgyy[myp][x];
815  }
816  }
817  }
818 #ifdef TIME_MEASURE
819  timer.Stop();
820  cerr << "second convolution took "<<timer.GetRealTime()<<" us\n";
821  timer.Reset();
822  timer.Start();
823 #endif
824 
825 #ifdef BIAS_DEBUG
826  if (DebugLevelIsSet(D_ST_WRITE_DEBUG_IM)){
827  _WriteDebugImage(sgxx, "sgxx-7x7.mip");
828  _WriteDebugImage(sgxy, "sgxy-7x7.mip");
829  _WriteDebugImage(sgyy, "sgyy-7x7.mip");
830 
831  CalcStructureTensorValid(gx, gy, sgxx, sgxy, sgyy);
832  _WriteDebugImage(sgxx, "sgxx-st.mip");
833  _WriteDebugImage(sgxy, "sgxy-st.mip");
834  _WriteDebugImage(sgyy, "sgyy-st.mip");
835  }
836 #endif
837 
838 #ifdef TIME_MEASURE
839  timer.Stop();
840  cerr << "CalcStructureTensor7x7 without _GradientProducts took "<<timer.GetRealTime()<<" us\n";
841 #endif
842  return 0;
843 }
844 
845 //////////////////////////////////////////////////////////////////////////
846 // protected functions
847 //////////////////////////////////////////////////////////////////////////
848 
849 template <class InputStorageType, class OutputStorageType>
851 _AllocInternalMem(const int width, const int height)
852 {
853  if (!_gxx.IsEmpty()) _gxx.Release();
854  _gxx.Init(width, height, 1);
855  if (!_gxy.IsEmpty()) _gxy.Release();
856  _gxy.Init(width, height, 1);
857  if (!_gyy.IsEmpty()) _gyy.Release();
858  _gyy.Init(width, height, 1);
859  if (!_tsgxx.IsEmpty()) _tsgxx.Release();
860  _tsgxx.Init(width, height, 1);
861  if (!_tsgxy.IsEmpty()) _tsgxy.Release();
862  _tsgxy.Init(width, height, 1);
863  if (!_tsgyy.IsEmpty()) _tsgyy.Release();
864  _tsgyy.Init(width, height, 1);
865 
866  _InternalMemWidth=width;
867  _InternalMemHeight=height;
868 }
869 
870 template <class InputStorageType, class OutputStorageType>
873 {
874 #ifdef BIAS_DEBUG
875  //if (ImageIO::Save(name, im)!=0){
876  if (ImageIO::Save(name, im)!=0){
877  BIASERR("error wrting "<<name);
878  }
879 #endif
880 }
881 
882 
883 template <class InputStorageType, class OutputStorageType>
886  const Image<OutputStorageType>& gy)
887 {
888  const int width=gx.GetWidth();
889  const int height=gx.GetHeight();
890  if (width>_InternalMemWidth || height>_InternalMemHeight){
891  //_DeleteInternalMem();
892  _AllocInternalMem(width, height);
893  }
894 #ifdef BIAS_DEBUG
895  BIASASSERT(gx.SamePixelAndChannelCount(gy));
896  BIASASSERT(gx.NotBiggerPixelAndSameChannelCount(_gxx));
897  BIASASSERT(gx.NotBiggerPixelAndSameChannelCount(_gxy));
898  BIASASSERT(gx.NotBiggerPixelAndSameChannelCount(_gyy));
899 #endif
900 
901  int xtlx, xtly, xbrx, xbry, ytlx, ytly, ybrx, ybry;
902  gx.GetROI()->GetCorners(xtlx, xtly, xbrx, xbry);
903  gy.GetROI()->GetCorners(ytlx, ytly, ybrx, ybry);
904 
905  const int minx=((xtlx<ytlx)?(ytlx):(xtlx));
906  const int maxx=((xbrx<ybrx)?(xbrx):(ybrx));
907  const int miny=((xtly<ytly)?(ytly):(xtly));
908  const int maxy=((xbrx<ybry)?(xbry):(ybry));
909 
910  _gxx.GetROI()->SetCorners(minx, miny, maxx, maxy);
911  _gxy.GetROI()->SetCorners(minx, miny, maxx, maxy);
912  _gyy.GetROI()->SetCorners(minx, miny, maxx, maxy);
913 
914  /*
915  register OutputStorageType **igx=gx.GetImageDataArray();
916  register OutputStorageType **igy=gy.GetImageDataArray();
917  register OutputStorageType **igxx=_gxx.GetImageDataArray();
918  register OutputStorageType **igxy=_gxy.GetImageDataArray();
919  register OutputStorageType **igyy=_gyy.GetImageDataArray();
920  register int x, y;
921  register OutputStorageType mgx, mgy;
922  for (y=miny; y<maxy; y++){
923  for (x=minx; x<maxx; x++){
924  mgx=igx[y][x];
925  mgy=igy[y][x];
926  igxx[y][x]=mgx*mgx;
927  igxy[y][x]=mgx*mgy;
928  igyy[y][x]=mgy*mgy;
929  }
930  }
931 
932  */
933 
934 
935  register const int step=width-(maxx-minx);
936  const int offset=miny*width+minx;
937  register const OutputStorageType *igx=gx.GetImageData()+offset;
938  register const OutputStorageType *igy=gy.GetImageData()+offset;
939  register OutputStorageType *igxx=_gxx.GetImageData()+offset;
940  register OutputStorageType *igxy=_gxy.GetImageData()+offset;
941  register OutputStorageType *igyy=_gyy.GetImageData()+offset;
942 
943  register const OutputStorageType *end=gx.GetImageData()+(maxy-1)*width+maxx;
944  const OutputStorageType *lend=gx.GetImageData()+miny*width+maxx;
945 
946  while (igx<end){
947  while (igx<lend){
948  *igxx++=*igx * *igx;
949  *igyy++=*igy * *igy;
950  *igxy++=*igx++ * *igy++;
951  }
952  igx+=step;
953  igy+=step;
954  igxx+=step;
955  igxy+=step;
956  igyy+=step;
957  lend+=width;
958  }
959 
960 
961 }
962 
963 template <class InputStorageType, class OutputStorageType>
965 GetBordersValid_(int &border_x, int &border_y) const
966 {
967  _grad->GetBorders(border_x, border_y);
968  border_x += _HalfWinSize;
969  border_y += _HalfWinSize;
970 }
971 
972 //////////////////////////////////////////////////////////////////////////
973 // instantiation
974 //////////////////////////////////////////////////////////////////////////
975 namespace BIAS{
976 #define FILTER_INSTANTIATION_CLASS StructureTensor
977 #define FILTER_INSTANTIATION_NO_UNSIGNED_OUTPUT
978 #include "Filterinst.hh"
979 }
void Release()
reimplemented from ImageBase
Definition: Image.cpp:1579
gradient calculation with separated gauss masks
int SetCorners(unsigned UpperLeftX, unsigned UpperLeftY, unsigned LowerRightX, unsigned LowerRightY)
Sets a rectangular region of interest.
Definition: ROI.cpp:287
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
virtual int FilterInt(const Image< InputStorageType > &src, Image< OutputStorageType > &dst)
virtual void GetBordersValid_(int &border_x, int &border_y) const
bool IsEmpty() const
check if ImageData_ points to allocated image buffer or not
Definition: ImageBase.hh:245
virtual int FilterFloat(const Image< InputStorageType > &src, Image< OutputStorageType > &dst)
void _AllocInternalMem(const int width, const int height)
allocates _gxx, _gyy, _gxy, _tsgxx, _tsgxy, _tsgyy and sets InternalMemWidth and _InternalMemHeight ...
unsigned int GetWidth() const
Definition: ImageBase.hh:312
basic class for structure tensor calculation
int CalcStructureTensor(const Image< InputStorageType > &src, Image< OutputStorageType > &sgxx, Image< OutputStorageType > &sgxy, Image< OutputStorageType > &sgyy)
Calculates the gradients uses internal gradient class and calls CalcStructureTensor with these gradie...
int CalcStructureTensor3x3(const Image< OutputStorageType > &gx, const Image< OutputStorageType > &gy, Image< OutputStorageType > &sgxx, Image< OutputStorageType > &sgxy, Image< OutputStorageType > &sgyy)
Uses _HalfWinSize=1 Destination images must be initialized.
void _GradientProducts(const Image< OutputStorageType > &gx, const Image< OutputStorageType > &gy)
fills up _gxx, _gxy, _gyy from gx and gy
FilterNTo2N< InputStorageType, OutputStorageType > * _grad
pointer to gradient filter actually used:
ROI * GetROI()
Returns a pointer to the roi object.
Definition: ImageBase.hh:615
unsigned int GetHeight() const
Definition: ImageBase.hh:319
virtual int Filter(const Image< InputStorageType > &src, Image< OutputStorageType > &dst)
calls the apropriate non virtual member function returns an 3 channel image where ...
int CalcStructureTensorValid(const Image< OutputStorageType > &gx, const Image< OutputStorageType > &gy, Image< OutputStorageType > &sgxx, Image< OutputStorageType > &sgxy, Image< OutputStorageType > &sgyy)
For any _HalfWinSize.
bool SamePixelAndChannelCount(const ImageBase &Image) const
checks if data area has same &quot;size&quot; as Image of other type
Definition: ImageBase.hh:73
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
int CalcStructureTensor7x7(const Image< OutputStorageType > &gx, const Image< OutputStorageType > &gy, Image< OutputStorageType > &sgxx, Image< OutputStorageType > &sgxy, Image< OutputStorageType > &sgyy)
Uses _HalfWinSize=3 Destination images must be initialized.
double GetRealTime() const
return real time (=wall time clock) in usec JW For Win32: real-time is measured differently from user...
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
const StorageType * GetImageData() const
overloaded GetImageData() from ImageBase
Definition: Image.hh:137
void _WriteDebugImage(Image< OutputStorageType > &im, std::string name)
for debugging only
bool NotBiggerPixelAndSameChannelCount(const ImageBase &Image) const
checks if data area has bigger or the same &quot;size&quot; as Image of other type
Definition: ImageBase.hh:78
class TimeMeasure contains functions for timing real time and cpu time.
Definition: TimeMeasure.hh:111
int CalcStructureTensor5x5(const Image< OutputStorageType > &gx, const Image< OutputStorageType > &gy, Image< OutputStorageType > &sgxx, Image< OutputStorageType > &sgxy, Image< OutputStorageType > &sgyy)
Uses _HalfWinSize=2 Destination images must be initialized.
base class for simple n-&gt;3n filter implementations
Definition: FilterNTo3N.hh:43
const StorageType ** GetImageDataArray() const
overloaded GetImageDataArray() from ImageBase
Definition: Image.hh:153