iLab Neuromorphic Robotics Toolkit  0.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MiscFiltersImpl.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> inline
38 nrt::centerSurround(nrt::Image<PixType> const center, nrt::Image<PixType> const surround, bool const rectify)
39 {
40  int const lw = center.width(), lh = center.height();
41  int const sw = surround.width(), sh = surround.height();
42 
43  if (sw > lw || sh > lh)
44  NRT_FATAL("Center (" << center.dims() << ") must be larger than surround (" << surround.dims() << ")");
45 
46  int scalex = lw / sw, remx = lw - 1 - (lw % sw);
47  int scaley = lh / sh, remy = lh - 1 - (lh % sh);
48 
49  // result has the size of the larger image:
51 
52  // scan large image and subtract corresponding pixel from small image:
53  int ci = 0, cj = 0;
54  typename nrt::Image<PixType>::const_iterator lptr = center.const_begin();
55  typename nrt::Image<PixType>::const_iterator sptr = surround.const_begin();
56  typename nrt::Image<DestType>::iterator dptr = result.begin();
57 
58  if (rectify == true) // compute abs(hires - lowres):
59  {
60  for (int j = 0; j < lh; ++j)
61  {
62  for (int i = 0; i < lw; ++i)
63  {
64  for(size_t chanIdx = 0; chanIdx < nrt::pixel_traits<PixType>::num_channels; ++chanIdx)
65  dptr->channels[chanIdx] = std::abs(lptr->channels[chanIdx] - sptr->channels[chanIdx]);
66 
67  dptr++; lptr++;
68 
69  if ((++ci) == scalex && i != remx) { ci = 0; ++sptr; }
70  }
71  if (ci) { ci = 0; ++sptr; } // in case the reduction is not round
72  if ((++cj) == scaley && j != remy) cj = 0; else sptr -= sw;
73  }
74  }
75  else // compute hires - lowres, clamped to 0:
76  {
77  for (int j = 0; j < lh; ++j)
78  {
79  for (int i = 0; i < lw; ++i)
80  {
81  for(size_t chanIdx = 0; chanIdx < nrt::pixel_traits<PixType>::num_channels; ++chanIdx)
82  dptr->channels[chanIdx] =
83  std::max(typename nrt::pixel_traits<PixType>::pod_type(0),
84  lptr->channels[chanIdx] - sptr->channels[chanIdx]);
85 
86  dptr++; lptr++;
87 
88  if ((++ci) == scalex && i != remx) { ci = 0; ++sptr; }
89  }
90  if (ci) { ci = 0; ++sptr; } // in case the reduction is not round
91  if ((++cj) == scaley && j != remy) cj = 0; else sptr -= sw;
92  }
93  }
94 
95  // attenuate borders:
96  int bord = std::max(result.dims().width(), result.dims().height());
97  inplaceAttenuateBorders(result, bord / 20);
98 
99  return result;
100 }
101 
102 // ######################################################################
103 namespace miscfilters_details
104 {
105  pthread_once_t trigtab_init_once = PTHREAD_ONCE_INIT;
106 
107  int const TRIGTAB_SIZE = 256;
108 
109  double sintab[TRIGTAB_SIZE];
110  double costab[TRIGTAB_SIZE];
111 
112  void trigtab_init()
113  {
114  for (int i = 0; i < 256; ++i)
115  {
116  double const arg = (2.0*M_PI*i) / 256.0;
117 
118 #if defined(HAVE_ASM_FSINCOS)
119  asm("fsincos"
120  :"=t"(costab[i]), "=u"(sintab[i])
121  :"0"(arg));
122 #elif defined(HAVE_SINCOS)
123  sincos(arg, &sintab[i], &costab[i]);
124 #else
125  sintab[i] = sin(arg);
126  costab[i] = cos(arg);
127 #endif
128  }
129  }
130 }
131 
132 // ######################################################################
133 template <NRT_PROMOTE_PIX_NO_DEFAULT> inline
135 nrt::orientedFilter(nrt::Image<PixType> src, float const k, float const theta,
136  float const intensity, bool const fastMath)
137 {
138  double kx = double(k) * cos((theta + 90.0) * M_PI / 180.0);
139  double ky = double(k) * sin((theta + 90.0) * M_PI / 180.0);
140 
142 
143  nrt::Image<DestType> const source = nrt::Image<DestType>(src);
144  typename nrt::Image<DestType>::const_iterator sptr = source.const_begin();
145  typename nrt::Image<DestType>::iterator reptr = re.begin(), imptr = im.begin();
146 
147  // (x,y) = (0,0) at center of image:
148  int w2l = src.width() / 2, w2r = src.width() - w2l;
149  int h2l = src.height() / 2, h2r = src.height() - h2l;
150 
151  if (fastMath)
152  {
153  pthread_once(&miscfilters_details::trigtab_init_once, &miscfilters_details::trigtab_init);
154 
155  double const kx2 = (256.0 * kx) / (2.0*M_PI);
156  double const ky2 = (256.0 * ky) / (2.0*M_PI);
157 
158  for (int j = -h2l; j < h2r; ++j)
159  for (int i = -w2l; i < w2r; ++i)
160  {
161  double const arg2 = kx2 * i + ky2 * j;
162  int idx = int(arg2) % 256;
163  if (idx < 0) idx += 256;
164  DestType const val = (*sptr++) * intensity;
165 
166  double const sinarg = miscfilters_details::sintab[idx];
167  double const cosarg = miscfilters_details::costab[idx];
168 
169  *reptr++ = DestType(val * cosarg);
170  *imptr++ = DestType(val * sinarg);
171  }
172  }
173  else
174  {
175  for (int j = -h2l; j < h2r; ++j)
176  for (int i = -w2l; i < w2r; ++i)
177  {
178  double const arg = kx * i + ky * j;
179  DestType const val = (*sptr++) * intensity;
180 
181 #if defined(NRT_HAVE_ASM_FSINCOS)
182  // If we the asm "fsincos" instruction is available, then
183  // that will be the fastest way to compute paired sin and
184  // cos calls.
185  double sinarg, cosarg;
186  asm("fsincos"
187  :"=t"(cosarg), "=u"(sinarg)
188  :"0"(arg));
189 
190 #elif defined(NRT_HAVE_SINCOS)
191  // Otherwise we use a sincos() library function if it is
192  // available as an extension (it's present at least in GNU
193  // libc, maybe on other systems as well) then it is faster
194  // (~30% overall speedup of orientedFilter()) to compute
195  // both sin+cos at once rather than doing separate calls
196  // -- this allows the use of the x87 assembly instruction
197  // fsincos (see /usr/include/bits/mathinline.h for the
198  // glibc sincos()).
199  double sinarg, cosarg;
200  sincos(arg, &sinarg, &cosarg);
201 
202 #else
203  // Otherwise we resort to explicit separate sin() and
204  // cos() calls.
205  double sinarg = sin(arg);
206  double cosarg = cos(arg);
207 #endif
208 
209  *reptr++ = DestType(val * cosarg);
210  *imptr++ = DestType(val * sinarg);
211  }
212  }
213 
214  re = nrt::lowPass9<promo>(re);
215  im = nrt::lowPass9<promo>(im);
216 
217  typedef typename nrt::pixel_traits<DestType>::pod_type pod_type;
218 
219  return nrt::channel_transform(re, im, [](pod_type const& c1, pod_type const& c2)
220  {
221  return pod_type(sqrt(c1*c1 + c2*c2));
222  });
223 }
224 
225 // ######################################################################
226 template <NRT_PROMOTE_PIX_NO_DEFAULT>
227 std::tuple<nrt::Image<DestType>, nrt::Image<DestType> > const
228 nrt::sobelFilter3(Image<PixType> const src, ConvolutionBoundaryStrategy strategy)
229 {
230  // Compute the x gradient
231  nrt::Image<DestType> xGradVfilter = { {DestType(1)},
232  {DestType(2)},
233  {DestType(1)} };
234  nrt::Image<DestType> xGradHfilter = {{ DestType(1) , DestType(0), DestType(-1) }};
235  nrt::Image<DestType> const Gx = separableFilter<promo>(src, xGradHfilter, xGradVfilter, strategy);
236 
237  // Compute the y gradient
238  nrt::Image<DestType> yGradVfilter = xGradHfilter; yGradVfilter.reshape(xGradHfilter.height(), xGradHfilter.width());
239  nrt::Image<DestType> yGradHfilter = yGradVfilter; yGradHfilter.reshape(xGradVfilter.height(), xGradVfilter.width());
240  nrt::Image<DestType> const Gy = separableFilter(src, yGradHfilter, yGradVfilter, strategy);
241 
242  // Compute the magnitude and orientation
243  nrt::Image<DestType> magnitude(src.dims(), nrt::ImageInitPolicy::None);
244  nrt::Image<DestType> orientation(src.dims(), nrt::ImageInitPolicy::None);
245  typename nrt::Image<DestType>::pod_iterator magIt = magnitude.pod_begin();
246  typename nrt::Image<DestType>::pod_iterator oriIt = orientation.pod_begin();
250 
251  while(gyIt != gyEnd)
252  {
253  typename nrt::pixel_traits<DestType>::pod_type const gx = *gxIt;
254  typename nrt::pixel_traits<DestType>::pod_type const gy = *gyIt;
255  *magIt = sqrt(gx*gx + gy*gy);
256  *oriIt = atan2(gy,gx);
257 
258  ++magIt; ++oriIt; ++gxIt; ++gyIt;
259  }
260 
261  return std::make_tuple(magnitude, orientation);
262 }
263 
264 
265 // ######################################################################
266 template <NRT_PROMOTE_PIX_NO_DEFAULT>
267 std::tuple<nrt::Image<DestType>, nrt::Image<DestType>> const
268 nrt::sobelFilter5(Image<PixType> const src, ConvolutionBoundaryStrategy strategy)
269 {
270  // Compute the x gradient
271  nrt::Image<DestType> xGradVfilter = { {DestType(1)},
272  {DestType(4)},
273  {DestType(6)},
274  {DestType(4)},
275  {DestType(1)} };
276  nrt::Image<DestType> xGradHfilter = {{ DestType(1), DestType(2), DestType(0), DestType(-2), DestType(-1) }};
277  nrt::Image<DestType> const Gx = separableFilter<promo>(src, xGradHfilter, xGradVfilter, strategy);
278 
279  // Compute the y gradient
280  nrt::Image<DestType> yGradVfilter = xGradHfilter; yGradVfilter.reshape(xGradHfilter.height(), xGradHfilter.width());
281  nrt::Image<DestType> yGradHfilter = xGradVfilter; yGradHfilter.reshape(xGradVfilter.height(), xGradVfilter.width());
282  nrt::Image<DestType> const Gy = separableFilter(src, yGradHfilter, yGradVfilter, strategy);
283 
284  // Compute the magnitude and orientation
285  nrt::Image<DestType> magnitude(src.dims(), nrt::ImageInitPolicy::None);
286  nrt::Image<DestType> orientation(src.dims(), nrt::ImageInitPolicy::None);
287  typename nrt::Image<DestType>::pod_iterator magIt = magnitude.pod_begin();
288  typename nrt::Image<DestType>::pod_iterator oriIt = orientation.pod_begin();
292 
293  while(gyIt != gyEnd)
294  {
295  typename nrt::pixel_traits<DestType>::pod_type const gx = *gxIt;
296  typename nrt::pixel_traits<DestType>::pod_type const gy = *gyIt;
297  *magIt = sqrt(gx*gx + gy*gy);
298  *oriIt = atan2(gy,gx);
299 
300  ++magIt; ++oriIt; ++gxIt; ++gyIt;
301  }
302 
303  return std::make_tuple(magnitude, orientation);
304 }
305 
306 // ######################################################################
307 template <NRT_PROMOTE_PIX_NO_DEFAULT>
310 {
312  int const imw = src.width(), imh = src.height();
314  typename nrt::Image<DestType>::iterator ii = res.begin(), prev_line = ii;
315 
316  *ii++ = *in++;
317 
318  for (int x = 1; x < imw; ++x) { *ii = *in + *(ii-1); ++ii; ++in; }
319  for (int y = 1; y < imh; ++y)
320  {
321  DestType s(0);
322  for (int x = 0; x < imw; ++x) { s += *in; *ii = s + *prev_line; ++ii; ++in; ++prev_line; }
323  }
324 
325  return res;
326 }
327 
328 // ######################################################################
329 template <NRT_PROMOTE_PIX_NO_DEFAULT>
332 {
334  int const imw = src.width(), imh = src.height();
336  typename nrt::Image<DestType>::iterator ii = res.begin(), prev_line = ii;
337 
338  *ii++ = (*in) * (*in); ++in;
339 
340  for (int x = 1; x < imw; ++x) { *ii = (*in) * (*in) + *(ii-1); ++ii; ++in; }
341  for (int y = 1; y < imh; ++y)
342  {
343  DestType s(0);
344  for (int x = 0; x < imw; ++x) { s += (*in) * (*in); *ii = s + *prev_line; ++ii; ++in; ++prev_line; }
345  }
346 
347  return res;
348 }
349 
350 // ######################################################################
351 namespace nrt {
352  namespace detail {
353  namespace miscfilters {
354  // helper function for median3x()/median3y()/median3()
355  template <class T> inline
356  T median3pod(T a, T b, T c)
357  {
358  if (b < a) std::swap(a, b);
359  if (c < b) std::swap(b, c);
360  if (b < a) std::swap(a, b);
361 
362  return b;
363  }
364 
365  // helper function for median3x()/median3y()/median3()
366  template <class T, template <typename> class PixType> inline
367  PixType<T> median3scalar(PixType<T> const & a, PixType<T> const & b, PixType<T> const & c)
368  {
369  PixType<T> result;
370  T const * aptr = a.begin(); T const * bptr = b.begin(); T const * cptr = c.begin();
371  T * rptr = result.begin();
372 
373  for (size_t i = 0; i < PixType<T>::numChannels; ++i) *rptr++ = median3pod(*aptr++, *bptr++, *cptr++);
374  return result;
375  }
376  }
377  }
378 }
379 
380 // ######################################################################
381 template <class PixType> inline
383 {
384  if (in.width() < 3) return in;
385  int const w = in.width(), h = in.height();
387 
389  typename nrt::Image<PixType>::iterator dptr = result.begin();
390 
391  for (int y = 0; y < h; ++y)
392  {
393  // first point is untouched:
394  *dptr++ = *sptr++;
395 
396  // Bulk of the row:
397  for (int x = 1; x < w-1; ++x)
398  {
399  *dptr++ = nrt::detail::miscfilters::median3scalar(sptr[-1], sptr[0], sptr[1]);
400  ++sptr;
401  }
402 
403  // Last point is untouched:
404  *dptr++ = *sptr++;
405  }
406 
407  return result;
408 }
409 
410 // ######################################################################
411 template <class PixType> inline
413 {
414  if (in.height() < 3) return in;
415  int const w = in.width(), h = in.height();
417 
419  typename nrt::Image<PixType>::iterator dptr = result.begin();
420 
421  // First row is untouched:
422  for (int x = 0; x < w; ++x) *dptr++ = *sptr++;
423 
424  // Bulk:
425  for (int y = 1; y < h-1; ++y)
426  for (int x = 0; x < w; ++x)
427  {
428  *dptr++ = nrt::detail::miscfilters::median3scalar(sptr[-w], sptr[0], sptr[w]);
429  ++sptr;
430  }
431 
432  // Last row is untouched:
433  for (int x = 0; x < w; ++x) *dptr++ = *sptr++;
434 
435  return result;
436 }
437 
438 // ######################################################################
439 template <class PixType> inline
440 nrt::Image<PixType> nrt::median3(nrt::Image<PixType> const & in, bool go_x, bool go_y)
441 {
442  nrt::Image<PixType> result = in;
443  if (go_x) result = median3x(result);
444  if (go_y) result = median3y(result);
445  return result;
446 }
447 
448 
449