iLab Neuromorphic Robotics Toolkit  0.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RescaleImpl.H
Go to the documentation of this file.
1 /*! @file
2  @author David J. Berg <dberg@usc.edu>
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 #ifndef INCLUDE_NRT_IMAGEPROC_RESHAPING_DETAILS_RESCALEIMPL_H
37 #define INCLUDE_NRT_IMAGEPROC_RESHAPING_DETAILS_RESCALEIMPL_H
38 
39 // ######################################################################
40 template <NRT_PROMOTE_PIX_NO_DEFAULT> inline
42 {
43  // Is the input empty?
44  if (src.size() == 0) return nrt::Image<DestType>(dims, DestType(0));
45 
46  // if we dont have dims, then empty image
47  if (dims.empty()) return nrt::Image<DestType>();
48 
49  // Promote the source:
50  nrt::Image<DestType> const source = nrt::Image<DestType>(src);
51 
52  // check if same size already:
53  if (dims == source.dims()) return source;
54 
55  int32 const orig_w = src.width(), orig_h = src.height();
56  int32 const new_w = dims.width(), new_h = dims.height();
57 
58  float const sw = float(orig_w) / float(new_w), sh = float(orig_h) / float(new_h);
59 
61  typename nrt::Image<DestType>::iterator dptr = result.begin();
62  typename nrt::Image<DestType>::const_iterator sptr = source.const_begin();
63 
64  int32 const oh1 = orig_h - 1, ow1 = orig_w - 1;
65 
66  // code inspired from one of the Graphics Gems books:
67  /* (1) (x,y) are the original coords corresponding to scaled coords (i,j)
68  (2) (x0,y0) are the greatest lower bound integral coords from (x,y)
69  (3) (x1,y1) are the least upper bound integral coords from (x,y)
70  (4) d00, d10, d01, d11 are the values of the original image at the corners of the rect (x0,y0),(x1,y1)
71  (5) the value in the scaled image is computed from bilinear interpolation among d00,d10,d01,d11 */
72  for (int32 j = 0; j < new_h; ++j)
73  {
74  float const y = j * sh;
75  int32 const y0 = int32(y);
76  int32 const y1 = std::min(y0 + 1, oh1);
77  float const fy = y - float(y0);
78 
79  int32 const wy0 = orig_w * y0;
80  int32 const wy1 = orig_w * y1;
81 
82  for (int32 i = 0; i < new_w; ++i)
83  {
84  float const x = i * sw;
85  int32 const x0 = int32(x);
86  int32 const x1 = std::min(x0 + 1, ow1);
87  float const fx = x - float(x0);
88 
90  d00( sptr[x0 + wy0] ), d10( sptr[x1 + wy0] ),
91  d01( sptr[x0 + wy1] ), d11( sptr[x1 + wy1] ),
92  dx0( d00 + (d10 - d00) * fx ),
93  dx1( d01 + (d11 - d01) * fx );
94 
95  DestType val(dx0 + (dx1 - dx0) * fy ); // no need to clamp
96 
97  *dptr++ = val;
98  }
99  }
100  return result;
101 }
102 
103 // ######################################################################
104 namespace
105 {
106  template <typename promo>
107  struct RescaleBilinearVisitor : public boost::static_visitor<nrt::GenericImage const>
108  {
109  inline RescaleBilinearVisitor(nrt::Dims<nrt::int32> const dims_) : dims(dims_) { }
110 
111  template <class PixType> inline
112  nrt::GenericImage const operator()(nrt::Image<PixType> const im) const
113  { return nrt::GenericImage(nrt::rescaleBilinear<promo>(im, dims)); }
114 
115  private:
117  };
118 }
119 
120 template <typename promo> inline
122 { return src.apply_visitor(RescaleBilinearVisitor<promo>(dims)); }
123 
124 
125 // ######################################################################
126 template <NRT_PROMOTE_PIX_NO_DEFAULT> inline
128 {
129  assert(factor > 0);
130 
131  if (src.width() <= 1) return src; // do not go smaller than 1 pixel wide
132 
133  nrt::Image<DestType> const source = nrt::Image<DestType>(src);
134  int const h = source.height();
135  int w2 = source.width() / factor;
136  int skip = source.width() % factor;
137  if (w2 == 0)
138  {
139  w2 = 1;
140  skip = source.width() - factor;
141  }
142 
144 
145  typename nrt::Image<DestType>::const_iterator sptr = source.const_begin();
146  typename nrt::Image<DestType>::iterator dptr = result.begin();
147 
148  for (int j = 0; j < h; ++j)
149  {
150  for (int i = 0; i < w2; ++i)
151  {
152  *dptr++ = *sptr; // copy one point
153  sptr += factor; // skip a few points
154  }
155  sptr += skip;
156  }
157 
158  return result;
159 }
160 
161 // ######################################################################
162 template <NRT_PROMOTE_PIX_NO_DEFAULT> inline
164 {
165  assert(factor > 0);
166 
167  if (src.height() <= 1) return src; // do not go smaller than 1 pixel high
168 
169  nrt::Image<DestType> const source = nrt::Image<DestType>(src);
170  int const w = source.width();
171  int h2 = source.height() / factor;
172  if (h2 == 0) h2 = 1;
173 
175 
176  typename nrt::Image<DestType>::const_iterator sptr = source.const_begin();
177  typename nrt::Image<DestType>::iterator dptr = result.begin();
178  int skip = w * factor;
179 
180  for (int j = 0; j < h2; ++j)
181  {
182  memcpy(dptr, sptr, w*sizeof(*dptr));
183  dptr += w; sptr += skip;
184  }
185 
186  return result;
187 }
188 
189 // ######################################################################
190 template <NRT_PROMOTE_PIX_NO_DEFAULT> inline
191 nrt::Image<DestType> nrt::decimateXY(nrt::Image<PixType> const & src, int const xfactor, int const yfactor_raw)
192 {
193  int const yfactor = yfactor_raw >= 0 ? yfactor_raw : xfactor;
194 
196 
197  // do not go smaller than 1x1:
198  if (src.width() <= 1 && src.height() <= 1) return source;
199 
200  if (source.width() == 1) // only thinout vertic
201  {
202  int h2 = source.height() / yfactor;
203 
205 
206  for (int j = 0; j < h2; ++j)
207  result(0,j) = source(0, j * yfactor);
208 
209  return result;
210  }
211  if (source.height() == 1)
212  {
213  int w2 = source.width() / xfactor;
214 
216 
217  for (int i = 0; i < w2; ++i)
218  result(i, 0) = source(i * xfactor, 0);
219 
220  return result;
221  }
222 
223  int const w2 = src.width() / xfactor;
224  int const h2 = src.height() / yfactor;
225 
227 
228  typename nrt::Image<DestType>::const_iterator sptr = source.const_begin();
229  typename nrt::Image<DestType>::iterator dptr = result.begin();
230  int const skip = source.width() % xfactor + source.width() * (yfactor - 1);
231 
232  for (int j = 0; j < h2; ++j)
233  {
234  for (int i = 0; i < w2; ++i)
235  {
236  *dptr++ = *sptr; // copy one pixel
237  sptr += xfactor; // skip some pixels
238  }
239  sptr += skip; // skip to start of next line
240  }
241 
242  return result;
243 }
244 
245 // ######################################################################
246 template <NRT_PROMOTE_PIX_NO_DEFAULT> inline
248  int const filterWidth)
249 {
250  nrt::int32 const new_w = dims.width();
251  nrt::int32 const new_h = dims.height();
252 
253  if (int(src.width()) == new_w && int(src.height()) == new_h) return src;
254 
256  assert(source.width() / new_w > 1 && source.height() / new_h > 1);
257 
258  int const wdepth = int(0.5+log(double(source.width() / new_w)) / M_LN2);
259  int const hdepth = int(0.5+log(double(source.height() / new_h)) / M_LN2);
260 
261  // TODO: Fix this by downsizing to max(width,height) and then bilinearly interpolating
262  if (wdepth != hdepth) NRT_FATAL("arrays must have same proportions");
263 
264  nrt::Image<DestType> result = src;
265 
266  for (int i = 0; i < wdepth; ++i)
267  {
268  result = decimateX<promo>(nrt::lowPassX<float>(filterWidth, result));
269  result = decimateY<promo>(nrt::lowPassY<float>(filterWidth, result));
270  }
271 
272  return result;
273 }
274 
275 // ######################################################################
276 template <NRT_PROMOTE_PIX_NO_DEFAULT> inline
278 {
279  assert(factor > 0);
280 
281  // if factor == 1 then we don't decimate:
282  if (factor == 1) return lowPass5x<promo>(src);
283 
284  // now factor is guaranteed to be >= 2
285  int const w = src.width(), h = src.height();
286  int w2 = w / factor;
287  if (w2 == 0) w2 = 1; // do not go smaller than 1 pixel
288 
289  // promote the source image to float if necessary, so that we do the promotion only once for all, rather than many
290  // times as we access the pixels of the image; if no promotion is necessary, "source" will just point to the original
291  // data of "src" through the copy-on-write/ref-counting behavior of Image:
292  nrt::Image<DestType> const source = nrt::Image<DestType>(src);
293  if (w < 2) return source; // nothing to smooth
294 
296  typename nrt::Image<DestType>::const_iterator sptr = source.const_begin();
297  typename nrt::Image<DestType>::iterator dptr = result.begin();
298 
299  if (w == 2) //////////////////////////////////////////////////
300  for (int j = 0; j < h; ++ j)
301  {
302  // leftmost point [ (6^) 4 ] / 10
303  *dptr++ = sptr[0] * (6.0F / 10.0F) + sptr[1] * (4.0F / 10.0F);
304 
305  sptr += 2; // sptr back to same position as dptr
306  }
307  else if (w == 3) //////////////////////////////////////////////////
308  for (int j = 0; j < h; ++ j)
309  {
310  // need the left most point in any case
311  // leftmost point [ (6^) 4 1 ] / 11
312  *dptr++ = sptr[0] * (6.0F / 11.0F) +
313  sptr[1] * (4.0F / 11.0F) +
314  sptr[2] * (1.0F / 11.0F);
315 
316  // since factor >= 2 we don't need any other point.
317 
318  sptr += 3; // sptr back to same position as dptr
319  }
320  else ////////////////////////////// general case for width() >= 4
321  // *** unfolded version (all particular cases treated) for
322  // max speed. do not use access overloading since can
323  // access the C array directly. the bet is: float
324  // computations are faster than int (true for most float
325  // copros), and int2float convert is fast -> do not have to
326  // divide by the coeff sum since we directly use float
327  // coeffs.
328  // notations: in () is the position of dest ptr, and ^ is src ptr
329  // ########## horizontal pass
330  for (int j = 0; j < h; ++ j)
331  {
332  int i1 = 0, i2 = 0;
333  typename nrt::Image<DestType>::const_iterator sptr2 = sptr;
334 
335 
336  // leftmost point [ (6^) 4 1 ] / 11
337  *dptr++ = sptr2[0] * (6.0F / 11.0F) +
338  sptr2[1] * (4.0F / 11.0F) +
339  sptr2[2] * (1.0F / 11.0F);
340  ++i2;
341  i1 += factor;
342 
343  // skip second point since factor >= 2:
344  sptr2 += (factor-2);
345 
346  // rest of the line except last 2 points [ 1^ 4 (6) 4 1 ] / 16.0
347  while ((i1 < (w-2)) && (i2 < w2))
348  {
349  *dptr++ = (sptr2[0] + sptr2[4]) * (1.0F / 16.0F) +
350  (sptr2[1] + sptr2[3]) * (4.0F / 16.0F) +
351  sptr2[2] * (6.0F / 16.0F);
352  i1 += factor; sptr2 += factor;
353  ++i2;
354  }
355 
356  // need special case for second to last point?
357  if ((i2 < w2) && (i1 == (w-2)))
358  {
359  sptr2 = sptr + w - 4;
360  // before last point [ 1^ 4 (6) 4 ] / 15
361  *dptr++ = sptr2[0] * (1.0F / 15.0F) +
362  (sptr2[1] + sptr2[3]) * (4.0F / 15.0F) +
363  sptr2[2] * (6.0F / 15.0F);
364  i1 += factor;
365  ++i2;
366  }
367 
368  // need special case for last point?
369  if ((i2 < w2) && (i1 == (w-1)))
370  {
371  sptr2 = sptr + w - 3;
372  // last point [ 1^ 4 (6) ] / 11
373  *dptr++ = sptr2[0] * (1.0F / 11.0F) +
374  sptr2[1] * (4.0F / 11.0F) +
375  sptr2[2] * (6.0F / 11.0F);
376  ++i2;
377  }
378  sptr += w;
379  }
380  return result;
381 }
382 
383 // ######################################################################
384 template <NRT_PROMOTE_PIX_NO_DEFAULT> inline
386 {
387  assert(factor > 0);
388 
389  // if factor == 1 then we don't decimate:
390  if (factor == 1) return lowPass5y<promo>(src);
391 
392  // now factor is guaranteed to be >= 2
393  int const w = src.width(), h = src.height();
394  int h2 = h / factor;
395  if (h2 == 0) h2 = 1; // do not go smaller than 1 pixel
396 
397  // promote the source image to float if necessary, so that we do the promotion only once for all, rather than many
398  // times as we access the pixels of the image; if no promotion is necessary, "source" will just point to the original
399  // data of "src" through the copy-on-write/ref-counting behavior of Image:
400  nrt::Image<DestType> const source = nrt::Image<DestType>(src);
401  if (h < 2) return source; // nothing to smooth
402 
404  typename nrt::Image<DestType>::const_iterator sptr = source.const_begin();
405  typename nrt::Image<DestType>::iterator dptr = result.begin();
406 
407  // ########## vertical pass (even though we scan horiz for speedup)
408  int const w2 = w * 2, w3 = w * 3, w4 = w * 4; // speedup
409 
410  if (h == 2) //////////////////////////////////////////////////
411  {
412  // topmost points ( [ (6^) 4 ] / 10 )^T
413  for (int i = 0; i < w; ++ i)
414  {
415  *dptr++ = sptr[0] * (6.0F / 10.0F) + sptr[w] * (4.0F / 10.0F);
416  ++ sptr;
417  }
418  sptr -= w; // go back to top-left
419  }
420  else if (h == 3) //////////////////////////////////////////////////
421  {
422  // topmost points ( [ (6^) 4 1 ] / 11 )^T
423  for (int i = 0; i < w; ++ i)
424  {
425  *dptr++ = sptr[ 0] * (6.0F / 11.0F) +
426  sptr[ w] * (4.0F / 11.0F) +
427  sptr[w2] * (1.0F / 11.0F);
428  ++ sptr;
429  }
430  sptr -= w; // go back to top-left
431 
432  // since factor >= 2 we don't need any other point.
433  }
434  else ///////////////////////////////// general case for height >= 4
435  {
436  int i1 = 0, i2 = 0;
437  int const skip = (factor - 1) * w;
438 
439  // topmost points ( [ (6^) 4 1 ] / 11 )^T
440  for (int i = 0; i < w; ++ i)
441  {
442  *dptr++ = sptr[ 0] * (6.0F / 11.0F) +
443  sptr[ w] * (4.0F / 11.0F) +
444  sptr[w2] * (1.0F / 11.0F);
445  ++ sptr;
446  }
447  sptr -= w; // go back to top-left
448  ++ i2;
449  i1 += factor;
450 
451  // second point skipped since factor >= 2
452  sptr += (skip - w);
453 
454  // rest of the column except last 2 points ( [ 1^ 4 (6) 4 1 ] / 16 )T
455  while((i1 < (h-2)) && (i2 < h2))
456  {
457  for (int i = 0; i < w; ++ i)
458  {
459  *dptr++ = (sptr[ 0] + sptr[w4]) * (1.0F / 16.0F) +
460  (sptr[ w] + sptr[w3]) * (4.0F / 16.0F) +
461  sptr[w2] * (6.0F / 16.0F);
462  ++ sptr;
463  }
464  sptr += skip;
465  i1 += factor;
466  ++ i2;
467  }
468 
469  // need special case for second to last point?
470  if ((i2 < h2) && (i1 == (h-2)))
471  {
472  sptr = source.const_end() - w4;
473  // before last points ( [ 1^ 4 (6) 4 ] / 15 )T
474  for (int i = 0; i < w; ++ i)
475  {
476  *dptr++ = sptr[ 0] * (1.0F / 15.0F) +
477  (sptr[ w] + sptr[w3]) * (4.0F / 15.0F) +
478  sptr[w2] * (6.0F / 15.0F);
479  ++ sptr;
480  }
481  i1 += factor;
482  ++i2;
483  }
484 
485  // need special case for last point?
486  if ((i2 < h2) && (i1 == (h-1)))
487  {
488  sptr = source.const_end() - w3;
489  // last points ( [ 1^ 4 (6) ] / 11 )T
490  for (int i = 0; i < w; ++ i)
491  {
492  *dptr++ = sptr[ 0] * (1.0F / 11.0F) +
493  sptr[ w] * (4.0F / 11.0F) +
494  sptr[w2] * (6.0F / 11.0F);
495  ++ sptr;
496  }
497  }
498  }
499  return result;
500 }
501 
502 
503 #endif // INCLUDE_NRT_IMAGEPROC_RESHAPING_DETAILS_RESCALEIMPL_H
504