Basic Image AlgorithmS Library  2.8.0
 All Classes Namespaces Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
GaussThreshold.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 "GaussThreshold.hh"
26 #include <Base/Image/ImageConvert.hh>
27 #include <Base/Image/ImageIO.hh>
28 
29 using namespace BIAS;
30 using namespace std;
31 
32 #define MIN_NORM 0.6
33 
34 //////////////////////////////////////////////////////////////////////////
35 // implementation
36 //////////////////////////////////////////////////////////////////////////
37 
38 template <class InputStorageType, class OutputStorageType>
41  : Convolution<InputStorageType, OutputStorageType>()
42 {
43  _GaussSigma = 0.7;
44  _GaussRatio = 0.01;
45  _LastRatio = _LastSigma = -1.0;
46  _GaussIO = NULL;
47  _GaussOO = NULL;
48  _Threshold = 0.0;
49 }
50 
51 template <class InputStorageType, class OutputStorageType>
54  : _GaussSigma(other._GaussSigma),
55  _GaussRatio(other._GaussRatio),
56  _LastSigma(other._LastSigma),
57  _LastRatio(other._LastRatio)
58 {
59  _Threshold = 0.0;
60 }
61 
62 template <class StorageType, class KernelType>
64 {
65  if (_GaussIO!=NULL) delete _GaussIO;
66  if (_GaussOO!=NULL) delete _GaussOO;
67 }
68 
69 ////////////////////////////////////////////////////////////////////////////
70 
71 template <class InputStorageType, class OutputStorageType>
74 {
75  BIASCDOUT(D_FILTERBASE_CALLSTACK, "Gauss::Filter\n");
76  if ((_GaussSigma != _LastSigma)||(_GaussRatio != _LastRatio))
77  _CalculateKernels(_GaussSigma, _GaussRatio);
78 
79  if (_fm.IsSize1x1()){ // do not filter if filter mask is 1x1
80  return ImageConvert::ConvertST(src, dst, dst.GetStorageType());
81  }
82  if ((src.GetChannelCount()==1) && (this->GetBorderHandling() == TBH_valid))
83  return Filter7x7GreyIgnoreBelowThreshold(src, dst, _Threshold);
84  BIASERR("Not yet implemented");
85  BIASABORT;
87 }
88 
89 
90 ////////////////////////////////////////////////////////////////////////////
91 
92 template <class InputStorageType, class OutputStorageType>
94 _CalculateKernels(double Sigma, double Ratio)
95 {
96  const int hws = (int)(floor(sqrt(-2.0*Sigma*Sigma*log(Ratio))));
97  const int size = (hws<<1)+1;
98 
99  Vector<FM_FLOAT> H(size), V(size);
100  BIASCDOUT(D_CONV_KERNEL, "Sigma = "<<Sigma<<"\tRatio = "<<Ratio<<endl);
101  BIASCDOUT(D_CONV_KERNEL, "size = "<<size<<"\thws = "<<hws<<endl);
102 
103  double sum=0.0;
104  int i, ip, im;
105  if (Sigma==0.0){
106  H[0]=1.0; V[0]=1.0;
107  } else {
108  double fac=1.0/(2.0*Sigma*Sigma);
109  register double g;
110  // calculate the filter coefficients as gauss function
111  for (i=0; i<=hws; i++){
112  ip = hws+i;
113  im = hws-i;
114  g = expf((float)(-(double)(i*i)*fac));
115  V[ip] = V[im] = H[im] = H[ip] = (FM_FLOAT)g;
116  //BIASCDOUT(D_CONV_KERNEL, i<<"v:"<<this->_VVec[im]<<" "<<endl);
117  sum += g*2.0;
118  }
119  sum -= V[hws]; // added twice
120  BIASCDOUT(D_CONV_KERNEL, "sum = "<<sum<<endl);
121 
122  // normalize filter, so that the sum is 1.0
123  sum = 1.0/sum;
124  BIASCDOUT(D_CONV_KERNEL, "sum = "<<sum<<endl);
125  for (i=0; i<=hws; i++){
126  ip=hws+i;
127  im=hws-i;
128  H[ip] = H[im] = V[ip] = V[im] = V[im] * (FM_FLOAT)sum;
129  BIASCDOUT(D_CONV_KERNEL, H[im]<<" "<<endl);
130  }
131  }
132  _fm.Init(H,V);
133  // compute the int approximation
134  int thebits = _fm.ComputeIntPrecisionBits(sizeof(InputStorageType)*8,
135  sizeof(CONV_INT)*8-1);
136  _fm.CreateIntFilter(thebits,thebits,thebits);
137 
138  _LastSigma = Sigma;
139  _LastRatio = Ratio;
140 
141 #ifdef BIAS_DEBUG
142  if (this->DebugLevelIsSet(D_CONV_KERNEL)){
143  sum = 0.0;
144  for (i=-hws; i<=hws; i++){
145  sum += H[i+hws];
146  BIASCDOUT(D_CONV_KERNEL, H[i+hws]<<" ");
147  }
148  BIASCDOUT(D_CONV_KERNEL, " sum = "<<sum<<endl);
149  }
150 #endif
151 
152  //#define COMPARE
153 #ifdef COMPARE
154  fprintf(stdout, "gauss: \n");
155  for (i=-hws; i<=hws; i++){
156  fprintf(stdout, "%f ",HVec[i+hws]);
157  }
158  fprintf(stdout, "\n");
159 #endif
160 }
161 
162 template <class InputStorageType, class OutputStorageType>
166  const double thresh)
167 {
168  BIASCDOUT(D_FILTERBASE_CALLSTACK, "GaussThreshold::Filter7x7\n");
169  SetHalfWinSize(3, false);
170  if ((_GaussSigma != _LastSigma)||(_GaussRatio != _LastRatio))
171  _CalculateKernels(_GaussSigma, _GaussRatio);
172 
173  const FM_FLOAT h0 = (const FM_FLOAT)((*_fm.GetSepfh())[3]);
174  const FM_FLOAT h1 = (const FM_FLOAT)((*_fm.GetSepfh())[2]);
175  const FM_FLOAT h2 = (const FM_FLOAT)((*_fm.GetSepfh())[1]);
176  const FM_FLOAT h3 = (const FM_FLOAT)((*_fm.GetSepfh())[0]);
177  const FM_FLOAT v0 = (const FM_FLOAT)((*_fm.GetSepfv())[3]);
178  const FM_FLOAT v1 = (const FM_FLOAT)((*_fm.GetSepfv())[2]);
179  const FM_FLOAT v2 = (const FM_FLOAT)((*_fm.GetSepfv())[1]);
180  const FM_FLOAT v3 = (const FM_FLOAT)((*_fm.GetSepfv())[0]);
181 
182 #ifdef BIAS_DEBUG
183  if (src.GetChannelCount()!=1) {
184  BIASERR("image with multiple channels. stopping.");
185  BIASABORT;
186  }
187 #endif
188  if(!dst.SamePixelAndChannelCount(src)) {
189  if (!dst.IsEmpty()) dst.Release();
190  dst.Init(src.GetWidth(), src.GetHeight(), 1);
191  }
192  if(!_tmpFloat.SamePixelAndChannelCount(src)) {
193  if (!_tmpFloat.IsEmpty()) _tmpFloat.Release();
194  _tmpFloat.Init(src.GetWidth(),src.GetHeight(),1,
196  }
197 
198  // horizontal convolution ignoring wrap-arounds
199  const InputStorageType* pS = src.GetImageData();
200  CONV_FLOAT* pD = _tmpFloat.GetImageData();
201  // copy first three bytes
202  *pD++ = CONV_FLOAT(thresh);
203  pS++;
204  *pD++ = CONV_FLOAT(thresh);
205  pS++;
206  *pD++ = CONV_FLOAT(thresh);
207 
208  const int width = src.GetWidth(), height = src.GetHeight();
209  const int width2 = 2*width, width3=3*width;
210  const InputStorageType* pEndH = src.GetImageData() + width*height-3;
211 
212  while (++pS<pEndH) {
213  register CONV_FLOAT normalization = 0;
214  register CONV_FLOAT accumulator = 0;
215  if (pS[-3]>thresh) {
216  normalization += h3;
217  accumulator += h3*CONV_FLOAT(pS[-3]);
218  }
219  if (pS[-2]>thresh) {
220  normalization += h2;
221  accumulator += h2*CONV_FLOAT(pS[-2]);
222  }
223  if (pS[-1]>thresh) {
224  normalization += h1;
225  accumulator += h1*CONV_FLOAT(pS[-1]);
226  }
227  if (pS[0]>thresh) {
228  normalization += h0;
229  accumulator += h0*CONV_FLOAT(pS[0]);
230  }
231  if (pS[1]>thresh) {
232  normalization += h1;
233  accumulator += h1*CONV_FLOAT(pS[1]);
234  }
235 
236  if (pS[2]>thresh) {
237  normalization += h2;
238  accumulator += h2*CONV_FLOAT(pS[2]);
239  }
240  if (pS[3]>thresh) {
241  normalization += h3;
242  accumulator += h3*CONV_FLOAT(pS[3]);
243  }
244  if (normalization>MIN_NORM)
245  *pD++ = CONV_FLOAT(accumulator / normalization);
246  else
247  *pD++ = CONV_FLOAT(thresh);
248  }
249  // copy last three bytes
250  *pD++ = CONV_FLOAT(thresh);pS++;
251  *pD++ = CONV_FLOAT(thresh);pS++;
252  *pD = CONV_FLOAT(thresh);pS++;
253 
254  const CONV_FLOAT cfthresh = CONV_FLOAT(thresh);
255 
256  // vertical convolution
257  const CONV_FLOAT* pEndV = _tmpFloat.GetImageData() + width*(height-3);
258  CONV_FLOAT* pI = _tmpFloat.GetImageData();
259  OutputStorageType* pO = dst.GetImageData();
260  const CONV_FLOAT* pEndFirst3 = _tmpFloat.GetImageData() + 3*width;
261  const CONV_FLOAT* pEndLast3 = _tmpFloat.GetImageData() + width*height;
262  // copy first three lines
263  while (pI<pEndFirst3) {
264  *pO++ = OutputStorageType(thresh);
265  pI++;
266  }
267  // convolve main
268  pI--;
269  while (++pI<pEndV) {
270  register CONV_FLOAT normalization = 0;
271  register CONV_FLOAT accumulator = 0;
272 
273  if (pI[-width3]>cfthresh) {
274  normalization += v3;
275  accumulator += v3*CONV_FLOAT(pI[-width3]);
276  }
277  if (pI[-width2]>cfthresh) {
278  normalization += v2;
279  accumulator += v2*CONV_FLOAT(pI[-width2]);
280  }
281  if (pI[-width]>cfthresh) {
282  normalization += v1;
283  accumulator += v1*CONV_FLOAT(pI[-width]);
284  }
285  if (pI[0]>cfthresh) {
286  normalization += v0;
287  accumulator += v0*CONV_FLOAT(pI[0]);
288  }
289  if (pI[width]>cfthresh) {
290  normalization += v1;
291  accumulator += v1*CONV_FLOAT(pI[width]);
292  }
293 
294  if (pI[width2]>cfthresh) {
295  normalization += v2;
296  accumulator += v2*CONV_FLOAT(pI[width2]);
297  }
298  if (pI[width3]>cfthresh) {
299  normalization += v3;
300  accumulator += v3*CONV_FLOAT(pI[width3]);
301  }
302 
303  if (normalization>MIN_NORM)
304  *pO++ = OutputStorageType(accumulator / normalization);
305  else
306  *pO++ = OutputStorageType(thresh);
307 
308  }
309  // copy last three lines
310  while (pI<pEndLast3) {
311  *pO++ = OutputStorageType(thresh);
312  pI++;
313  }
314 
315  // update roi
316  unsigned minx, miny, maxx, maxy;
317  int hws = 3;
318 
319  src.GetROI()->GetCorners(minx, miny, maxx, maxy);
320  minx += hws;
321  miny += hws;
322  maxx -= hws;
323  maxy -= hws;
324  if (minx>maxx) minx = maxx = (minx+maxx)/2;
325  if (miny>maxy) miny = maxy = (miny+maxy)/2;
326  dst.GetROI()->SetCorners(minx, miny, maxx, maxy);
327 
328  return 0;
329 }
330 
331 
332 template <class InputStorageType, class OutputStorageType>
336  const double thresh)
337 {
338  BIASCDOUT(D_FILTERBASE_CALLSTACK, "GaussThreshold::Filter7x7 only below threshold\n");
339  SetHalfWinSize(3, false);
340  if ((_GaussSigma != _LastSigma)||(_GaussRatio != _LastRatio))
341  _CalculateKernels(_GaussSigma, _GaussRatio);
342 
343  const FM_FLOAT h0 = (const FM_FLOAT)((*_fm.GetSepfh())[3]);
344  const FM_FLOAT h1 = (const FM_FLOAT)((*_fm.GetSepfh())[2]);
345  const FM_FLOAT h2 = (const FM_FLOAT)((*_fm.GetSepfh())[1]);
346  const FM_FLOAT h3 = (const FM_FLOAT)((*_fm.GetSepfh())[0]);
347  const FM_FLOAT v0 = (const FM_FLOAT)((*_fm.GetSepfv())[3]);
348  const FM_FLOAT v1 = (const FM_FLOAT)((*_fm.GetSepfv())[2]);
349  const FM_FLOAT v2 = (const FM_FLOAT)((*_fm.GetSepfv())[1]);
350  const FM_FLOAT v3 = (const FM_FLOAT)((*_fm.GetSepfv())[0]);
351 
352 #ifdef BIAS_DEBUG
353  if (src.GetChannelCount()!=1) {
354  BIASERR("image with multiple channels. stopping.");
355  BIASABORT;
356  }
357 #endif
358  if(!dst.SamePixelAndChannelCount(src)) {
359  if (!dst.IsEmpty()) dst.Release();
360  dst.Init(src.GetWidth(), src.GetHeight(), 1);
361  }
362  if(!_tmpFloat.SamePixelAndChannelCount(src)) {
363  if (!_tmpFloat.IsEmpty()) _tmpFloat.Release();
364  _tmpFloat.Init(src.GetWidth(),src.GetHeight(),1,
366  }
367 
368  // horizontal convolution ignoring wrap-arounds
369  const InputStorageType* pS = src.GetImageData();
370  CONV_FLOAT* pD = _tmpFloat.GetImageData();
371  // copy first three bytes
372  *pD++ = CONV_FLOAT(thresh);
373  pS++;
374  *pD++ = CONV_FLOAT(thresh);
375  pS++;
376  *pD++ = CONV_FLOAT(thresh);
377 
378  const int width = src.GetWidth(), height = src.GetHeight();
379  const int width2 = 2*width, width3=3*width;
380  const InputStorageType* pEndH = src.GetImageData() + width*height-3;
381 
382  while (++pS<pEndH) {
383  register CONV_FLOAT normalization = 0;
384  register CONV_FLOAT accumulator = 0;
385  if (pS[-3]>thresh) {
386  normalization += h3;
387  accumulator += h3*CONV_FLOAT(pS[-3]);
388  }
389  if (pS[-2]>thresh) {
390  normalization += h2;
391  accumulator += h2*CONV_FLOAT(pS[-2]);
392  }
393  if (pS[-1]>thresh) {
394  normalization += h1;
395  accumulator += h1*CONV_FLOAT(pS[-1]);
396  }
397  if (pS[0]>thresh) {
398  normalization += h0;
399  accumulator += h0*CONV_FLOAT(pS[0]);
400  }
401  if (pS[1]>thresh) {
402  normalization += h1;
403  accumulator += h1*CONV_FLOAT(pS[1]);
404  }
405 
406  if (pS[2]>thresh) {
407  normalization += h2;
408  accumulator += h2*CONV_FLOAT(pS[2]);
409  }
410  if (pS[3]>thresh) {
411  normalization += h3;
412  accumulator += h3*CONV_FLOAT(pS[3]);
413  }
414 
415  // only fill pixels with a value below threshold
416  if (*pS <= OutputStorageType(thresh) && normalization>MIN_NORM)
417  *pD++ = CONV_FLOAT(accumulator / normalization);
418  else
419  *pD++ = CONV_FLOAT(*pS);
420 
421  }
422  // copy last three bytes
423  *pD++ = CONV_FLOAT(thresh);pS++;
424  *pD++ = CONV_FLOAT(thresh);pS++;
425  *pD = CONV_FLOAT(thresh);pS++;
426 
427  const CONV_FLOAT cfthresh = CONV_FLOAT(thresh);
428 
429  // vertical convolution
430  const CONV_FLOAT* pEndV = _tmpFloat.GetImageData() + width*(height-3);
431  CONV_FLOAT* pI = _tmpFloat.GetImageData();
432  OutputStorageType* pO = dst.GetImageData();
433  const CONV_FLOAT* pEndFirst3 = _tmpFloat.GetImageData() + 3*width;
434  const CONV_FLOAT* pEndLast3 = _tmpFloat.GetImageData() + width*height;
435 
436  // copy first three lines of original image to output image
437  while (pI<pEndFirst3) {
438  *pO++ = OutputStorageType(thresh);
439  pI++;
440  }
441  // convolve main
442  pI--;
443  while (++pI<pEndV) {
444  register CONV_FLOAT normalization = 0;
445  register CONV_FLOAT accumulator = 0;
446 
447  if (pI[-width3]>cfthresh) {
448  normalization += v3;
449  accumulator += v3*CONV_FLOAT(pI[-width3]);
450  }
451  if (pI[-width2]>cfthresh) {
452  normalization += v2;
453  accumulator += v2*CONV_FLOAT(pI[-width2]);
454  }
455  if (pI[-width]>cfthresh) {
456  normalization += v1;
457  accumulator += v1*CONV_FLOAT(pI[-width]);
458  }
459  if (pI[0]>cfthresh) {
460  normalization += v0;
461  accumulator += v0*CONV_FLOAT(pI[0]);
462  }
463  if (pI[width]>cfthresh) {
464  normalization += v1;
465  accumulator += v1*CONV_FLOAT(pI[width]);
466  }
467 
468  if (pI[width2]>cfthresh) {
469  normalization += v2;
470  accumulator += v2*CONV_FLOAT(pI[width2]);
471  }
472  if (pI[width3]>cfthresh) {
473  normalization += v3;
474  accumulator += v3*CONV_FLOAT(pI[width3]);
475  }
476 
477  // only fill pixels which are below threshold
478  if (*pI <= OutputStorageType(thresh) && normalization > MIN_NORM)
479  *pO++ = OutputStorageType(accumulator / normalization);
480  else
481  *pO++ = OutputStorageType(*pI);
482  }
483 
484  // copy last three lines of original image to output image
485  while (pI < pEndLast3) {
486  *pO++ = OutputStorageType(thresh);
487  pI++;
488  }
489 
490  // update roi
491  unsigned minx, miny, maxx, maxy;
492  int hws = 3;
493 
494  src.GetROI()->GetCorners(minx, miny, maxx, maxy);
495  minx += hws;
496  miny += hws;
497  maxx -= hws;
498  maxy -= hws;
499  if (minx>maxx) minx = maxx = (minx+maxx)/2;
500  if (miny>maxy) miny = maxy = (miny+maxy)/2;
501  dst.GetROI()->SetCorners(minx, miny, maxx, maxy);
502 
503  return 0;
504 }
505 
506 
507 //////////////////////////////////////////////////////////////////////////
508 // instantiation
509 //////////////////////////////////////////////////////////////////////////
510 namespace BIAS{
511 #define FILTER_INSTANTIATION_CLASS GaussThreshold
512 #include "Filterinst.hh"
513 }
void Release()
reimplemented from ImageBase
Definition: Image.cpp:1579
generic convolution class.
Definition: Convolution.hh:66
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
bool IsEmpty() const
check if ImageData_ points to allocated image buffer or not
Definition: ImageBase.hh:245
bool IsInterleaved() const
Definition: ImageBase.hh:491
double _GaussSigma
the parameter sigma of gaussian kernel
double _GaussRatio
minimum ratio of 1D kernel center and border, ignore smaller entries
virtual int Filter(const Image< InputStorageType > &src, Image< OutputStorageType > &dst)
decides on the image types which FilterFunction should be called
float image storage type
Definition: ImageBase.hh:118
unsigned int GetWidth() const
Definition: ImageBase.hh:312
Gauss< InputStorageType, OutputStorageType > * _GaussIO
ROI * GetROI()
Returns a pointer to the roi object.
Definition: ImageBase.hh:615
static int ConvertST(const BIAS::ImageBase &source, BIAS::ImageBase &dest, ImageBase::EStorageType targetST)
Function to convert the storage type of images e.g.
unsigned int GetChannelCount() const
returns the number of Color channels, e.g.
Definition: ImageBase.hh:382
smoothing with gaussian kernel using a threshold
unsigned int GetHeight() const
Definition: ImageBase.hh:319
bool SamePixelAndChannelCount(const ImageBase &Image) const
checks if data area has same &quot;size&quot; as Image of other type
Definition: ImageBase.hh:73
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 _CalculateKernels(double Sigma, double Ratio)
calculates the kernel
Gauss< OutputStorageType, OutputStorageType > * _GaussOO
to allow for iterated gauss convolution saving the intermediate image we need a different instance...
virtual int Filter(const Image< InputStorageType > &src, Image< OutputStorageType > &dst)
sets gauss kernel if params changed and calls convolution or fast grey implementation if possible ...
enum EStorageType GetStorageType() const
Definition: ImageBase.hh:414
int Filter7x7GreyOnlyBelowThreshold(const Image< InputStorageType > &src, Image< OutputStorageType > &dst, const double thresh)
7x7 gauss filtering, values below threshold are ignored and only pixels below threshold are filed ...
int Filter7x7GreyIgnoreBelowThreshold(const Image< InputStorageType > &src, Image< OutputStorageType > &dst, const double thresh)
7x7 gauss filtering, values below threshold are ignored useful for depth map filtering ...