iLab Neuromorphic Robotics Toolkit  0.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
SeparableFilterImpl.H
Go to the documentation of this file.
1 /*! @file
2  @author Randolph Voorhies
3  @copyright GNU Public License (GPL v3)
4  @section License
5  @verbatim
6  // ////////////////////////////////////////////////////////////////////////
7  // The iLab Neuromorphic Robotics Toolkit (NRT) //
8  // Copyright 2010-2012 by the University of Southern California (USC) //
9  // and the iLab at USC. //
10  // //
11  // iLab - University of Southern California //
12  // Hedco Neurociences Building, Room HNB-10 //
13  // Los Angeles, Ca 90089-2520 - USA //
14  // //
15  // See http://ilab.usc.edu for information about this project. //
16  // ////////////////////////////////////////////////////////////////////////
17  // This file is part of The iLab Neuromorphic Robotics Toolkit. //
18  // //
19  // The iLab Neuromorphic Robotics Toolkit is free software: you can //
20  // redistribute it and/or modify it under the terms of the GNU General //
21  // Public License as published by the Free Software Foundation, either //
22  // version 3 of the License, or (at your option) any later version. //
23  // //
24  // The iLab Neuromorphic Robotics Toolkit is distributed in the hope //
25  // that it will be useful, but WITHOUT ANY WARRANTY; without even the //
26  // implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR //
27  // PURPOSE. See the GNU General Public License for more details. //
28  // //
29  // You should have received a copy of the GNU General Public License //
30  // along with The iLab Neuromorphic Robotics Toolkit. If not, see //
31  // <http://www.gnu.org/licenses/>. //
32  // ////////////////////////////////////////////////////////////////////////
33  @endverbatim */
34 
35 
36 template <NRT_PROMOTE_PIX_NO_DEFAULT>
39  nrt::ConvolutionBoundaryStrategy const boundary)
40 {
41  // promote the source image to float if necessary, so that we do the promotion only once for all, rather than many
42  // times as we access the pixels of the image; if no promotion is necessary, "source" will just point to the original
43  // data of "src" through the copy-on-write/ref-counting behavior of Image:
44  nrt::Image<DestType> const source = nrt::Image<DestType>(src);
45 
46  assert(hFilt.width() == 1 || hFilt.height() == 1);
47 
48  int const hfs = int(std::max(hFilt.width(), hFilt.height()));
49  if (hfs == 0) return source; // no filter
50 
51  int const w = src.width(), h = src.height();
52 
54  typename nrt::Image<DestType>::const_iterator sptr = source.const_begin();
55  typename nrt::Image<DestType>::iterator dptr = result.begin();
56 
57  const int hfs2 = (hfs - 1) / 2;
58 
59  // flip the filter to accelerate convolution:
60  nrt::Image<DestType> hFilter(hFilt.dims());
61  typename nrt::Image<DestType>::iterator hFilterIt = hFilter.begin();
62  typename nrt::Image<DestType>::const_iterator hFiltIt = hFilt.const_begin();
63  DestType sumh(0.0F);
64 
65  for (int x = 0; x < hfs; ++x)
66  {
67  sumh += hFiltIt[x]; hFilterIt[hfs-1-x] = hFiltIt[x];
68  }
69 
70  // *** horizontal pass ***
71  if (w < hfs) // special function for very small images
72  {
73  switch (boundary)
74  {
76  {
77  for (int j = 0; j < h; ++j)
78  for (int i = 0; i < w; ++i)
79  {
80  DestType val = DestType(0);
81  for (int k = 0; k < hfs; ++k)
82  if (i + k - hfs2 >= 0 && i + k - hfs2 < w)
83  {
84  val += sptr[k - hfs2] * hFilterIt[k];
85  }
86  *dptr++ = val;
87  ++ sptr;
88  }
89  //TODO modify the way to access the value on the ptr
90  //Iterator it = hFilterIt; Iterator it2 = sptr - hfs2;
91  // for (int k = 0; k < hfs; k ++)
92  // if (i + k - hfs2 >= 0 && i + k - hfs2 < w)
93  // {
94  // val += *it2 * *it;
95  // }
96  // *dptr++ = val; sptr ++; it ++; it2 ++;
97  //this way is faster to compute
98  }
99  break;
101  {
102  for (int j = 0; j < h; ++j)
103  for (int i = 0; i < w; ++i)
104  {
105  DestType val(0); DestType sum(0);
106  for (int k = 0; k < hfs; ++k)
107  if (i + k - hfs2 >= 0 && i + k - hfs2 < w)
108  {
109  val += sptr[k - hfs2] * hFilterIt[k];
110  sum += hFilterIt[k];
111  }
112  *dptr++ = val * sumh / sum; ++ sptr;
113  }
114  }
115  break;
117  {
118  for (int j = 0; j < h; ++j)
119  for (int i = 0; i < w; ++i)
120  {
121  DestType val = DestType(0);
122  for (int k = 0; k < hfs; ++k)
123  if (i + k - hfs2 < 0)
124  {
125  val += sptr[-i] * hFilterIt[k];
126  }
127  else if (i + k - hfs2 >= w)
128  {
129  val += sptr[w-1-i] * hFilterIt[k];
130  }
131  else
132  {
133  val += sptr[k - hfs2] * hFilterIt[k];
134  }
135  *dptr++ = val; ++sptr;
136  }
137  }
138  break;
139  }
140  }
141 
142  else // function for reasonably large images
143  for (int jj = 0; jj < h; ++jj)
144  {
145  // leftmost points:
146  switch (boundary)
147  {
149  {
150  for (int j = hfs2; j > 0; --j)
151  {
152  DestType val = DestType(0);
153  for (int k = j; k < hfs; ++k)
154  {
155  val += sptr[k - j] * hFilterIt[k];
156  }
157  *dptr++ = val;
158  }
159  }
160  break;
162  {
163  for (int j = hfs2; j > 0; --j)
164  {
165  DestType val(0); DestType sum(0);
166  for (int k = j; k < hfs; ++k)
167  {
168  val += sptr[k - j] * hFilterIt[k];
169  sum += hFilterIt[k];
170  }
171  *dptr++ = val * sumh / sum;
172  }
173  }
174  break;
176  {
177  for (int j = hfs2; j > 0; --j)
178  {
179  DestType val = DestType(0);
180  for (int k = 0; k < j; ++k)
181  {
182  val += sptr[0] * hFilterIt[k];
183  }
184  for (int k = j; k < hfs; ++k)
185  {
186  val += sptr[k - j] * hFilterIt[k];
187  }
188  *dptr++ = val;
189  }
190  }
191  break;
192  }
193  // bulk (far from the borders):
194  for (int i = 0; i < w - hfs + 1; ++i)
195  {
196  DestType val = DestType(0);
197  for (int k = 0; k < hfs; ++k) val += sptr[k] * hFilterIt[k];
198  *dptr++ = val; ++sptr;
199  }
200  // rightmost points:
201  switch (boundary)
202  {
204  {
205  for (int j = hfs - 1; j > hfs2; --j)
206  {
207  DestType val = DestType(0);
208  for (int k = 0; k < j; ++k)
209  {
210  val += sptr[k] * hFilterIt[k];
211  }
212  *dptr++ = val;
213  ++sptr;
214  }
215  }
216  break;
218  {
219  for (int j = hfs - 1; j > hfs2; --j)
220  {
221  DestType val(0); DestType sum(0);
222  for (int k = 0; k < j; ++k)
223  {
224  val += sptr[k] * hFilterIt[k];
225  sum += hFilterIt[k];
226  }
227  *dptr++ = val * sumh / sum;
228  ++sptr;
229  }
230  }
231  break;
233  {
234  for (int j = hfs - 1; j > hfs2; --j)
235  {
236  DestType val = DestType(0);
237  for (int k = 0; k < j; ++k)
238  {
239  val += sptr[k] * hFilterIt[k];
240  }
241  for (int k = j; k < hfs; ++k)
242  {
243  val += sptr[j-1] * hFilterIt[k];
244  }
245  *dptr++ = val;
246  ++sptr;
247  }
248  }
249  break;
250  }
251  sptr += hfs2; // sptr back to same as dptr (start of next line)
252  }
253  return result;
254 }
255 
256 // ######################################################################
257 template <NRT_PROMOTE_PIX_NO_DEFAULT>
260  nrt::ConvolutionBoundaryStrategy const boundary)
261 {
262  // promote the source image to float if necessary, so that we do the promotion only once for all, rather than many
263  // times as we access the pixels of the image; if no promotion is necessary, "source" will just point to the original
264  // data of "src" through the copy-on-write/ref-counting behavior of Image:
265  nrt::Image<DestType> const source = nrt::Image<DestType>(src);
266  assert(vFilt.width() == 1 || vFilt.height() == 1);
267  int const vfs = int(std::max(vFilt.width(), vFilt.height()));
268  if (vfs == 0) return source; // no filter
269 
270  int const w = int(src.width()), h = int(src.height());
271 
273  typename nrt::Image<DestType>::const_iterator sptr = source.const_begin();
274  typename nrt::Image<DestType>::iterator dptr = result.begin();
275 
276  const int vfs2 = (vfs - 1) / 2;
277 
278  // flip the filter to accelerate convolution:
279  nrt::Image<DestType> vFilter(vFilt.dims());
280  typename nrt::Image<DestType>::iterator vFilterIt = vFilter.begin();
281  typename nrt::Image<DestType>::const_iterator vFiltIt = vFilt.const_begin();
282  DestType sumv(0);
283  for (int x = 0; x < vfs; ++x)
284  { sumv += vFiltIt[x]; vFilterIt[vfs-1-x] = vFiltIt[x]; }
285 
286  int ww[vfs];
287  for (int i = 0; i < vfs; i ++) ww[i] = w * i; // speedup precompute
288 
289  if (h < vfs) // special function for very small images
290  {
291  switch (boundary)
292  {
294  {
295  for (int j = 0; j < h; ++j)
296  for (int i = 0; i < w; ++i)
297  {
298  DestType val = DestType(0);
299  for (int k = 0; k < vfs; ++k)
300  if (j + k - vfs2 >= 0 && j + k - vfs2 < h)
301  {
302  val += sptr[w * (k - vfs2)] * vFilterIt[k];
303  }
304  *dptr++ = val; ++sptr;
305  }
306  }
307  break;
309  {
310  for (int j = 0; j < h; ++j)
311  for (int i = 0; i < w; ++i)
312  {
313  DestType val(0); DestType sum(0);
314  for (int k = 0; k < vfs; ++k)
315  if (j + k - vfs2 >= 0 && j + k - vfs2 < h)
316  {
317  val += sptr[w * (k - vfs2)] * vFilterIt[k];
318  sum += vFilterIt[k];
319  }
320  *dptr++ = val * sumv / sum; ++sptr;
321  }
322  }
323  break;
325  {
326  for (int j = 0; j < h; ++j)
327  for (int i = 0; i < w; ++i)
328  {
329  DestType val = DestType(0);
330  for (int k = 0; k < vfs; ++k)
331  if (j + k - vfs2 < 0)
332  {
333  val += sptr[w * (-j)] * vFilterIt[k];
334  }
335  else if (j + k - vfs2 >= h)
336  {
337  val += sptr[w * (h-1-j)] * vFilterIt[k];
338  }
339  else
340  {
341  val += sptr[w * (k - vfs2)] * vFilterIt[k];
342  }
343  *dptr++ = val; ++sptr;
344  }
345  }
346  break;
347  }
348  }
349  else // function for reasonably large images
350  {
351  // top points:
352  switch (boundary)
353  {
355  {
356  for (int j = vfs2; j > 0; --j)
357  {
358  for (int i = 0; i < w; ++i) // scan all points of given horiz
359  {
360  DestType val = DestType(0);
361  for (int k = j; k < vfs; ++k)
362  {
363  val += sptr[ww[k - j]] * vFilterIt[k];
364  }
365  *dptr++ = val;
366  ++sptr;
367  }
368  sptr -= w; // back to top-left corner
369  }
370  }
371  break;
373  {
374  for (int j = vfs2; j > 0; --j)
375  {
376  for (int i = 0; i < w; ++i) // scan all points of given horiz
377  {
378  DestType val(0); DestType sum(0);
379  for (int k = j; k < vfs; ++k)
380  {
381  val += sptr[ww[k - j]] * vFilterIt[k];
382  sum += vFilterIt[k];
383  }
384  *dptr++ = val * sumv / sum;
385  ++sptr;
386  }
387  sptr -= w; // back to top-left corner
388  }
389  }
390  break;
392  {
393  for (int j = vfs2; j > 0; --j)
394  {
395  for (int i = 0; i < w; ++i) // scan all points of given horiz
396  {
397  DestType val = DestType(0);
398  for (int k = 0; k < j; ++k)
399  {
400  val += sptr[ww[0]] * vFilterIt[k];
401  }
402  for (int k = j; k < vfs; ++k)
403  {
404  val += sptr[ww[k - j]] * vFilterIt[k];
405  }
406  *dptr++ = val;
407  sptr++;
408  }
409  sptr -= w; // back to top-left corner
410  }
411  }
412  break;
413  }
414  // bulk (far from edges):
415  for (int j = 0; j < h - vfs + 1; ++j)
416  for (int i = 0; i < w; ++i)
417  {
418  DestType val = DestType(0);
419  for (int k = 0; k < vfs; ++k) val += sptr[ww[k]] * vFilterIt[k];
420  *dptr++ = val; ++sptr;
421  }
422  // bottommost points:
423  switch (boundary)
424  {
426  {
427  for (int j = vfs - 1; j > vfs2; --j)
428  for (int i = 0; i < w; ++i)
429  {
430  DestType val = DestType(0);
431  for (int k = 0; k < j; ++k)
432  {
433  val += sptr[ww[k]] * vFilterIt[k];
434  }
435  *dptr++ = val;
436  ++sptr;
437  }
438  }
439  break;
441  {
442  for (int j = vfs - 1; j > vfs2; --j)
443  for (int i = 0; i < w; ++i)
444  {
445  DestType val(0); DestType sum(0);
446  for (int k = 0; k < j; ++k)
447  {
448  val += sptr[ww[k]] * vFilterIt[k];
449  sum += vFilterIt[k];
450  }
451  *dptr++ = val * sumv / sum;
452  ++sptr;
453  }
454  }
455  break;
457  {
458  for (int j = vfs - 1; j > vfs2; --j)
459  for (int i = 0; i < w; ++i)
460  {
461  DestType val = DestType(0);
462  for (int k = 0; k < j; ++k)
463  {
464  val += sptr[ww[k]] * vFilterIt[k];
465  }
466  for (int k = j; k < vfs; ++k)
467  {
468  val += sptr[ww[j-1]] * vFilterIt[k];
469  }
470  *dptr++ = val;
471  ++sptr;
472  }
473  }
474  break;
475  }
476  }
477  return result;
478 }
479 
480 // ######################################################################
481 template <NRT_PROMOTE_PIX_NO_DEFAULT>
484  nrt::Image<DestType> const vKernel, nrt::ConvolutionBoundaryStrategy const boundary)
485 {
487 
488  assert(hKernel.width() == 1 || hKernel.height() == 1);
489  assert(vKernel.width() == 1 || vKernel.height() == 1);
490 
491  return xFilter<DestType>(yFilter<DestType>(result, vKernel, boundary), hKernel, boundary);
492 }
493 
494 // ######################################################################
495 template <NRT_PROMOTE_PIX_NO_DEFAULT>
498  nrt::ConvolutionBoundaryStrategy const boundary)
499 {
500  nrt::Image<DestType> const result = nrt::Image<DestType>(src);
501 
502  if(kernel.width() > 0 && kernel.height() == 1)
503  return xFilter<DestType>(result, kernel, boundary);
504 
505  if(kernel.width() == 1 && kernel.height() > 0)
506  return yFilter<DestType>(result, kernel, boundary);
507 
508  if(kernel.width() == 0 && kernel.height() == 0)
509  return result;
510 
511  NRT_FATAL(kernel.dims() << " kernel passed to separableFilter (kernels must be 1-D)");
512  return result; // never happens
513 }