[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

splineimageview.hxx VIGRA

1 /************************************************************************/
2 /* */
3 /* Copyright 1998-2004 by Ullrich Koethe */
4 /* */
5 /* This file is part of the VIGRA computer vision library. */
6 /* The VIGRA Website is */
7 /* http://hci.iwr.uni-heidelberg.de/vigra/ */
8 /* Please direct questions, bug reports, and contributions to */
9 /* ullrich.koethe@iwr.uni-heidelberg.de or */
10 /* vigra@informatik.uni-hamburg.de */
11 /* */
12 /* Permission is hereby granted, free of charge, to any person */
13 /* obtaining a copy of this software and associated documentation */
14 /* files (the "Software"), to deal in the Software without */
15 /* restriction, including without limitation the rights to use, */
16 /* copy, modify, merge, publish, distribute, sublicense, and/or */
17 /* sell copies of the Software, and to permit persons to whom the */
18 /* Software is furnished to do so, subject to the following */
19 /* conditions: */
20 /* */
21 /* The above copyright notice and this permission notice shall be */
22 /* included in all copies or substantial portions of the */
23 /* Software. */
24 /* */
25 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
26 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
27 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
28 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
29 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
30 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
31 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
32 /* OTHER DEALINGS IN THE SOFTWARE. */
33 /* */
34 /************************************************************************/
35 
36 #ifndef VIGRA_SPLINEIMAGEVIEW_HXX
37 #define VIGRA_SPLINEIMAGEVIEW_HXX
38 
39 #include "mathutil.hxx"
40 #include "recursiveconvolution.hxx"
41 #include "splines.hxx"
42 #include "array_vector.hxx"
43 #include "basicimage.hxx"
44 #include "copyimage.hxx"
45 #include "tinyvector.hxx"
46 #include "fixedpoint.hxx"
47 #include "multi_array.hxx"
48 
49 namespace vigra {
50 
51 /********************************************************/
52 /* */
53 /* SplineImageView */
54 /* */
55 /********************************************************/
56 /** \brief Create a continuous view onto a discrete image using splines.
57 
58  This class is very useful if image values or derivatives at arbitrary
59  real-valued coordinates are needed. Access at such coordinates is implemented by
60  interpolating the given discrete image values with a spline of the
61  specified <tt>ORDER</TT>. Continuous derivatives are available up to
62  degree <tt>ORDER-1</TT>. If the requested coordinates are near the image border,
63  reflective boundary conditions are applied. In principle, this class can also be used
64  for image resizing, but here the functions from the <tt>resize...</tt> family are
65  more efficient, since they exploit the regularity of the sampling grid.
66 
67  The <tt>SplineImageView</tt> template is explicitly specialized to make it as efficient as possible.
68  In particular, unnecessary copying of the image is avoided when the iterators passed
69  in the constructor originate from a \ref vigra::BasicImage. In addition, these specializations
70  provide function <tt>unchecked(...)</tt> that do not perform bounds checking. If the original image
71  is not a variant of \ref vigra::BasicImage, one can customize the internal representation by
72  using \ref vigra::SplineImageView0 or \ref vigra::SplineImageView1.
73 
74  <b>Usage:</b>
75 
76  <b>\#include</b> <vigra/splineimageview.hxx><br>
77  Namespace: vigra
78 
79  \code
80  BImage img(w,h);
81  ... // fill img
82 
83  // construct spline view for quadratic interpolation
84  SplineImageView<2, double> spi2(srcImageRange(img));
85 
86  double x = ..., y = ...;
87  double v2 = spi2(x, y);
88 
89  // construct spline view for linear interpolation
90  SplineImageView<1, UInt32> spi1(srcImageRange(img));
91 
92  UInt32 v1 = spi1(x, y);
93 
94  FixedPoint<16, 15> fx(...), fy(...);
95  UInt32 vf = spi1.unchecked(fx, fy); // caller is sure that (fx, fy) are valid coordinates
96  \endcode
97 */
98 template <int ORDER, class VALUETYPE>
99 class SplineImageView
100 {
101  typedef typename NumericTraits<VALUETYPE>::RealPromote InternalValue;
102 
103  public:
104 
105  /** The view's value type (return type of access and derivative functions).
106  */
107  typedef VALUETYPE value_type;
108 
109  /** The view's squared norm type (return type of g2() etc.).
110  */
111  typedef typename NormTraits<VALUETYPE>::SquaredNormType SquaredNormType;
112 
113  /** The view's size type.
114  */
115  typedef Size2D size_type;
116 
117  /** The view's difference type.
118  */
119  typedef TinyVector<double, 2> difference_type;
120 
121  /** The order of the spline used.
122  */
123  enum StaticOrder { order = ORDER };
124 
125  /** The type of the internal image holding the spline coefficients.
126  */
127  typedef BasicImage<InternalValue> InternalImage;
128 
129  private:
130  typedef typename InternalImage::traverser InternalTraverser;
131  typedef typename InternalTraverser::row_iterator InternalRowIterator;
132  typedef typename InternalTraverser::column_iterator InternalColumnIterator;
133  typedef BSpline<ORDER, double> Spline;
134 
135  enum { ksize_ = ORDER + 1, kcenter_ = ORDER / 2 };
136 
137  public:
138 
139  /** Construct SplineImageView for a 2D MultiArrayView.
140 
141  If <tt>skipPrefiltering = true</tt> (default: <tt>false</tt>), the recursive
142  prefilter of the cardinal spline function is not applied, resulting
143  in an approximating (smoothing) rather than interpolating spline. This is
144  especially useful if customized prefilters are to be applied.
145  */
146  template <class U, class S>
147  SplineImageView(MultiArrayView<2, U, S> const & s, bool skipPrefiltering = false)
148  : w_(s.shape(0)), h_(s.shape(1)), w1_(w_-1), h1_(h_-1),
149  x0_(kcenter_), x1_(w_ - kcenter_ - 2), y0_(kcenter_), y1_(h_ - kcenter_ - 2),
150  image_(w_, h_),
151  x_(-1.0), y_(-1.0),
152  u_(-1.0), v_(-1.0)
153  {
154  copyImage(srcImageRange(s), destImage(image_));
155  if(!skipPrefiltering)
156  init();
157  }
158 
159  /** Construct SplineImageView for an image given by \ref ImageIterators and \ref DataAccessors.
160 
161  If <tt>skipPrefiltering = true</tt> (default: <tt>false</tt>), the recursive
162  prefilter of the cardinal spline function is not applied, resulting
163  in an approximating (smoothing) rather than interpolating spline. This is
164  especially useful if customized prefilters are to be applied.
165  */
166  template <class SrcIterator, class SrcAccessor>
167  SplineImageView(SrcIterator is, SrcIterator iend, SrcAccessor sa, bool skipPrefiltering = false)
168  : w_(iend.x - is.x), h_(iend.y - is.y), w1_(w_-1), h1_(h_-1),
169  x0_(kcenter_), x1_(w_ - kcenter_ - 2), y0_(kcenter_), y1_(h_ - kcenter_ - 2),
170  image_(w_, h_),
171  x_(-1.0), y_(-1.0),
172  u_(-1.0), v_(-1.0)
173  {
174  copyImage(srcIterRange(is, iend, sa), destImage(image_));
175  if(!skipPrefiltering)
176  init();
177  }
178 
179  /** Construct SplineImageView for an image given by \ref ArgumentObjectFactories.
180 
181  If <tt>skipPrefiltering = true</tt> (default: <tt>false</tt>), the recursive
182  prefilter of the cardinal spline function is not applied, resulting
183  in an approximating (smoothing) rather than interpolating spline. This is
184  especially useful if customized prefilters are to be applied.
185  */
186  template <class SrcIterator, class SrcAccessor>
187  SplineImageView(triple<SrcIterator, SrcIterator, SrcAccessor> s, bool skipPrefiltering = false)
188  : w_(s.second.x - s.first.x), h_(s.second.y - s.first.y), w1_(w_-1), h1_(h_-1),
189  x0_(kcenter_), x1_(w_ - kcenter_ - 2), y0_(kcenter_), y1_(h_ - kcenter_ - 2),
190  image_(w_, h_),
191  x_(-1.0), y_(-1.0),
192  u_(-1.0), v_(-1.0)
193  {
194  copyImage(srcIterRange(s.first, s.second, s.third), destImage(image_));
195  if(!skipPrefiltering)
196  init();
197  }
198 
199  /** Access interpolated function at real-valued coordinate <tt>(x, y)</tt>.
200  If <tt>(x, y)</tt> is near the image border or outside the image, the value
201  is calculated with reflective boundary conditions. An exception is thrown if the
202  coordinate is outside the first reflection.
203  */
204  value_type operator()(double x, double y) const;
205 
206  /** Access derivative of order <tt>(dx, dy)</tt> at real-valued coordinate <tt>(x, y)</tt>.
207  If <tt>(x, y)</tt> is near the image border or outside the image, the value
208  is calculated with reflective boundary conditions. An exception is thrown if the
209  coordinate is outside the first reflection.
210  */
211  value_type operator()(double x, double y, unsigned int dx, unsigned int dy) const;
212 
213  /** Access 1st derivative in x-direction at real-valued coordinate <tt>(x, y)</tt>.
214  Equivalent to <tt>splineView(x, y, 1, 0)</tt>.
215  */
216  value_type dx(double x, double y) const
217  { return operator()(x, y, 1, 0); }
218 
219  /** Access 1st derivative in y-direction at real-valued coordinate <tt>(x, y)</tt>.
220  Equivalent to <tt>splineView(x, y, 0, 1)</tt>.
221  */
222  value_type dy(double x, double y) const
223  { return operator()(x, y, 0, 1); }
224 
225  /** Access 2nd derivative in x-direction at real-valued coordinate <tt>(x, y)</tt>.
226  Equivalent to <tt>splineView(x, y, 2, 0)</tt>.
227  */
228  value_type dxx(double x, double y) const
229  { return operator()(x, y, 2, 0); }
230 
231  /** Access mixed 2nd derivative at real-valued coordinate <tt>(x, y)</tt>.
232  Equivalent to <tt>splineView(x, y, 1, 1)</tt>.
233  */
234  value_type dxy(double x, double y) const
235  { return operator()(x, y, 1, 1); }
236 
237  /** Access 2nd derivative in y-direction at real-valued coordinate <tt>(x, y)</tt>.
238  Equivalent to <tt>splineView(x, y, 0, 2)</tt>.
239  */
240  value_type dyy(double x, double y) const
241  { return operator()(x, y, 0, 2); }
242 
243  /** Access 3rd derivative in x-direction at real-valued coordinate <tt>(x, y)</tt>.
244  Equivalent to <tt>splineView(x, y, 3, 0)</tt>.
245  */
246  value_type dx3(double x, double y) const
247  { return operator()(x, y, 3, 0); }
248 
249  /** Access 3rd derivative in y-direction at real-valued coordinate <tt>(x, y)</tt>.
250  Equivalent to <tt>splineView(x, y, 0, 3)</tt>.
251  */
252  value_type dy3(double x, double y) const
253  { return operator()(x, y, 0, 3); }
254 
255  /** Access mixed 3rd derivative dxxy at real-valued coordinate <tt>(x, y)</tt>.
256  Equivalent to <tt>splineView(x, y, 2, 1)</tt>.
257  */
258  value_type dxxy(double x, double y) const
259  { return operator()(x, y, 2, 1); }
260 
261  /** Access mixed 3rd derivative dxyy at real-valued coordinate <tt>(x, y)</tt>.
262  Equivalent to <tt>splineView(x, y, 1, 2)</tt>.
263  */
264  value_type dxyy(double x, double y) const
265  { return operator()(x, y, 1, 2); }
266 
267  /** Access interpolated function at real-valued coordinate <tt>d</tt>.
268  Equivalent to <tt>splineView(d[0], d[1])</tt>.
269  */
270  value_type operator()(difference_type const & d) const
271  { return operator()(d[0], d[1]); }
272 
273  /** Access derivative of order <tt>(dx, dy)</tt> at real-valued coordinate <tt>d</tt>.
274  Equivalent to <tt>splineView(d[0], d[1], dx, dy)</tt>.
275  */
276  value_type operator()(difference_type const & d, unsigned int dx, unsigned int dy) const
277  { return operator()(d[0], d[1], dx, dy); }
278 
279  /** Access 1st derivative in x-direction at real-valued coordinate <tt>d</tt>.
280  Equivalent to <tt>splineView.dx(d[0], d[1])</tt>.
281  */
282  value_type dx(difference_type const & d) const
283  { return dx(d[0], d[1]); }
284 
285  /** Access 1st derivative in y-direction at real-valued coordinate <tt>d</tt>.
286  Equivalent to <tt>splineView.dy(d[0], d[1])</tt>.
287  */
288  value_type dy(difference_type const & d) const
289  { return dy(d[0], d[1]); }
290 
291  /** Access 2nd derivative in x-direction at real-valued coordinate <tt>d</tt>.
292  Equivalent to <tt>splineView.dxx(d[0], d[1])</tt>.
293  */
294  value_type dxx(difference_type const & d) const
295  { return dxx(d[0], d[1]); }
296 
297  /** Access mixed 2nd derivative at real-valued coordinate <tt>d</tt>.
298  Equivalent to <tt>splineView.dxy(d[0], d[1])</tt>.
299  */
300  value_type dxy(difference_type const & d) const
301  { return dxy(d[0], d[1]); }
302 
303  /** Access 2nd derivative in y-direction at real-valued coordinate <tt>d</tt>.
304  Equivalent to <tt>splineView.dyy(d[0], d[1])</tt>.
305  */
306  value_type dyy(difference_type const & d) const
307  { return dyy(d[0], d[1]); }
308 
309  /** Access 3rd derivative in x-direction at real-valued coordinate <tt>d</tt>.
310  Equivalent to <tt>splineView.dx3(d[0], d[1])</tt>.
311  */
312  value_type dx3(difference_type const & d) const
313  { return dx3(d[0], d[1]); }
314 
315  /** Access 3rd derivative in y-direction at real-valued coordinate <tt>d</tt>.
316  Equivalent to <tt>splineView.dy3(d[0], d[1])</tt>.
317  */
318  value_type dy3(difference_type const & d) const
319  { return dy3(d[0], d[1]); }
320 
321  /** Access mixed 3rd derivative dxxy at real-valued coordinate <tt>d</tt>.
322  Equivalent to <tt>splineView.dxxy(d[0], d[1])</tt>.
323  */
324  value_type dxxy(difference_type const & d) const
325  { return dxxy(d[0], d[1]); }
326 
327  /** Access mixed 3rd derivative dxyy at real-valued coordinate <tt>d</tt>.
328  Equivalent to <tt>splineView.dxyy(d[0], d[1])</tt>.
329  */
330  value_type dxyy(difference_type const & d) const
331  { return dxyy(d[0], d[1]); }
332 
333  /** Access gradient squared magnitude at real-valued coordinate <tt>(x, y)</tt>.
334  */
335  SquaredNormType g2(double x, double y) const;
336 
337  /** Access 1st derivative in x-direction of gradient squared magnitude
338  at real-valued coordinate <tt>(x, y)</tt>.
339  */
340  SquaredNormType g2x(double x, double y) const;
341 
342  /** Access 1st derivative in y-direction of gradient squared magnitude
343  at real-valued coordinate <tt>(x, y)</tt>.
344  */
345  SquaredNormType g2y(double x, double y) const;
346 
347  /** Access 2nd derivative in x-direction of gradient squared magnitude
348  at real-valued coordinate <tt>(x, y)</tt>.
349  */
350  SquaredNormType g2xx(double x, double y) const;
351 
352  /** Access mixed 2nd derivative of gradient squared magnitude
353  at real-valued coordinate <tt>(x, y)</tt>.
354  */
355  SquaredNormType g2xy(double x, double y) const;
356 
357  /** Access 2nd derivative in y-direction of gradient squared magnitude
358  at real-valued coordinate <tt>(x, y)</tt>.
359  */
360  SquaredNormType g2yy(double x, double y) const;
361 
362  /** Access gradient squared magnitude at real-valued coordinate <tt>d</tt>.
363  */
364  SquaredNormType g2(difference_type const & d) const
365  { return g2(d[0], d[1]); }
366 
367  /** Access 1st derivative in x-direction of gradient squared magnitude
368  at real-valued coordinate <tt>d</tt>.
369  */
370  SquaredNormType g2x(difference_type const & d) const
371  { return g2x(d[0], d[1]); }
372 
373  /** Access 1st derivative in y-direction of gradient squared magnitude
374  at real-valued coordinate <tt>d</tt>.
375  */
376  SquaredNormType g2y(difference_type const & d) const
377  { return g2y(d[0], d[1]); }
378 
379  /** Access 2nd derivative in x-direction of gradient squared magnitude
380  at real-valued coordinate <tt>d</tt>.
381  */
382  SquaredNormType g2xx(difference_type const & d) const
383  { return g2xx(d[0], d[1]); }
384 
385  /** Access mixed 2nd derivative of gradient squared magnitude
386  at real-valued coordinate <tt>d</tt>.
387  */
388  SquaredNormType g2xy(difference_type const & d) const
389  { return g2xy(d[0], d[1]); }
390 
391  /** Access 2nd derivative in y-direction of gradient squared magnitude
392  at real-valued coordinate <tt>d</tt>.
393  */
394  SquaredNormType g2yy(difference_type const & d) const
395  { return g2yy(d[0], d[1]); }
396 
397  /** The width of the image.
398  <tt>0 <= x <= width()-1</tt> is required for all access functions.
399  */
400  unsigned int width() const
401  { return w_; }
402 
403  /** The height of the image.
404  <tt>0 <= y <= height()-1</tt> is required for all access functions.
405  */
406  unsigned int height() const
407  { return h_; }
408 
409  /** The size of the image.
410  <tt>0 <= x <= size().x-1</tt> and <tt>0 <= y <= size().y-1</tt>
411  are required for all access functions.
412  */
413  size_type size() const
414  { return size_type(w_, h_); }
415 
416  /** The shape of the image.
417  Same as size(), except for the return type.
418  */
419  TinyVector<unsigned int, 2> shape() const
420  { return TinyVector<unsigned int, 2>(w_, h_); }
421 
422  /** The internal image holding the spline coefficients.
423  */
424  InternalImage const & image() const
425  {
426  return image_;
427  }
428 
429  /** Get the array of polynomial coefficients for the facet containing
430  the point <tt>(x, y)</tt>. The array <tt>res</tt> must have
431  dimension <tt>(ORDER+1)x(ORDER+1)</tt>. From these coefficients, the
432  value of the interpolated function can be calculated by the following
433  algorithm
434 
435  \code
436  SplineImageView<ORDER, float> view(...);
437  double x = ..., y = ...;
438  double dx, dy;
439 
440  // calculate the local facet coordinates of x and y
441  if(ORDER % 2)
442  {
443  // odd order => facet coordinates between 0 and 1
444  dx = x - floor(x);
445  dy = y - floor(y);
446  }
447  else
448  {
449  // even order => facet coordinates between -0.5 and 0.5
450  dx = x - floor(x + 0.5);
451  dy = y - floor(y + 0.5);
452  }
453 
454  BasicImage<float> coefficients;
455  view.coefficientArray(x, y, coefficients);
456 
457  float f_x_y = 0.0;
458  for(int ny = 0; ny < ORDER + 1; ++ny)
459  for(int nx = 0; nx < ORDER + 1; ++nx)
460  f_x_y += pow(dx, nx) * pow(dy, ny) * coefficients(nx, ny);
461 
462  assert(abs(f_x_y - view(x, y)) < 1e-6);
463  \endcode
464  */
465  template <class Array>
466  void coefficientArray(double x, double y, Array & res) const;
467 
468  /** Check if x is in the original image range.
469  Equivalent to <tt>0 <= x <= width()-1</tt>.
470  */
471  bool isInsideX(double x) const
472  {
473  return x >= 0.0 && x <= width()-1.0;
474  }
475 
476  /** Check if y is in the original image range.
477  Equivalent to <tt>0 <= y <= height()-1</tt>.
478  */
479  bool isInsideY(double y) const
480  {
481  return y >= 0.0 && y <= height()-1.0;
482  }
483 
484  /** Check if x and y are in the original image range.
485  Equivalent to <tt>0 <= x <= width()-1</tt> and <tt>0 <= y <= height()-1</tt>.
486  */
487  bool isInside(double x, double y) const
488  {
489  return isInsideX(x) && isInsideY(y);
490  }
491 
492  /** Check if x and y are in the valid range. Points outside the original image range are computed
493  by reflective boundary conditions, but only within the first reflection.
494  Equivalent to <tt>-width() + ORDER/2 + 2 < x < 2*width() - ORDER/2 - 2</tt> and
495  <tt>-height() + ORDER/2 + 2 < y < 2*height() - ORDER/2 - 2</tt>.
496  */
497  bool isValid(double x, double y) const
498  {
499  return x < w1_ + x1_ && x > -x1_ && y < h1_ + y1_ && y > -y1_;
500  }
501 
502  /** Check whether the points <tt>(x0, y0)</tt> and <tt>(x1, y1)</tt> are in
503  the same spline facet. For odd order splines, facets span the range
504  <tt>(floor(x), floor(x)+1) x (floor(y), floor(y)+1)</tt> (i.e. we have
505  integer facet borders), whereas even order splines have facet between
506  half integer values
507  <tt>(floor(x)-0.5, floor(x)+0.5) x (floor(y)-0.5, floor(y)+0.5)</tt>.
508  */
509  bool sameFacet(double x0, double y0, double x1, double y1) const
510  {
511  x0 = VIGRA_CSTD::floor((ORDER % 2) ? x0 : x0 + 0.5);
512  y0 = VIGRA_CSTD::floor((ORDER % 2) ? y0 : y0 + 0.5);
513  x1 = VIGRA_CSTD::floor((ORDER % 2) ? x1 : x1 + 0.5);
514  y1 = VIGRA_CSTD::floor((ORDER % 2) ? y1 : y1 + 0.5);
515  return x0 == x1 && y0 == y1;
516  }
517 
518  protected:
519 
520  void init();
521  void calculateIndices(double x, double y) const;
522  void coefficients(double t, double * const & c) const;
523  void derivCoefficients(double t, unsigned int d, double * const & c) const;
524  value_type convolve() const;
525 
526  unsigned int w_, h_;
527  int w1_, h1_;
528  double x0_, x1_, y0_, y1_;
529  InternalImage image_;
530  Spline k_;
531  mutable double x_, y_, u_, v_, kx_[ksize_], ky_[ksize_];
532  mutable int ix_[ksize_], iy_[ksize_];
533 };
534 
535 template <int ORDER, class VALUETYPE>
536 void SplineImageView<ORDER, VALUETYPE>::init()
537 {
538  ArrayVector<double> const & b = k_.prefilterCoefficients();
539 
540  for(unsigned int i=0; i<b.size(); ++i)
541  {
542  recursiveFilterX(srcImageRange(image_), destImage(image_), b[i], BORDER_TREATMENT_REFLECT);
543  recursiveFilterY(srcImageRange(image_), destImage(image_), b[i], BORDER_TREATMENT_REFLECT);
544  }
545 }
546 
547 namespace detail
548 {
549 
550 template <int i>
551 struct SplineImageViewUnrollLoop1
552 {
553  template <class Array>
554  static void exec(int c0, Array c)
555  {
556  SplineImageViewUnrollLoop1<i-1>::exec(c0, c);
557  c[i] = c0 + i;
558  }
559 };
560 
561 template <>
562 struct SplineImageViewUnrollLoop1<0>
563 {
564  template <class Array>
565  static void exec(int c0, Array c)
566  {
567  c[0] = c0;
568  }
569 };
570 
571 template <int i, class ValueType>
572 struct SplineImageViewUnrollLoop2
573 {
574  template <class Array1, class RowIterator, class Array2>
575  static ValueType
576  exec(Array1 k, RowIterator r, Array2 x)
577  {
578  return ValueType(k[i] * r[x[i]]) + SplineImageViewUnrollLoop2<i-1, ValueType>::exec(k, r, x);
579  }
580 };
581 
582 template <class ValueType>
583 struct SplineImageViewUnrollLoop2<0, ValueType>
584 {
585  template <class Array1, class RowIterator, class Array2>
586  static ValueType
587  exec(Array1 k, RowIterator r, Array2 x)
588  {
589  return ValueType(k[0] * r[x[0]]);
590  }
591 };
592 
593 } // namespace detail
594 
595 template <int ORDER, class VALUETYPE>
596 void
597 SplineImageView<ORDER, VALUETYPE>::calculateIndices(double x, double y) const
598 {
599  if(x == x_ && y == y_)
600  return; // still in cache
601 
602  if(x > x0_ && x < x1_ && y > y0_ && y < y1_)
603  {
604  detail::SplineImageViewUnrollLoop1<ORDER>::exec(
605  (ORDER % 2) ? int(x - kcenter_) : int(x + 0.5 - kcenter_), ix_);
606  detail::SplineImageViewUnrollLoop1<ORDER>::exec(
607  (ORDER % 2) ? int(y - kcenter_) : int(y + 0.5 - kcenter_), iy_);
608 
609  u_ = x - ix_[kcenter_];
610  v_ = y - iy_[kcenter_];
611  }
612  else
613  {
614  vigra_precondition(isValid(x,y),
615  "SplineImageView::calculateIndices(): coordinates out of range.");
616 
617  int xCenter = (ORDER % 2) ?
618  (int)VIGRA_CSTD::floor(x) :
619  (int)VIGRA_CSTD::floor(x + 0.5);
620  int yCenter = (ORDER % 2) ?
621  (int)VIGRA_CSTD::floor(y) :
622  (int)VIGRA_CSTD::floor(y + 0.5);
623 
624  if(x >= x1_)
625  {
626  for(int i = 0; i < ksize_; ++i)
627  ix_[i] = w1_ - vigra::abs(w1_ - xCenter - (i - kcenter_));
628  }
629  else
630  {
631  for(int i = 0; i < ksize_; ++i)
632  ix_[i] = vigra::abs(xCenter - (kcenter_ - i));
633  }
634  if(y >= y1_)
635  {
636  for(int i = 0; i < ksize_; ++i)
637  iy_[i] = h1_ - vigra::abs(h1_ - yCenter - (i - kcenter_));
638  }
639  else
640  {
641  for(int i = 0; i < ksize_; ++i)
642  iy_[i] = vigra::abs(yCenter - (kcenter_ - i));
643  }
644  u_ = x - xCenter;
645  v_ = y - yCenter;
646  }
647  x_ = x;
648  y_ = y;
649 }
650 
651 template <int ORDER, class VALUETYPE>
652 void SplineImageView<ORDER, VALUETYPE>::coefficients(double t, double * const & c) const
653 {
654  t += kcenter_;
655  for(int i = 0; i<ksize_; ++i)
656  c[i] = k_(t-i);
657 }
658 
659 template <int ORDER, class VALUETYPE>
660 void SplineImageView<ORDER, VALUETYPE>::derivCoefficients(double t,
661  unsigned int d, double * const & c) const
662 {
663  t += kcenter_;
664  for(int i = 0; i<ksize_; ++i)
665  c[i] = k_(t-i, d);
666 }
667 
668 template <int ORDER, class VALUETYPE>
669 VALUETYPE SplineImageView<ORDER, VALUETYPE>::convolve() const
670 {
671  typedef typename NumericTraits<VALUETYPE>::RealPromote RealPromote;
672  RealPromote sum;
673  sum = RealPromote(
674  ky_[0]*detail::SplineImageViewUnrollLoop2<ORDER, RealPromote>::exec(kx_, image_.rowBegin(iy_[0]), ix_));
675 
676  for(int j=1; j<ksize_; ++j)
677  {
678  sum += RealPromote(
679  ky_[j]*detail::SplineImageViewUnrollLoop2<ORDER, RealPromote>::exec(kx_, image_.rowBegin(iy_[j]), ix_));
680  }
681  return detail::RequiresExplicitCast<VALUETYPE>::cast(sum);
682 }
683 
684 template <int ORDER, class VALUETYPE>
685 template <class Array>
686 void
687 SplineImageView<ORDER, VALUETYPE>::coefficientArray(double x, double y, Array & res) const
688 {
689  typedef typename Array::value_type ResType;
690  typename Spline::WeightMatrix const & weights = Spline::weights();
691  ResType tmp[ksize_][ksize_];
692 
693  calculateIndices(x, y);
694  for(int j=0; j<ksize_; ++j)
695  {
696  for(int i=0; i<ksize_; ++i)
697  {
698  tmp[i][j] = ResType();
699  for(int k=0; k<ksize_; ++k)
700  {
701  tmp[i][j] += weights[i][k]*image_(ix_[k], iy_[j]);
702  }
703  }
704  }
705  for(int j=0; j<ksize_; ++j)
706  {
707  for(int i=0; i<ksize_; ++i)
708  {
709  res(i,j) = ResType();
710  for(int k=0; k<ksize_; ++k)
711  {
712  res(i,j) += weights[j][k]*tmp[i][k];
713  }
714  }
715  }
716 }
717 
718 template <int ORDER, class VALUETYPE>
719 VALUETYPE SplineImageView<ORDER, VALUETYPE>::operator()(double x, double y) const
720 {
721  calculateIndices(x, y);
722  coefficients(u_, kx_);
723  coefficients(v_, ky_);
724  return convolve();
725 }
726 
727 template <int ORDER, class VALUETYPE>
728 VALUETYPE SplineImageView<ORDER, VALUETYPE>::operator()(double x, double y,
729  unsigned int dx, unsigned int dy) const
730 {
731  calculateIndices(x, y);
732  derivCoefficients(u_, dx, kx_);
733  derivCoefficients(v_, dy, ky_);
734  return convolve();
735 }
736 
737 template <int ORDER, class VALUETYPE>
738 typename SplineImageView<ORDER, VALUETYPE>::SquaredNormType
739 SplineImageView<ORDER, VALUETYPE>::g2(double x, double y) const
740 {
741  return squaredNorm(dx(x,y)) + squaredNorm(dy(x,y));
742 }
743 
744 template <int ORDER, class VALUETYPE>
745 typename SplineImageView<ORDER, VALUETYPE>::SquaredNormType
746 SplineImageView<ORDER, VALUETYPE>::g2x(double x, double y) const
747 {
748  return SquaredNormType(2.0)*(dot(dx(x,y), dxx(x,y)) + dot(dy(x,y), dxy(x,y)));
749 }
750 
751 template <int ORDER, class VALUETYPE>
752 typename SplineImageView<ORDER, VALUETYPE>::SquaredNormType
753 SplineImageView<ORDER, VALUETYPE>::g2y(double x, double y) const
754 {
755  return SquaredNormType(2.0)*(dot(dx(x,y), dxy(x,y)) + dot(dy(x,y), dyy(x,y)));
756 }
757 
758 template <int ORDER, class VALUETYPE>
759 typename SplineImageView<ORDER, VALUETYPE>::SquaredNormType
760 SplineImageView<ORDER, VALUETYPE>::g2xx(double x, double y) const
761 {
762  return SquaredNormType(2.0)*(squaredNorm(dxx(x,y)) + dot(dx(x,y), dx3(x,y)) +
763  squaredNorm(dxy(x,y)) + dot(dy(x,y), dxxy(x,y)));
764 }
765 
766 template <int ORDER, class VALUETYPE>
767 typename SplineImageView<ORDER, VALUETYPE>::SquaredNormType
768 SplineImageView<ORDER, VALUETYPE>::g2yy(double x, double y) const
769 {
770  return SquaredNormType(2.0)*(squaredNorm(dxy(x,y)) + dot(dx(x,y), dxyy(x,y)) +
771  squaredNorm(dyy(x,y)) + dot(dy(x,y), dy3(x,y)));
772 }
773 
774 template <int ORDER, class VALUETYPE>
775 typename SplineImageView<ORDER, VALUETYPE>::SquaredNormType
776 SplineImageView<ORDER, VALUETYPE>::g2xy(double x, double y) const
777 {
778  return SquaredNormType(2.0)*(dot(dx(x,y), dxxy(x,y)) + dot(dy(x,y), dxyy(x,y)) +
779  dot(dxy(x,y), dxx(x,y) + dyy(x,y)));
780 }
781 
782 /********************************************************/
783 /* */
784 /* SplineImageView0 */
785 /* */
786 /********************************************************/
787 template <class VALUETYPE, class INTERNAL_INDEXER>
788 class SplineImageView0Base
789 {
790  typedef typename INTERNAL_INDEXER::value_type InternalValue;
791  public:
792  typedef VALUETYPE value_type;
793  typedef typename NormTraits<VALUETYPE>::SquaredNormType SquaredNormType;
794  typedef Size2D size_type;
795  typedef TinyVector<double, 2> difference_type;
796  enum StaticOrder { order = 0 };
797 
798  public:
799 
800  SplineImageView0Base(unsigned int w, unsigned int h)
801  : w_(w), h_(h)
802  {}
803 
804  SplineImageView0Base(int w, int h, INTERNAL_INDEXER i)
805  : w_(w), h_(h), internalIndexer_(i)
806  {}
807 
808  template <unsigned IntBits1, unsigned FractionalBits1,
809  unsigned IntBits2, unsigned FractionalBits2>
810  value_type unchecked(FixedPoint<IntBits1, FractionalBits1> x,
811  FixedPoint<IntBits2, FractionalBits2> y) const
812  {
813  return internalIndexer_(round(x), round(y));
814  }
815 
816  template <unsigned IntBits1, unsigned FractionalBits1,
817  unsigned IntBits2, unsigned FractionalBits2>
818  value_type unchecked(FixedPoint<IntBits1, FractionalBits1> x,
819  FixedPoint<IntBits2, FractionalBits2> y,
820  unsigned int dx, unsigned int dy) const
821  {
822  if((dx != 0) || (dy != 0))
823  return NumericTraits<VALUETYPE>::zero();
824  return unchecked(x, y);
825  }
826 
827  value_type unchecked(double x, double y) const
828  {
829  return internalIndexer_((int)(x + 0.5), (int)(y + 0.5));
830  }
831 
832  value_type unchecked(double x, double y, unsigned int dx, unsigned int dy) const
833  {
834  if((dx != 0) || (dy != 0))
835  return NumericTraits<VALUETYPE>::zero();
836  return unchecked(x, y);
837  }
838 
839  value_type operator()(double x, double y) const
840  {
841  int ix, iy;
842  if(x < 0.0)
843  {
844  ix = (int)(-x + 0.5);
845  vigra_precondition(ix <= (int)w_ - 1,
846  "SplineImageView::operator(): coordinates out of range.");
847  }
848  else
849  {
850  ix = (int)(x + 0.5);
851  if(ix >= (int)w_)
852  {
853  ix = 2*w_-2-ix;
854  vigra_precondition(ix >= 0,
855  "SplineImageView::operator(): coordinates out of range.");
856  }
857  }
858  if(y < 0.0)
859  {
860  iy = (int)(-y + 0.5);
861  vigra_precondition(iy <= (int)h_ - 1,
862  "SplineImageView::operator(): coordinates out of range.");
863  }
864  else
865  {
866  iy = (int)(y + 0.5);
867  if(iy >= (int)h_)
868  {
869  iy = 2*h_-2-iy;
870  vigra_precondition(iy >= 0,
871  "SplineImageView::operator(): coordinates out of range.");
872  }
873  }
874  return internalIndexer_(ix, iy);
875  }
876 
877  value_type operator()(double x, double y, unsigned int dx, unsigned int dy) const
878  {
879  if((dx != 0) || (dy != 0))
880  return NumericTraits<VALUETYPE>::zero();
881  return operator()(x, y);
882  }
883 
884  value_type dx(double /*x*/, double /*y*/) const
885  { return NumericTraits<VALUETYPE>::zero(); }
886 
887  value_type dy(double /*x*/, double /*y*/) const
888  { return NumericTraits<VALUETYPE>::zero(); }
889 
890  value_type dxx(double /*x*/, double /*y*/) const
891  { return NumericTraits<VALUETYPE>::zero(); }
892 
893  value_type dxy(double /*x*/, double /*y*/) const
894  { return NumericTraits<VALUETYPE>::zero(); }
895 
896  value_type dyy(double /*x*/, double /*y*/) const
897  { return NumericTraits<VALUETYPE>::zero(); }
898 
899  value_type dx3(double /*x*/, double /*y*/) const
900  { return NumericTraits<VALUETYPE>::zero(); }
901 
902  value_type dy3(double /*x*/, double /*y*/) const
903  { return NumericTraits<VALUETYPE>::zero(); }
904 
905  value_type dxxy(double /*x*/, double /*y*/) const
906  { return NumericTraits<VALUETYPE>::zero(); }
907 
908  value_type dxyy(double /*x*/, double /*y*/) const
909  { return NumericTraits<VALUETYPE>::zero(); }
910 
911  value_type operator()(difference_type const & d) const
912  { return operator()(d[0], d[1]); }
913 
914  value_type operator()(difference_type const & d, unsigned int dx, unsigned int dy) const
915  { return operator()(d[0], d[1], dx, dy); }
916 
917  value_type dx(difference_type const & /*d*/) const
918  { return NumericTraits<VALUETYPE>::zero(); }
919 
920  value_type dy(difference_type const & /*d*/) const
921  { return NumericTraits<VALUETYPE>::zero(); }
922 
923  value_type dxx(difference_type const & /*d*/) const
924  { return NumericTraits<VALUETYPE>::zero(); }
925 
926  value_type dxy(difference_type const & /*d*/) const
927  { return NumericTraits<VALUETYPE>::zero(); }
928 
929  value_type dyy(difference_type const & /*d*/) const
930  { return NumericTraits<VALUETYPE>::zero(); }
931 
932  value_type dx3(difference_type const & /*d*/) const
933  { return NumericTraits<VALUETYPE>::zero(); }
934 
935  value_type dy3(difference_type const & /*d*/) const
936  { return NumericTraits<VALUETYPE>::zero(); }
937 
938  value_type dxxy(difference_type const & /*d*/) const
939  { return NumericTraits<VALUETYPE>::zero(); }
940 
941  value_type dxyy(difference_type const & /*d*/) const
942  { return NumericTraits<VALUETYPE>::zero(); }
943 
944  SquaredNormType g2(double /*x*/, double /*y*/) const
945  { return NumericTraits<SquaredNormType>::zero(); }
946 
947  SquaredNormType g2x(double /*x*/, double /*y*/) const
948  { return NumericTraits<SquaredNormType>::zero(); }
949 
950  SquaredNormType g2y(double /*x*/, double /*y*/) const
951  { return NumericTraits<SquaredNormType>::zero(); }
952 
953  SquaredNormType g2xx(double /*x*/, double /*y*/) const
954  { return NumericTraits<SquaredNormType>::zero(); }
955 
956  SquaredNormType g2xy(double /*x*/, double /*y*/) const
957  { return NumericTraits<SquaredNormType>::zero(); }
958 
959  SquaredNormType g2yy(double /*x*/, double /*y*/) const
960  { return NumericTraits<SquaredNormType>::zero(); }
961 
962  SquaredNormType g2(difference_type const & /*d*/) const
963  { return NumericTraits<SquaredNormType>::zero(); }
964 
965  SquaredNormType g2x(difference_type const & /*d*/) const
966  { return NumericTraits<SquaredNormType>::zero(); }
967 
968  SquaredNormType g2y(difference_type const & /*d*/) const
969  { return NumericTraits<SquaredNormType>::zero(); }
970 
971  SquaredNormType g2xx(difference_type const & /*d*/) const
972  { return NumericTraits<SquaredNormType>::zero(); }
973 
974  SquaredNormType g2xy(difference_type const & /*d*/) const
975  { return NumericTraits<SquaredNormType>::zero(); }
976 
977  SquaredNormType g2yy(difference_type const & /*d*/) const
978  { return NumericTraits<SquaredNormType>::zero(); }
979 
980  unsigned int width() const
981  { return w_; }
982 
983  unsigned int height() const
984  { return h_; }
985 
986  size_type size() const
987  { return size_type(w_, h_); }
988 
989  TinyVector<unsigned int, 2> shape() const
990  { return TinyVector<unsigned int, 2>(w_, h_); }
991 
992  template <class Array>
993  void coefficientArray(double x, double y, Array & res) const
994  {
995  res(0, 0) = operator()(x,y);
996  }
997 
998  bool isInsideX(double x) const
999  {
1000  return x >= 0.0 && x <= width() - 1.0;
1001  }
1002 
1003  bool isInsideY(double y) const
1004  {
1005  return y >= 0.0 && y <= height() - 1.0;
1006  }
1007 
1008  bool isInside(double x, double y) const
1009  {
1010  return isInsideX(x) && isInsideY(y);
1011  }
1012 
1013  bool isValid(double x, double y) const
1014  {
1015  return x < 2.0*w_-2.0 && x > 1.0-w_ && y < 2.0*h_-2.0 && y > 1.0-h_;
1016  }
1017 
1018  bool sameFacet(double x0, double y0, double x1, double y1) const
1019  {
1020  x0 = VIGRA_CSTD::floor(x0 + 0.5);
1021  y0 = VIGRA_CSTD::floor(y0 + 0.5);
1022  x1 = VIGRA_CSTD::floor(x1 + 0.5);
1023  y1 = VIGRA_CSTD::floor(y1 + 0.5);
1024  return x0 == x1 && y0 == y1;
1025  }
1026 
1027  protected:
1028  unsigned int w_, h_;
1029  INTERNAL_INDEXER internalIndexer_;
1030 };
1031 
1032 /** \brief Create an image view for nearest-neighbor interpolation.
1033 
1034  This class behaves like \ref vigra::SplineImageView&lt;0, ...&gt;, but one can pass
1035  an additional template argument that determined the internal representation of the image.
1036  If this is equal to the argument type passed in the constructor, the image is not copied.
1037  By default, this works for \ref vigra::BasicImage, \ref vigra::BasicImageView,
1038  \ref vigra::MultiArray&lt;2, ...&gt;, and \ref vigra::MultiArrayView&lt;2, ...&gt;.
1039 
1040 */
1041 template <class VALUETYPE, class INTERNAL_TRAVERSER = typename BasicImage<VALUETYPE>::const_traverser>
1042 class SplineImageView0
1043 : public SplineImageView0Base<VALUETYPE, INTERNAL_TRAVERSER>
1044 {
1045  typedef SplineImageView0Base<VALUETYPE, INTERNAL_TRAVERSER> Base;
1046  public:
1047  typedef typename Base::value_type value_type;
1048  typedef typename Base::SquaredNormType SquaredNormType;
1049  typedef typename Base::size_type size_type;
1050  typedef typename Base::difference_type difference_type;
1051  enum StaticOrder { order = Base::order };
1052  typedef BasicImage<VALUETYPE> InternalImage;
1053 
1054  protected:
1055  typedef typename IteratorTraits<INTERNAL_TRAVERSER>::mutable_iterator InternalTraverser;
1056  typedef typename IteratorTraits<InternalTraverser>::DefaultAccessor InternalAccessor;
1057  typedef typename IteratorTraits<INTERNAL_TRAVERSER>::const_iterator InternalConstTraverser;
1058  typedef typename IteratorTraits<InternalConstTraverser>::DefaultAccessor InternalConstAccessor;
1059 
1060  public:
1061 
1062  /* when traverser and accessor types passed to the constructor are the same as the corresponding
1063  internal types, we need not copy the image (speed up)
1064  */
1065  SplineImageView0(InternalTraverser is, InternalTraverser iend, InternalAccessor /*sa*/)
1066  : Base(iend.x - is.x, iend.y - is.y, is)
1067  {}
1068 
1069  SplineImageView0(triple<InternalTraverser, InternalTraverser, InternalAccessor> s)
1070  : Base(s.second.x - s.first.x, s.second.y - s.first.y, s.first)
1071  {}
1072 
1073  SplineImageView0(InternalConstTraverser is, InternalConstTraverser iend, InternalConstAccessor /*sa*/)
1074  : Base(iend.x - is.x, iend.y - is.y, is)
1075  {}
1076 
1077  SplineImageView0(triple<InternalConstTraverser, InternalConstTraverser, InternalConstAccessor> s)
1078  : Base(s.second.x - s.first.x, s.second.y - s.first.y, s.first)
1079  {}
1080 
1081  template<class T, class SU>
1082  SplineImageView0(MultiArrayView<2, T, SU> const & i)
1083  : Base(i.shape(0), i.shape(1)),
1084  image_(i.shape(0), i.shape(1))
1085  {
1086  for(unsigned int y=0; y<this->height(); ++y)
1087  for(unsigned int x=0; x<this->width(); ++x)
1088  image_(x,y) = detail::RequiresExplicitCast<VALUETYPE>::cast(i(x,y));
1089  this->internalIndexer_ = image_.upperLeft();
1090  }
1091 
1092  template <class SrcIterator, class SrcAccessor>
1093  SplineImageView0(SrcIterator is, SrcIterator iend, SrcAccessor sa)
1094  : Base(iend.x - is.x, iend.y - is.y),
1095  image_(iend - is)
1096  {
1097  copyImage(srcIterRange(is, iend, sa), destImage(image_));
1098  this->internalIndexer_ = image_.upperLeft();
1099  }
1100 
1101  template <class SrcIterator, class SrcAccessor>
1102  SplineImageView0(triple<SrcIterator, SrcIterator, SrcAccessor> s)
1103  : Base(s.second.x - s.first.x, s.second.y - s.first.y),
1104  image_(s.second - s.first)
1105  {
1106  copyImage(s, destImage(image_));
1107  this->internalIndexer_ = image_.upperLeft();
1108  }
1109 
1110  InternalImage const & image() const
1111  { return image_; }
1112 
1113  protected:
1114  InternalImage image_;
1115 };
1116 
1117 template <class VALUETYPE, class StridedOrUnstrided>
1118 class SplineImageView0<VALUETYPE, MultiArrayView<2, VALUETYPE, StridedOrUnstrided> >
1119 : public SplineImageView0Base<VALUETYPE, MultiArrayView<2, VALUETYPE, StridedOrUnstrided> >
1120 {
1121  typedef SplineImageView0Base<VALUETYPE, MultiArrayView<2, VALUETYPE, StridedOrUnstrided> > Base;
1122  public:
1123  typedef typename Base::value_type value_type;
1124  typedef typename Base::SquaredNormType SquaredNormType;
1125  typedef typename Base::size_type size_type;
1126  typedef typename Base::difference_type difference_type;
1127  enum StaticOrder { order = Base::order };
1128  typedef BasicImage<VALUETYPE> InternalImage;
1129 
1130  protected:
1131  typedef MultiArrayView<2, VALUETYPE, StridedOrUnstrided> InternalIndexer;
1132 
1133  public:
1134 
1135  /* when traverser and accessor types passed to the constructor are the same as the corresponding
1136  internal types, we need not copy the image (speed up)
1137  */
1138  SplineImageView0(InternalIndexer const & i)
1139  : Base(i.shape(0), i.shape(1), i)
1140  {}
1141 
1142  template<class T, class SU>
1143  SplineImageView0(MultiArrayView<2, T, SU> const & i)
1144  : Base(i.shape(0), i.shape(1)),
1145  image_(i.shape(0), i.shape(1))
1146  {
1147  for(unsigned int y=0; y<this->height(); ++y)
1148  for(unsigned int x=0; x<this->width(); ++x)
1149  image_(x,y) = detail::RequiresExplicitCast<VALUETYPE>::cast(i(x,y));
1150  this->internalIndexer_ = InternalIndexer(typename InternalIndexer::difference_type(this->width(), this->height()),
1151  image_.data());
1152  }
1153 
1154  template <class SrcIterator, class SrcAccessor>
1155  SplineImageView0(SrcIterator is, SrcIterator iend, SrcAccessor sa)
1156  : Base(iend.x - is.x, iend.y - is.y),
1157  image_(iend-is)
1158  {
1159  copyImage(srcIterRange(is, iend, sa), destImage(image_));
1160  this->internalIndexer_ = InternalIndexer(typename InternalIndexer::difference_type(this->width(), this->height()),
1161  image_.data());
1162  }
1163 
1164  template <class SrcIterator, class SrcAccessor>
1165  SplineImageView0(triple<SrcIterator, SrcIterator, SrcAccessor> s)
1166  : Base(s.second.x - s.first.x, s.second.y - s.first.y),
1167  image_(s.second - s.first)
1168  {
1169  copyImage(s, destImage(image_));
1170  this->internalIndexer_ = InternalIndexer(typename InternalIndexer::difference_type(this->width(), this->height()),
1171  image_.data());
1172  }
1173 
1174  InternalImage const & image() const
1175  { return image_; }
1176 
1177  protected:
1178  InternalImage image_;
1179 };
1180 
1181 template <class VALUETYPE>
1182 class SplineImageView<0, VALUETYPE>
1183 : public SplineImageView0<VALUETYPE>
1184 {
1185  typedef SplineImageView0<VALUETYPE> Base;
1186  public:
1187  typedef typename Base::value_type value_type;
1188  typedef typename Base::SquaredNormType SquaredNormType;
1189  typedef typename Base::size_type size_type;
1190  typedef typename Base::difference_type difference_type;
1191  enum StaticOrder { order = Base::order };
1192  typedef typename Base::InternalImage InternalImage;
1193 
1194  protected:
1195  typedef typename Base::InternalTraverser InternalTraverser;
1196  typedef typename Base::InternalAccessor InternalAccessor;
1197  typedef typename Base::InternalConstTraverser InternalConstTraverser;
1198  typedef typename Base::InternalConstAccessor InternalConstAccessor;
1199 
1200 public:
1201 
1202  /* when traverser and accessor types passed to the constructor are the same as the corresponding
1203  internal types, we need not copy the image (speed up)
1204  */
1205  SplineImageView(InternalTraverser is, InternalTraverser iend, InternalAccessor sa, bool /* unused */ = false)
1206  : Base(is, iend, sa)
1207  {}
1208 
1209  SplineImageView(triple<InternalTraverser, InternalTraverser, InternalAccessor> s, bool /* unused */ = false)
1210  : Base(s)
1211  {}
1212 
1213  SplineImageView(InternalConstTraverser is, InternalConstTraverser iend, InternalConstAccessor sa, bool /* unused */ = false)
1214  : Base(is, iend, sa)
1215  {}
1216 
1217  SplineImageView(triple<InternalConstTraverser, InternalConstTraverser, InternalConstAccessor> s, bool /* unused */ = false)
1218  : Base(s)
1219  {}
1220 
1221  template <class SrcIterator, class SrcAccessor>
1222  SplineImageView(SrcIterator is, SrcIterator iend, SrcAccessor sa, bool /* unused */ = false)
1223  : Base(is, iend, sa)
1224  {
1225  copyImage(srcIterRange(is, iend, sa), destImage(this->image_));
1226  }
1227 
1228  template <class SrcIterator, class SrcAccessor>
1229  SplineImageView(triple<SrcIterator, SrcIterator, SrcAccessor> s, bool /* unused */ = false)
1230  : Base(s)
1231  {
1232  copyImage(s, destImage(this->image_));
1233  }
1234 };
1235 
1236 /********************************************************/
1237 /* */
1238 /* SplineImageView1 */
1239 /* */
1240 /********************************************************/
1241 template <class VALUETYPE, class INTERNAL_INDEXER>
1242 class SplineImageView1Base
1243 {
1244  typedef typename INTERNAL_INDEXER::value_type InternalValue;
1245  public:
1246  typedef VALUETYPE value_type;
1247  typedef Size2D size_type;
1248  typedef typename NormTraits<VALUETYPE>::SquaredNormType SquaredNormType;
1249  typedef TinyVector<double, 2> difference_type;
1250  enum StaticOrder { order = 1 };
1251 
1252  public:
1253 
1254  SplineImageView1Base(unsigned int w, unsigned int h)
1255  : w_(w), h_(h)
1256  {}
1257 
1258  SplineImageView1Base(int w, int h, INTERNAL_INDEXER i)
1259  : w_(w), h_(h), internalIndexer_(i)
1260  {}
1261 
1262  template <unsigned IntBits1, unsigned FractionalBits1,
1263  unsigned IntBits2, unsigned FractionalBits2>
1264  value_type unchecked(FixedPoint<IntBits1, FractionalBits1> x,
1265  FixedPoint<IntBits2, FractionalBits2> y) const
1266  {
1267  int ix = floor(x);
1268  FixedPoint<0, FractionalBits1> tx = frac(x - FixedPoint<IntBits1, FractionalBits1>(ix));
1269  FixedPoint<0, FractionalBits1> dtx = dual_frac(tx);
1270  if(ix == (int)w_ - 1)
1271  {
1272  --ix;
1273  tx.value = FixedPoint<0, FractionalBits1>::ONE;
1274  dtx.value = 0;
1275  }
1276  int iy = floor(y);
1277  FixedPoint<0, FractionalBits2> ty = frac(y - FixedPoint<IntBits2, FractionalBits2>(iy));
1278  FixedPoint<0, FractionalBits2> dty = dual_frac(ty);
1279  if(iy == (int)h_ - 1)
1280  {
1281  --iy;
1282  ty.value = FixedPoint<0, FractionalBits2>::ONE;
1283  dty.value = 0;
1284  }
1285  return fixed_point_cast<value_type>(
1286  dty*(dtx*fixedPoint(internalIndexer_(ix,iy)) +
1287  tx*fixedPoint(internalIndexer_(ix+1,iy))) +
1288  ty *(dtx*fixedPoint(internalIndexer_(ix,iy+1)) +
1289  tx*fixedPoint(internalIndexer_(ix+1,iy+1))));
1290  }
1291 
1292  template <unsigned IntBits1, unsigned FractionalBits1,
1293  unsigned IntBits2, unsigned FractionalBits2>
1294  value_type unchecked(FixedPoint<IntBits1, FractionalBits1> x,
1295  FixedPoint<IntBits2, FractionalBits2> y,
1296  unsigned int dx, unsigned int dy) const
1297  {
1298  int ix = floor(x);
1299  FixedPoint<0, FractionalBits1> tx = frac(x - FixedPoint<IntBits1, FractionalBits1>(ix));
1300  FixedPoint<0, FractionalBits1> dtx = dual_frac(tx);
1301  if(ix == (int)w_ - 1)
1302  {
1303  --ix;
1304  tx.value = FixedPoint<0, FractionalBits1>::ONE;
1305  dtx.value = 0;
1306  }
1307  int iy = floor(y);
1308  FixedPoint<0, FractionalBits2> ty = frac(y - FixedPoint<IntBits2, FractionalBits2>(iy));
1309  FixedPoint<0, FractionalBits2> dty = dual_frac(ty);
1310  if(iy == (int)h_ - 1)
1311  {
1312  --iy;
1313  ty.value = FixedPoint<0, FractionalBits2>::ONE;
1314  dty.value = 0;
1315  }
1316  switch(dx)
1317  {
1318  case 0:
1319  switch(dy)
1320  {
1321  case 0:
1322  return fixed_point_cast<value_type>(
1323  dty*(dtx*fixedPoint(internalIndexer_(ix,iy)) +
1324  tx*fixedPoint(internalIndexer_(ix+1,iy))) +
1325  ty *(dtx*fixedPoint(internalIndexer_(ix,iy+1)) +
1326  tx*fixedPoint(internalIndexer_(ix+1,iy+1))));
1327  case 1:
1328  return fixed_point_cast<value_type>(
1329  (dtx*fixedPoint(internalIndexer_(ix,iy+1)) + tx*fixedPoint(internalIndexer_(ix+1,iy+1))) -
1330  (dtx*fixedPoint(internalIndexer_(ix,iy)) + tx*fixedPoint(internalIndexer_(ix+1,iy))));
1331  default:
1332  return NumericTraits<VALUETYPE>::zero();
1333  }
1334  case 1:
1335  switch(dy)
1336  {
1337  case 0:
1338  return fixed_point_cast<value_type>(
1339  dty*(fixedPoint(internalIndexer_(ix+1,iy)) - fixedPoint(internalIndexer_(ix,iy))) +
1340  ty *(fixedPoint(internalIndexer_(ix+1,iy+1)) - fixedPoint(internalIndexer_(ix,iy+1))));
1341  case 1:
1342  return detail::RequiresExplicitCast<value_type>::cast(
1343  (internalIndexer_(ix+1,iy+1) - internalIndexer_(ix,iy+1)) -
1344  (internalIndexer_(ix+1,iy) - internalIndexer_(ix,iy)));
1345  default:
1346  return NumericTraits<VALUETYPE>::zero();
1347  }
1348  default:
1349  return NumericTraits<VALUETYPE>::zero();
1350  }
1351  }
1352 
1353  value_type unchecked(double x, double y) const
1354  {
1355  int ix = (int)std::floor(x);
1356  if(ix == (int)w_ - 1)
1357  --ix;
1358  double tx = x - ix;
1359  int iy = (int)std::floor(y);
1360  if(iy == (int)h_ - 1)
1361  --iy;
1362  double ty = y - iy;
1363  return NumericTraits<value_type>::fromRealPromote(
1364  (1.0-ty)*((1.0-tx)*internalIndexer_(ix,iy) + tx*internalIndexer_(ix+1,iy)) +
1365  ty *((1.0-tx)*internalIndexer_(ix,iy+1) + tx*internalIndexer_(ix+1,iy+1)));
1366  }
1367 
1368  value_type unchecked(double x, double y, unsigned int dx, unsigned int dy) const
1369  {
1370  int ix = (int)std::floor(x);
1371  if(ix == (int)w_ - 1)
1372  --ix;
1373  double tx = x - ix;
1374  int iy = (int)std::floor(y);
1375  if(iy == (int)h_ - 1)
1376  --iy;
1377  double ty = y - iy;
1378  switch(dx)
1379  {
1380  case 0:
1381  switch(dy)
1382  {
1383  case 0:
1384  return detail::RequiresExplicitCast<value_type>::cast(
1385  (1.0-ty)*((1.0-tx)*internalIndexer_(ix,iy) + tx*internalIndexer_(ix+1,iy)) +
1386  ty *((1.0-tx)*internalIndexer_(ix,iy+1) + tx*internalIndexer_(ix+1,iy+1)));
1387  case 1:
1388  return detail::RequiresExplicitCast<value_type>::cast(
1389  ((1.0-tx)*internalIndexer_(ix,iy+1) + tx*internalIndexer_(ix+1,iy+1)) -
1390  ((1.0-tx)*internalIndexer_(ix,iy) + tx*internalIndexer_(ix+1,iy)));
1391  default:
1392  return NumericTraits<VALUETYPE>::zero();
1393  }
1394  case 1:
1395  switch(dy)
1396  {
1397  case 0:
1398  return detail::RequiresExplicitCast<value_type>::cast(
1399  (1.0-ty)*(internalIndexer_(ix+1,iy) - internalIndexer_(ix,iy)) +
1400  ty *(internalIndexer_(ix+1,iy+1) - internalIndexer_(ix,iy+1)));
1401  case 1:
1402  return detail::RequiresExplicitCast<value_type>::cast(
1403  (internalIndexer_(ix+1,iy+1) - internalIndexer_(ix,iy+1)) -
1404  (internalIndexer_(ix+1,iy) - internalIndexer_(ix,iy)));
1405  default:
1406  return NumericTraits<VALUETYPE>::zero();
1407  }
1408  default:
1409  return NumericTraits<VALUETYPE>::zero();
1410  }
1411  }
1412 
1413  value_type operator()(double x, double y) const
1414  {
1415  return operator()(x, y, 0, 0);
1416  }
1417 
1418  value_type operator()(double x, double y, unsigned int dx, unsigned int dy) const
1419  {
1420  value_type mul = NumericTraits<value_type>::one();
1421  if(x < 0.0)
1422  {
1423  x = -x;
1424  vigra_precondition(x <= w_ - 1.0,
1425  "SplineImageView::operator(): coordinates out of range.");
1426  if(dx % 2)
1427  mul = -mul;
1428  }
1429  else if(x > w_ - 1.0)
1430  {
1431  x = 2.0*w_-2.0-x;
1432  vigra_precondition(x >= 0.0,
1433  "SplineImageView::operator(): coordinates out of range.");
1434  if(dx % 2)
1435  mul = -mul;
1436  }
1437  if(y < 0.0)
1438  {
1439  y = -y;
1440  vigra_precondition(y <= h_ - 1.0,
1441  "SplineImageView::operator(): coordinates out of range.");
1442  if(dy % 2)
1443  mul = -mul;
1444  }
1445  else if(y > h_ - 1.0)
1446  {
1447  y = 2.0*h_-2.0-y;
1448  vigra_precondition(y >= 0.0,
1449  "SplineImageView::operator(): coordinates out of range.");
1450  if(dy % 2)
1451  mul = -mul;
1452  }
1453  return mul*unchecked(x, y, dx, dy);
1454  }
1455 
1456  value_type dx(double x, double y) const
1457  { return operator()(x, y, 1, 0); }
1458 
1459  value_type dy(double x, double y) const
1460  { return operator()(x, y, 0, 1); }
1461 
1462  value_type dxx(double /*x*/, double /*y*/) const
1463  { return NumericTraits<VALUETYPE>::zero(); }
1464 
1465  value_type dxy(double x, double y) const
1466  { return operator()(x, y, 1, 1); }
1467 
1468  value_type dyy(double /*x*/, double /*y*/) const
1469  { return NumericTraits<VALUETYPE>::zero(); }
1470 
1471  value_type dx3(double /*x*/, double /*y*/) const
1472  { return NumericTraits<VALUETYPE>::zero(); }
1473 
1474  value_type dy3(double /*x*/, double /*y*/) const
1475  { return NumericTraits<VALUETYPE>::zero(); }
1476 
1477  value_type dxxy(double /*x*/, double /*y*/) const
1478  { return NumericTraits<VALUETYPE>::zero(); }
1479 
1480  value_type dxyy(double /*x*/, double /*y*/) const
1481  { return NumericTraits<VALUETYPE>::zero(); }
1482 
1483  value_type operator()(difference_type const & d) const
1484  { return operator()(d[0], d[1]); }
1485 
1486  value_type operator()(difference_type const & d, unsigned int dx, unsigned int dy) const
1487  { return operator()(d[0], d[1], dx, dy); }
1488 
1489  value_type dx(difference_type const & d) const
1490  { return operator()(d[0], d[1], 1, 0); }
1491 
1492  value_type dy(difference_type const & d) const
1493  { return operator()(d[0], d[1], 0, 1); }
1494 
1495  value_type dxx(difference_type const & /*d*/) const
1496  { return NumericTraits<VALUETYPE>::zero(); }
1497 
1498  value_type dxy(difference_type const & d) const
1499  { return operator()(d[0], d[1], 1, 1); }
1500 
1501  value_type dyy(difference_type const & /*d*/) const
1502  { return NumericTraits<VALUETYPE>::zero(); }
1503 
1504  value_type dx3(difference_type const & /*d*/) const
1505  { return NumericTraits<VALUETYPE>::zero(); }
1506 
1507  value_type dy3(difference_type const & /*d*/) const
1508  { return NumericTraits<VALUETYPE>::zero(); }
1509 
1510  value_type dxxy(difference_type const & /*d*/) const
1511  { return NumericTraits<VALUETYPE>::zero(); }
1512 
1513  value_type dxyy(difference_type const & /*d*/) const
1514  { return NumericTraits<VALUETYPE>::zero(); }
1515 
1516  SquaredNormType g2(double x, double y) const
1517  { return squaredNorm(dx(x,y)) + squaredNorm(dy(x,y)); }
1518 
1519  SquaredNormType g2x(double /*x*/, double /*y*/) const
1520  { return NumericTraits<SquaredNormType>::zero(); }
1521 
1522  SquaredNormType g2y(double /*x*/, double /*y*/) const
1523  { return NumericTraits<SquaredNormType>::zero(); }
1524 
1525  SquaredNormType g2xx(double /*x*/, double /*y*/) const
1526  { return NumericTraits<SquaredNormType>::zero(); }
1527 
1528  SquaredNormType g2xy(double /*x*/, double /*y*/) const
1529  { return NumericTraits<SquaredNormType>::zero(); }
1530 
1531  SquaredNormType g2yy(double /*x*/, double /*y*/) const
1532  { return NumericTraits<SquaredNormType>::zero(); }
1533 
1534  SquaredNormType g2(difference_type const & d) const
1535  { return g2(d[0], d[1]); }
1536 
1537  SquaredNormType g2x(difference_type const & /*d*/) const
1538  { return NumericTraits<SquaredNormType>::zero(); }
1539 
1540  SquaredNormType g2y(difference_type const & /*d*/) const
1541  { return NumericTraits<SquaredNormType>::zero(); }
1542 
1543  SquaredNormType g2xx(difference_type const & /*d*/) const
1544  { return NumericTraits<SquaredNormType>::zero(); }
1545 
1546  SquaredNormType g2xy(difference_type const & /*d*/) const
1547  { return NumericTraits<SquaredNormType>::zero(); }
1548 
1549  SquaredNormType g2yy(difference_type const & /*d*/) const
1550  { return NumericTraits<SquaredNormType>::zero(); }
1551 
1552  unsigned int width() const
1553  { return w_; }
1554 
1555  unsigned int height() const
1556  { return h_; }
1557 
1558  size_type size() const
1559  { return size_type(w_, h_); }
1560 
1561  TinyVector<unsigned int, 2> shape() const
1562  { return TinyVector<unsigned int, 2>(w_, h_); }
1563 
1564  template <class Array>
1565  void coefficientArray(double x, double y, Array & res) const;
1566 
1567  void calculateIndices(double x, double y, int & ix, int & iy, int & ix1, int & iy1) const;
1568 
1569  bool isInsideX(double x) const
1570  {
1571  return x >= 0.0 && x <= width() - 1.0;
1572  }
1573 
1574  bool isInsideY(double y) const
1575  {
1576  return y >= 0.0 && y <= height() - 1.0;
1577  }
1578 
1579  bool isInside(double x, double y) const
1580  {
1581  return isInsideX(x) && isInsideY(y);
1582  }
1583 
1584  bool isValid(double x, double y) const
1585  {
1586  return x < 2.0*w_-2.0 && x > 1.0-w_ && y < 2.0*h_-2.0 && y > 1.0-h_;
1587  }
1588 
1589  bool sameFacet(double x0, double y0, double x1, double y1) const
1590  {
1591  x0 = VIGRA_CSTD::floor(x0);
1592  y0 = VIGRA_CSTD::floor(y0);
1593  x1 = VIGRA_CSTD::floor(x1);
1594  y1 = VIGRA_CSTD::floor(y1);
1595  return x0 == x1 && y0 == y1;
1596  }
1597 
1598  protected:
1599  unsigned int w_, h_;
1600  INTERNAL_INDEXER internalIndexer_;
1601 };
1602 
1603 template <class VALUETYPE, class INTERNAL_INDEXER>
1604 template <class Array>
1605 void SplineImageView1Base<VALUETYPE, INTERNAL_INDEXER>::coefficientArray(double x, double y, Array & res) const
1606 {
1607  int ix, iy, ix1, iy1;
1608  calculateIndices(x, y, ix, iy, ix1, iy1);
1609  res(0,0) = internalIndexer_(ix,iy);
1610  res(1,0) = internalIndexer_(ix1,iy) - internalIndexer_(ix,iy);
1611  res(0,1) = internalIndexer_(ix,iy1) - internalIndexer_(ix,iy);
1612  res(1,1) = internalIndexer_(ix,iy) - internalIndexer_(ix1,iy) -
1613  internalIndexer_(ix,iy1) + internalIndexer_(ix1,iy1);
1614 }
1615 
1616 template <class VALUETYPE, class INTERNAL_INDEXER>
1617 void SplineImageView1Base<VALUETYPE, INTERNAL_INDEXER>::calculateIndices(double x, double y, int & ix, int & iy, int & ix1, int & iy1) const
1618 {
1619  if(x < 0.0)
1620  {
1621  x = -x;
1622  vigra_precondition(x <= w_ - 1.0,
1623  "SplineImageView::calculateIndices(): coordinates out of range.");
1624  ix = (int)VIGRA_CSTD::ceil(x);
1625  ix1 = ix - 1;
1626  }
1627  else if(x >= w_ - 1.0)
1628  {
1629  x = 2.0*w_-2.0-x;
1630  vigra_precondition(x > 0.0,
1631  "SplineImageView::calculateIndices(): coordinates out of range.");
1632  ix = (int)VIGRA_CSTD::ceil(x);
1633  ix1 = ix - 1;
1634  }
1635  else
1636  {
1637  ix = (int)VIGRA_CSTD::floor(x);
1638  ix1 = ix + 1;
1639  }
1640  if(y < 0.0)
1641  {
1642  y = -y;
1643  vigra_precondition(y <= h_ - 1.0,
1644  "SplineImageView::calculateIndices(): coordinates out of range.");
1645  iy = (int)VIGRA_CSTD::ceil(y);
1646  iy1 = iy - 1;
1647  }
1648  else if(y >= h_ - 1.0)
1649  {
1650  y = 2.0*h_-2.0-y;
1651  vigra_precondition(y > 0.0,
1652  "SplineImageView::calculateIndices(): coordinates out of range.");
1653  iy = (int)VIGRA_CSTD::ceil(y);
1654  iy1 = iy - 1;
1655  }
1656  else
1657  {
1658  iy = (int)VIGRA_CSTD::floor(y);
1659  iy1 = iy + 1;
1660  }
1661 }
1662 
1663 /** \brief Create an image view for bi-linear interpolation.
1664 
1665  This class behaves like \ref vigra::SplineImageView&lt;1, ...&gt;, but one can pass
1666  an additional template argument that determined the internal representation of the image.
1667  If this is equal to the argument type passed in the constructor, the image is not copied.
1668  By default, this works for \ref vigra::BasicImage, \ref vigra::BasicImageView,
1669  \ref vigra::MultiArray&lt;2, ...&gt;, and \ref vigra::MultiArrayView&lt;2, ...&gt;.
1670 
1671  In addition to the function provided by \ref vigra::SplineImageView, there are functions
1672  <tt>unchecked(x,y)</tt> and <tt>unchecked(x,y, xorder, yorder)</tt> which improve speed by
1673  not applying bounds checking and reflective border treatment (<tt>isInside(x, y)</tt> must
1674  be <tt>true</tt>), but otherwise behave identically to their checked counterparts.
1675  In addition, <tt>x</tt> and <tt>y</tt> can have type \ref vigra::FixedPoint instead of
1676  <tt>double</tt>.
1677 */
1678 template <class VALUETYPE, class INTERNAL_TRAVERSER = typename BasicImage<VALUETYPE>::const_traverser>
1679 class SplineImageView1
1680 : public SplineImageView1Base<VALUETYPE, INTERNAL_TRAVERSER>
1681 {
1682  typedef SplineImageView1Base<VALUETYPE, INTERNAL_TRAVERSER> Base;
1683  public:
1684  typedef typename Base::value_type value_type;
1685  typedef typename Base::SquaredNormType SquaredNormType;
1686  typedef typename Base::size_type size_type;
1687  typedef typename Base::difference_type difference_type;
1688  enum StaticOrder { order = Base::order };
1689  typedef BasicImage<VALUETYPE> InternalImage;
1690 
1691  protected:
1692  typedef typename IteratorTraits<INTERNAL_TRAVERSER>::mutable_iterator InternalTraverser;
1693  typedef typename IteratorTraits<InternalTraverser>::DefaultAccessor InternalAccessor;
1694  typedef typename IteratorTraits<INTERNAL_TRAVERSER>::const_iterator InternalConstTraverser;
1695  typedef typename IteratorTraits<InternalConstTraverser>::DefaultAccessor InternalConstAccessor;
1696 
1697  public:
1698 
1699  /* when traverser and accessor types passed to the constructor are the same as the corresponding
1700  internal types, we need not copy the image (speed up)
1701  */
1702  SplineImageView1(InternalTraverser is, InternalTraverser iend, InternalAccessor /*sa*/)
1703  : Base(iend.x - is.x, iend.y - is.y, is)
1704  {}
1705 
1706  SplineImageView1(triple<InternalTraverser, InternalTraverser, InternalAccessor> s)
1707  : Base(s.second.x - s.first.x, s.second.y - s.first.y, s.first)
1708  {}
1709 
1710  SplineImageView1(InternalConstTraverser is, InternalConstTraverser iend, InternalConstAccessor /*sa*/)
1711  : Base(iend.x - is.x, iend.y - is.y, is)
1712  {}
1713 
1714  SplineImageView1(triple<InternalConstTraverser, InternalConstTraverser, InternalConstAccessor> s)
1715  : Base(s.second.x - s.first.x, s.second.y - s.first.y, s.first)
1716  {}
1717 
1718  template<class T, class SU>
1719  SplineImageView1(MultiArrayView<2, T, SU> const & i)
1720  : Base(i.shape(0), i.shape(1)),
1721  image_(i.shape(0), i.shape(1))
1722  {
1723  for(unsigned int y=0; y<this->height(); ++y)
1724  for(unsigned int x=0; x<this->width(); ++x)
1725  image_(x,y) = detail::RequiresExplicitCast<VALUETYPE>::cast(i(x,y));
1726  this->internalIndexer_ = image_.upperLeft();
1727  }
1728 
1729  template <class SrcIterator, class SrcAccessor>
1730  SplineImageView1(SrcIterator is, SrcIterator iend, SrcAccessor sa)
1731  : Base(iend.x - is.x, iend.y - is.y),
1732  image_(iend - is)
1733  {
1734  copyImage(srcIterRange(is, iend, sa), destImage(image_));
1735  this->internalIndexer_ = image_.upperLeft();
1736  }
1737 
1738  template <class SrcIterator, class SrcAccessor>
1739  SplineImageView1(triple<SrcIterator, SrcIterator, SrcAccessor> s)
1740  : Base(s.second.x - s.first.x, s.second.y - s.first.y),
1741  image_(s.second - s.first)
1742  {
1743  copyImage(s, destImage(image_));
1744  this->internalIndexer_ = image_.upperLeft();
1745  }
1746 
1747  InternalImage const & image() const
1748  { return image_; }
1749 
1750  protected:
1751  InternalImage image_;
1752 };
1753 
1754 template <class VALUETYPE, class StridedOrUnstrided>
1755 class SplineImageView1<VALUETYPE, MultiArrayView<2, VALUETYPE, StridedOrUnstrided> >
1756 : public SplineImageView1Base<VALUETYPE, MultiArrayView<2, VALUETYPE, StridedOrUnstrided> >
1757 {
1758  typedef SplineImageView1Base<VALUETYPE, MultiArrayView<2, VALUETYPE, StridedOrUnstrided> > Base;
1759  public:
1760  typedef typename Base::value_type value_type;
1761  typedef typename Base::SquaredNormType SquaredNormType;
1762  typedef typename Base::size_type size_type;
1763  typedef typename Base::difference_type difference_type;
1764  enum StaticOrder { order = Base::order };
1765  typedef BasicImage<VALUETYPE> InternalImage;
1766 
1767  protected:
1768  typedef MultiArrayView<2, VALUETYPE, StridedOrUnstrided> InternalIndexer;
1769 
1770  public:
1771 
1772  /* when traverser and accessor types passed to the constructor are the same as the corresponding
1773  internal types, we need not copy the image (speed up)
1774  */
1775  SplineImageView1(InternalIndexer const & i)
1776  : Base(i.shape(0), i.shape(1), i)
1777  {}
1778 
1779  template<class T, class SU>
1780  SplineImageView1(MultiArrayView<2, T, SU> const & i)
1781  : Base(i.shape(0), i.shape(1)),
1782  image_(i.shape(0), i.shape(1))
1783  {
1784  copyImage(srcImageRange(i), destImage(image_));
1785  // for(unsigned int y=0; y<this->height(); ++y)
1786  // for(unsigned int x=0; x<this->width(); ++x)
1787  // image_(x,y) = detail::RequiresExplicitCast<VALUETYPE>::cast(i(x,y));
1788  this->internalIndexer_ = InternalIndexer(typename InternalIndexer::difference_type(this->width(), this->height()),
1789  image_.data());
1790  }
1791 
1792  template <class SrcIterator, class SrcAccessor>
1793  SplineImageView1(SrcIterator is, SrcIterator iend, SrcAccessor sa)
1794  : Base(iend.x - is.x, iend.y - is.y),
1795  image_(iend-is)
1796  {
1797  copyImage(srcIterRange(is, iend, sa), destImage(image_));
1798  this->internalIndexer_ = InternalIndexer(typename InternalIndexer::difference_type(this->width(), this->height()),
1799  image_.data());
1800  }
1801 
1802  template <class SrcIterator, class SrcAccessor>
1803  SplineImageView1(triple<SrcIterator, SrcIterator, SrcAccessor> s)
1804  : Base(s.second.x - s.first.x, s.second.y - s.first.y),
1805  image_(s.second - s.first)
1806  {
1807  copyImage(s, destImage(image_));
1808  this->internalIndexer_ = InternalIndexer(typename InternalIndexer::difference_type(this->width(), this->height()),
1809  image_.data());
1810  }
1811 
1812  InternalImage const & image() const
1813  { return image_; }
1814 
1815  protected:
1816  InternalImage image_;
1817 };
1818 
1819 template <class VALUETYPE>
1820 class SplineImageView<1, VALUETYPE>
1821 : public SplineImageView1<VALUETYPE>
1822 {
1823  typedef SplineImageView1<VALUETYPE> Base;
1824  public:
1825  typedef typename Base::value_type value_type;
1826  typedef typename Base::SquaredNormType SquaredNormType;
1827  typedef typename Base::size_type size_type;
1828  typedef typename Base::difference_type difference_type;
1829  enum StaticOrder { order = Base::order };
1830  typedef typename Base::InternalImage InternalImage;
1831 
1832  protected:
1833  typedef typename Base::InternalTraverser InternalTraverser;
1834  typedef typename Base::InternalAccessor InternalAccessor;
1835  typedef typename Base::InternalConstTraverser InternalConstTraverser;
1836  typedef typename Base::InternalConstAccessor InternalConstAccessor;
1837 
1838 public:
1839 
1840  /* when traverser and accessor types passed to the constructor are the same as the corresponding
1841  internal types, we need not copy the image (speed up)
1842  */
1843  SplineImageView(InternalTraverser is, InternalTraverser iend, InternalAccessor sa, bool /* unused */ = false)
1844  : Base(is, iend, sa)
1845  {}
1846 
1847  SplineImageView(triple<InternalTraverser, InternalTraverser, InternalAccessor> s, bool /* unused */ = false)
1848  : Base(s)
1849  {}
1850 
1851  SplineImageView(InternalConstTraverser is, InternalConstTraverser iend, InternalConstAccessor sa, bool /* unused */ = false)
1852  : Base(is, iend, sa)
1853  {}
1854 
1855  SplineImageView(triple<InternalConstTraverser, InternalConstTraverser, InternalConstAccessor> s, bool /* unused */ = false)
1856  : Base(s)
1857  {}
1858 
1859  template <class SrcIterator, class SrcAccessor>
1860  SplineImageView(SrcIterator is, SrcIterator iend, SrcAccessor sa, bool /* unused */ = false)
1861  : Base(is, iend, sa)
1862  {
1863  copyImage(srcIterRange(is, iend, sa), destImage(this->image_));
1864  }
1865 
1866  template <class SrcIterator, class SrcAccessor>
1867  SplineImageView(triple<SrcIterator, SrcIterator, SrcAccessor> s, bool /* unused */ = false)
1868  : Base(s)
1869  {
1870  copyImage(s, destImage(this->image_));
1871  }
1872 
1873  template<class T, class SU>
1874  SplineImageView(MultiArrayView<2, T, SU> const & i, bool /* unused */ = false)
1875  : Base(i)
1876  {}
1877 };
1878 
1879 } // namespace vigra
1880 
1881 
1882 #endif /* VIGRA_SPLINEIMAGEVIEW_HXX */
vigra::SplineImageView::g2
SquaredNormType g2(difference_type const &d) const
Definition: splineimageview.hxx:364
vigra::SplineImageView::g2
SquaredNormType g2(double x, double y) const
Definition: splineimageview.hxx:739
vigra::SplineImageView::isInside
bool isInside(double x, double y) const
Definition: splineimageview.hxx:487
vigra::SplineImageView::value_type
VALUETYPE value_type
Definition: splineimageview.hxx:107
vigra::SplineImageView::isInsideY
bool isInsideY(double y) const
Definition: splineimageview.hxx:479
vigra::dot
PromoteTraits< V1, V2 >::Promote dot(RGBValue< V1, RIDX1, GIDX1, BIDX1 > const &r1, RGBValue< V2, RIDX2, GIDX2, BIDX2 > const &r2)
dot product
Definition: rgbvalue.hxx:906
vigra::IteratorTraits
Export associated information for each image iterator.
Definition: iteratortraits.hxx:109
vigra::SplineImageView::width
unsigned int width() const
Definition: splineimageview.hxx:400
vigra::SplineImageView::operator()
value_type operator()(difference_type const &d) const
Definition: splineimageview.hxx:270
vigra::SplineImageView::g2y
SquaredNormType g2y(double x, double y) const
Definition: splineimageview.hxx:753
vigra::recursiveFilterY
void recursiveFilterY(...)
Performs 1 dimensional recursive filtering (1st and 2nd order) in y direction.
vigra::SplineImageView::g2x
SquaredNormType g2x(difference_type const &d) const
Definition: splineimageview.hxx:370
vigra::MultiArrayView::shape
const difference_type & shape() const
Definition: multi_array.hxx:1648
vigra::SplineImageView::coefficientArray
void coefficientArray(double x, double y, Array &res) const
Definition: splineimageview.hxx:687
vigra::SplineImageView::SplineImageView
SplineImageView(MultiArrayView< 2, U, S > const &s, bool skipPrefiltering=false)
Definition: splineimageview.hxx:147
vigra::SplineImageView::dyy
value_type dyy(difference_type const &d) const
Definition: splineimageview.hxx:306
vigra::SplineImageView::dxx
value_type dxx(difference_type const &d) const
Definition: splineimageview.hxx:294
vigra::squaredNorm
FFTWComplex< R >::SquaredNormType squaredNorm(const FFTWComplex< R > &a)
squared norm (= squared magnitude)
Definition: fftw3.hxx:1044
vigra::SplineImageView::g2y
SquaredNormType g2y(difference_type const &d) const
Definition: splineimageview.hxx:376
vigra::round
int round(FixedPoint< IntBits, FracBits > v)
rounding to the nearest integer.
Definition: fixedpoint.hxx:683
vigra::fixed_point_cast
TARGET fixed_point_cast(FixedPoint< IntBits, FracBits > v)
Definition: fixedpoint.hxx:486
vigra::SplineImageView::dyy
value_type dyy(double x, double y) const
Definition: splineimageview.hxx:240
vigra::SplineImageView::dxy
value_type dxy(difference_type const &d) const
Definition: splineimageview.hxx:300
vigra::SplineImageView::g2x
SquaredNormType g2x(double x, double y) const
Definition: splineimageview.hxx:746
vigra::SplineImageView::g2xy
SquaredNormType g2xy(difference_type const &d) const
Definition: splineimageview.hxx:388
vigra::SplineImageView::StaticOrder
StaticOrder
Definition: splineimageview.hxx:123
vigra::SplineImageView::operator()
value_type operator()(double x, double y) const
Definition: splineimageview.hxx:719
vigra::SplineImageView::g2xx
SquaredNormType g2xx(double x, double y) const
Definition: splineimageview.hxx:760
vigra::SplineImageView::dy3
value_type dy3(difference_type const &d) const
Definition: splineimageview.hxx:318
vigra::SplineImageView::isValid
bool isValid(double x, double y) const
Definition: splineimageview.hxx:497
vigra::dual_frac
FixedPoint< 0, FracBits > dual_frac(FixedPoint< IntBits, FracBits > v)
dual fractional part: 1 - frac(v).
Definition: fixedpoint.hxx:658
vigra::SplineImageView::dx
value_type dx(difference_type const &d) const
Definition: splineimageview.hxx:282
vigra::Size2D
Two dimensional size object.
Definition: diff2d.hxx:482
vigra::frac
FixedPoint< 0, FracBits > frac(FixedPoint< IntBits, FracBits > v)
fractional part.
Definition: fixedpoint.hxx:650
vigra::SplineImageView::g2yy
SquaredNormType g2yy(double x, double y) const
Definition: splineimageview.hxx:768
vigra::SplineImageView::g2xx
SquaredNormType g2xx(difference_type const &d) const
Definition: splineimageview.hxx:382
vigra::SplineImageView::dxxy
value_type dxxy(double x, double y) const
Definition: splineimageview.hxx:258
vigra::BasicImageIterator
Definition: basicimage.hxx:262
vigra::ArrayVector
Definition: array_vector.hxx:58
vigra::SplineImageView::dxx
value_type dxx(double x, double y) const
Definition: splineimageview.hxx:228
vigra::SplineImageView::image
InternalImage const & image() const
Definition: splineimageview.hxx:424
vigra::mul
void mul(FixedPoint< IntBits1, FracBits1 > l, FixedPoint< IntBits2, FracBits2 > r, FixedPoint< IntBits3, FracBits3 > &result)
multiplication with enforced result type.
Definition: fixedpoint.hxx:605
vigra::SplineImageView::operator()
value_type operator()(difference_type const &d, unsigned int dx, unsigned int dy) const
Definition: splineimageview.hxx:276
vigra::recursiveFilterX
void recursiveFilterX(...)
Performs 1 dimensional recursive filtering (1st and 2nd order) in x direction.
vigra::sum
NumericTraits< V >::Promote sum(TinyVectorBase< V, SIZE, D1, D2 > const &l)
sum of the vector's elements
Definition: tinyvector.hxx:2073
vigra::SplineImageView::SquaredNormType
NormTraits< VALUETYPE >::SquaredNormType SquaredNormType
Definition: splineimageview.hxx:111
vigra::SplineImageView::difference_type
TinyVector< double, 2 > difference_type
Definition: splineimageview.hxx:119
vigra::SplineImageView1
Create an image view for bi-linear interpolation.
Definition: splineimageview.hxx:1679
vigra::SplineImageView::size
size_type size() const
Definition: splineimageview.hxx:413
vigra::BasicImage::data
const_pointer data() const
Definition: basicimage.hxx:1059
vigra::SplineImageView::shape
TinyVector< unsigned int, 2 > shape() const
Definition: splineimageview.hxx:419
vigra::copyImage
void copyImage(...)
Copy source image into destination image.
vigra::SplineImageView::dxxy
value_type dxxy(difference_type const &d) const
Definition: splineimageview.hxx:324
vigra::SplineImageView::dx3
value_type dx3(double x, double y) const
Definition: splineimageview.hxx:246
vigra::SplineImageView::dx
value_type dx(double x, double y) const
Definition: splineimageview.hxx:216
vigra::SplineImageView::dy3
value_type dy3(double x, double y) const
Definition: splineimageview.hxx:252
vigra::SplineImageView::dx3
value_type dx3(difference_type const &d) const
Definition: splineimageview.hxx:312
vigra::SplineImageView::dxyy
value_type dxyy(difference_type const &d) const
Definition: splineimageview.hxx:330
vigra::SplineImageView::size_type
Size2D size_type
Definition: splineimageview.hxx:115
vigra::SplineImageView::g2yy
SquaredNormType g2yy(difference_type const &d) const
Definition: splineimageview.hxx:394
vigra::SplineImageView::SplineImageView
SplineImageView(triple< SrcIterator, SrcIterator, SrcAccessor > s, bool skipPrefiltering=false)
Definition: splineimageview.hxx:187
vigra::SplineImageView::sameFacet
bool sameFacet(double x0, double y0, double x1, double y1) const
Definition: splineimageview.hxx:509
vigra::abs
FFTWComplex< R >::NormType abs(const FFTWComplex< R > &a)
absolute value (= magnitude)
Definition: fftw3.hxx:1002
vigra::MultiArrayView
Base class for, and view to, vigra::MultiArray.
Definition: multi_array.hxx:704
vigra::SplineImageView::dy
value_type dy(difference_type const &d) const
Definition: splineimageview.hxx:288
vigra::SplineImageView::dxyy
value_type dxyy(double x, double y) const
Definition: splineimageview.hxx:264
vigra::SplineImageView::dy
value_type dy(double x, double y) const
Definition: splineimageview.hxx:222
vigra::SplineImageView::InternalImage
BasicImage< InternalValue > InternalImage
Definition: splineimageview.hxx:127
vigra::SplineImageView::g2xy
SquaredNormType g2xy(double x, double y) const
Definition: splineimageview.hxx:776
vigra::ceil
int ceil(FixedPoint< IntBits, FracBits > v)
rounding up.
Definition: fixedpoint.hxx:675
vigra::StandardAccessor
Encapsulate access to the values an iterator points to.
Definition: accessor.hxx:133
vigra::IteratorAdaptor
Quickly create 1-dimensional iterator adapters.
Definition: iteratoradapter.hxx:147
vigra::SplineImageView
Create a continuous view onto a discrete image using splines.
Definition: splineimageview.hxx:99
vigra::floor
int floor(FixedPoint< IntBits, FracBits > v)
rounding down.
Definition: fixedpoint.hxx:667
vigra::SplineImageView0
Create an image view for nearest-neighbor interpolation.
Definition: splineimageview.hxx:1042
vigra::SplineImageView::dxy
value_type dxy(double x, double y) const
Definition: splineimageview.hxx:234
vigra::SplineImageView::height
unsigned int height() const
Definition: splineimageview.hxx:406
vigra::SplineImageView::isInsideX
bool isInsideX(double x) const
Definition: splineimageview.hxx:471
vigra::BasicImage::upperLeft
traverser upperLeft()
Definition: basicimage.hxx:925
vigra::SplineImageView::SplineImageView
SplineImageView(SrcIterator is, SrcIterator iend, SrcAccessor sa, bool skipPrefiltering=false)
Definition: splineimageview.hxx:167

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.11.1 (Fri May 19 2017)