SLIP  1.4
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
bspline_interpolation.hpp
Go to the documentation of this file.
1 /*
2  * Copyright(c):
3  * Signal Image and Communications (SIC) Department
4  * http://www.sic.sp2mi.univ-poitiers.fr/
5  * - University of Poitiers, France http://www.univ-poitiers.fr
6  * - XLIM Institute UMR CNRS 7252 http://www.xlim.fr/
7  *
8  * and
9  *
10  * D2 Fluid, Thermic and Combustion
11  * - University of Poitiers, France http://www.univ-poitiers.fr
12  * - PPRIME Institute - UPR CNRS 3346 http://www.pprime.fr
13  * - ISAE-ENSMA http://www.ensma.fr
14  *
15  * Contributor(s):
16  * The SLIP team,
17  * Benoit Tremblais <tremblais_AT_sic.univ-poitiers.fr>,
18  * Laurent David <laurent.david_AT_lea.univ-poitiers.fr>,
19  * Ludovic Chatellier <ludovic.chatellier_AT_univ-poitiers.fr>,
20  * Lionel Thomas <lionel.thomas_AT_univ-poitiers.fr>,
21  * Denis Arrivault <arrivault_AT_sic.univ-poitiers.fr>,
22  * Julien Dombre <julien.dombre_AT_univ-poitiers.fr>.
23  *
24  * Description:
25  * The Simple Library of Image Processing (SLIP) is a new image processing
26  * library. It is written in the C++ language following as much as possible
27  * the ISO/ANSI C++ standard. It is consequently compatible with any system
28  * satisfying the ANSI C++ complience. It works on different Unix , Linux ,
29  * Mircrosoft Windows and Mac OS X plateforms. SLIP is a research library that
30  * was created by the Signal, Image and Communications (SIC) departement of
31  * the XLIM, UMR 7252 CNRS Institute in collaboration with the Fluids, Thermic
32  * and Combustion departement of the P', UPR 3346 CNRS Institute of the
33  * University of Poitiers.
34  *
35  * The SLIP Library source code has been registered to the APP (French Agency
36  * for the Protection of Programs) by the University of Poitiers and CNRS,
37  * under registration number IDDN.FR.001.300034.000.S.P.2010.000.21000.
38 
39  * http://www.sic.sp2mi.univ-poitiers.fr/slip/
40  *
41  * This software is governed by the CeCILL-C license under French law and
42  * abiding by the rules of distribution of free software. You can use,
43  * modify and/ or redistribute the software under the terms of the CeCILL-C
44  * license as circulated by CEA, CNRS and INRIA at the following URL
45  * http://www.cecill.info.
46  * As a counterpart to the access to the source code and rights to copy,
47  * modify and redistribute granted by the license, users are provided only
48  * with a limited warranty and the software's author, the holder of the
49  * economic rights, and the successive licensors have only limited
50  * liability.
51  *
52  * In this respect, the user's attention is drawn to the risks associated
53  * with loading, using, modifying and/or developing or reproducing the
54  * software by the user in light of its specific status of free software,
55  * that may mean that it is complicated to manipulate, and that also
56  * therefore means that it is reserved for developers and experienced
57  * professionals having in-depth computer knowledge. Users are therefore
58  * encouraged to load and test the software's suitability as regards their
59  * requirements in conditions enabling the security of their systems and/or
60  * data to be ensured and, more generally, to use and operate it in the
61  * same conditions as regards security.
62  *
63  * The fact that you are presently reading this means that you have had
64  * knowledge of the CeCILL-C license and that you accept its terms.
65  */
66 
67 
74 #ifndef SLIP_BSPLINE_INTERPOLATION_HPP
75 #define SLIP_BSPLINE_INTERPOLATION_HPP
76 
77 #include <cmath>
78 #include <limits>
79 #include "DPoint2d.hpp"
80 #include "DPoint3d.hpp"
81 #include "Array.hpp"
82 #include "Array2d.hpp"
83 #include "Array3d.hpp"
84 #include "macros.hpp"
85 #include "arithmetic_op.hpp"
86 #include "dynamic.hpp"
87 namespace slipalgo
88 {
90  /* @{ */
108 template <class Array>
109 inline
110 void BSplinePoles(const std::size_t degree,
111  Array &Poles)
112 {
113 
114  /* recover the poles from a lookup table */
115  switch (degree)
116  {
117  case 2L:
118  Poles.resize(1);
119  Poles[0] = std::sqrt(8.0) - 3.0;
120  break;
121  case 3L:
122  Poles.resize(1);
123  Poles[0] = std::sqrt(3.0) - 2.0;
124  break;
125  case 4L:
126  Poles.resize(2);
127  Poles[0] = std::sqrt(664.0 - std::sqrt(438976.0)) + std::sqrt(304.0) - 19.0;
128  Poles[1] = std::sqrt(664.0 + std::sqrt(438976.0)) - std::sqrt(304.0) - 19.0;
129  break;
130  case 5L:
131  Poles.resize(2);
132  Poles[0] = std::sqrt(135.0 / 2.0 - std::sqrt(17745.0 / 4.0)) + std::sqrt(105.0 / 4.0)
133  - 13.0 / 2.0;
134  Poles[1] = std::sqrt(135.0 / 2.0 + std::sqrt(17745.0 / 4.0)) - std::sqrt(105.0 / 4.0)
135  - 13.0 / 2.0;
136  break;
137  case 6L:
138  Poles.resize(3);
139  Poles[0] = -0.48829458930304475513011803888378906211227916123938;
140  Poles[1] = -0.081679271076237512597937765737059080653379610398148;
141  Poles[2] = -0.0014141518083258177510872439765585925278641690553467;
142  break;
143  case 7L:
144  Poles.resize(3);
145  Poles[0] = -0.53528043079643816554240378168164607183392315234269;
146  Poles[1] = -0.12255461519232669051527226435935734360548654942730;
147  Poles[2] = -0.0091486948096082769285930216516478534156925639545994;
148  break;
149  case 8L:
150  Poles.resize(4);
151  Poles[0] = -0.57468690924876543053013930412874542429066157804125;
152  Poles[1] = -0.16303526929728093524055189686073705223476814550830;
153  Poles[2] = -0.023632294694844850023403919296361320612665920854629;
154  Poles[3] = -0.00015382131064169091173935253018402160762964054070043;
155  break;
156  case 9L:
157  Poles.resize(4);
158  Poles[0] = -0.60799738916862577900772082395428976943963471853991;
159  Poles[1] = -0.20175052019315323879606468505597043468089886575747;
160  Poles[2] = -0.043222608540481752133321142979429688265852380231497;
161  Poles[3] = -0.0021213069031808184203048965578486234220548560988624;
162  break;
163  default:
164  std::cerr<<"Invalid spline degree"<<std::endl;;
165  exit(1);
166  }
167 }
168 
169 
183 template<typename InputIterator>
184 inline
185 double BSplineInitialCausalCoefficient(InputIterator infirst,
186  InputIterator inlast,
187  double z,
188  double Tolerance)
189 
190 {
191  std::size_t DataLength = inlast-infirst;
192 
193  /* this initialization corresponds to mirror boundaries */
194  std::size_t Horizon = DataLength;
195  if (Tolerance > 0.0)
196  {
197  Horizon = (long)std::ceil(std::log(Tolerance) / std::log(std::fabs(z)));
198  }
199 
200  double Sum, zn, z2n, iz;
201  if (Horizon < DataLength)
202  {
203  /* accelerated loop */
204  zn = z;
205  Sum = *infirst;
206  InputIterator it = infirst+1;
207  for (std::size_t n = 1L; n < Horizon; n++,it++)
208  {
209  Sum += zn * (*it);
210  zn *= z;
211  }
212  return(Sum);
213  }
214  else
215  {
216  /* full loop */
217  zn = z;
218  iz = 1.0 / z;
219  z2n = std::pow(z, (double)(DataLength - 1L));
220  Sum = (*infirst) + z2n * (*(inlast-1));
221  z2n *= z2n * iz;
222  for (InputIterator it=infirst+1; it!=inlast-1; it++)
223  {
224  Sum += (zn + z2n) * (*it);
225  zn *= z;
226  z2n *= iz;
227  }
228  return(Sum / (1.0 - zn * zn));
229  }
230 }
231 
243 template<typename InputIterator>
244 inline
245 double
246 BSplineInitialAntiCausalCoefficient(InputIterator outlast,
247  const double z)
248 
249 { /* this initialization corresponds to mirror boundaries */
250  return ((z / (z * z - 1.0)) * (z * (*(outlast-1)) + (*outlast)));
251 }
281 template<typename InputIterator,
282  typename OutputIterator,
283  typename InputIterator2>
284 inline
285 void BSplineComputeInterpolationCoefficients(InputIterator infirst,
286  InputIterator inlast,
287  OutputIterator outfirst,
288  InputIterator2 Poles_first,
289  InputIterator2 Poles_last,
290  double Tolerance = std::numeric_limits<double>::epsilon())
291 {
292  double Lambda = 1.0;
293  //std::size_t k;
294 
295  std::size_t NbPoles = static_cast<std::size_t>(Poles_last - Poles_first);
296 
297  std::size_t length = (inlast-infirst);
298 
299  /* special case required by mirror boundaries */
300  if (length == 1L)
301  {
302  return;
303  }
304  /* compute the overall gain */
305  InputIterator2 it_poles = Poles_first;
306  for (std::size_t k = 0L; k < NbPoles; ++k, ++it_poles)
307  {
308  Lambda = Lambda * (1.0 - *it_poles) * (1.0 - 1.0 / *it_poles);
309  }
310  InputIterator it;
311  OutputIterator it2;
312  /* apply the gain */
313  for (it=infirst,it2=outfirst;it!=inlast;it++,it2++)
314  {
315  *it2 = (*it)*Lambda;
316  }
317 
318  // loop over all poles
319  it_poles = Poles_first;
320  for (std::size_t k = 0L; k < NbPoles; ++k, ++it_poles)
321  {
322  // causal initialization
323  *outfirst = Lambda*slipalgo::BSplineInitialCausalCoefficient(infirst,inlast, *it_poles, Tolerance);
324  // causal recursion
325  for (it=infirst+1,it2=outfirst+1;it!=inlast;it++,it2++)
326  {
327  *it2 += *it_poles * (*(it2-1));
328  }
329  it2--;
330  // anticausal initialization
331  *it2 = slipalgo::BSplineInitialAntiCausalCoefficient(it2, *it_poles);
332  // anticausal recursion
333  for (it2--;it2!=outfirst-1;it2--)
334  {
335  *it2 = *it_poles * ((*(it2+1)) - *it2);
336  }
337 
338  }
339 }
340 
362 template<typename InputIterator,
363  typename OutputIterator,
364  class Array>
365 inline
366 void BSplineComputeInterpolationCoefficients(InputIterator infirst,
367  InputIterator inlast,
368  OutputIterator outfirst,
369  const Array& Poles,
370  double Tolerance)
371 {
372  slipalgo::BSplineComputeInterpolationCoefficients(infirst,inlast,
373  outfirst,
374  Poles.begin(),Poles.end(),
375  Tolerance);
376 }
377 
390  template<typename RandomAccessIterator>
391  inline
392  void BSplineComputeInterpolationIndexes(const std::size_t spline_degree,
393  double x,
394  RandomAccessIterator index_first,
395  RandomAccessIterator index_last)
396  {
397  assert(static_cast<std::size_t>(index_last - index_first) == (spline_degree + 1));
398  long splinedegree=static_cast<long>(spline_degree);
399 
400  /* compute the interpolation indexes */
401  long i = 0;
402  if (splinedegree & 1L)
403  {
404  i = static_cast<long>(std::floor(x)) - splinedegree / 2L;
405  }
406  else
407  {
408  i = static_cast<long>(std::floor(x + 0.5)) - spline_degree / 2L;
409  }
410  for (; index_first != index_last; ++index_first)
411  {
412  *index_first = i++;
413  }
414  }
415 
429  template<typename RandomAccessIterator,
430  typename Array>
431  inline
432  void BSplineComputeInterpolationWeights(const std::size_t spline_degree,
433  double x,
434  RandomAccessIterator index_first,
435  RandomAccessIterator index_last,
436  Array& Weights)
437  {
438  double w = 0.0;
439  double w2 = 0.0;
440  double w4 = 0.0;
441  double t = 0.0;
442  double t0 = 0.0;
443  double t1 = 0.0;
444  long SplineDegree=static_cast<long>(spline_degree);
445 
446  switch (SplineDegree)
447  {
448  case 2L:
449  /* x */
450  Weights.resize(3);
451  w = x - static_cast<double>(index_first[1]);
452  Weights[1] = 3.0 / 4.0 - w * w;
453  Weights[2] = (1.0 / 2.0) * (w - Weights[1] + 1.0);
454  Weights[0] = 1.0 - Weights[1] - Weights[2];
455  break;
456  case 3L:
457  /* x */
458  Weights.resize(4);
459  w = x - static_cast<double>(index_first[1]);
460  Weights[3] = (1.0 / 6.0) * w * w * w;
461  Weights[0] = (1.0 / 6.0) + (1.0 / 2.0) * w * (w - 1.0) - Weights[3];
462  Weights[2] = w + Weights[0] - 2.0 * Weights[3];
463  Weights[1] = 1.0 - Weights[0] - Weights[2] - Weights[3];
464  break;
465  case 4L:
466  /* x */
467  Weights.resize(5);
468  w = x - static_cast<double>(index_first[2]);
469  w2 = w * w;
470  t = (1.0 / 6.0) * w2;
471  Weights[0] = 1.0 / 2.0 - w;
472  Weights[0] *= Weights[0];
473  Weights[0] *= (1.0 / 24.0) * Weights[0];
474  t0 = w * (t - 11.0 / 24.0);
475  t1 = 19.0 / 96.0 + w2 * (1.0 / 4.0 - t);
476  Weights[1] = t1 + t0;
477  Weights[3] = t1 - t0;
478  Weights[4] = Weights[0] + t0 + (1.0 / 2.0) * w;
479  Weights[2] = 1.0 - Weights[0] - Weights[1] - Weights[3] - Weights[4];
480 
481  break;
482  case 5L:
483  /* x */
484  Weights.resize(6);
485  w = x - static_cast<double>(index_first[2]);
486  w2 = w * w;
487  Weights[5] = (1.0 / 120.0) * w * w2 * w2;
488  w2 -= w;
489  w4 = w2 * w2;
490  w -= 1.0 / 2.0;
491  t = w2 * (w2 - 3.0);
492  Weights[0] = (1.0 / 24.0) * (1.0 / 5.0 + w2 + w4) - Weights[5];
493  t0 = (1.0 / 24.0) * (w2 * (w2 - 5.0) + 46.0 / 5.0);
494  t1 = (-1.0 / 12.0) * w * (t + 4.0);
495  Weights[2] = t0 + t1;
496  Weights[3] = t0 - t1;
497  t0 = (1.0 / 16.0) * (9.0 / 5.0 - t);
498  t1 = (1.0 / 24.0) * w * (w4 - w2 - 5.0);
499  Weights[1] = t0 + t1;
500  Weights[4] = t0 - t1;
501  break;
502  case 6L:
503  /* x */
504  Weights.resize(7);
505  w = x - static_cast<double>(index_first[3]);
506  Weights[0] = 1.0 / 2.0 - w;
507  Weights[0] *= Weights[0] * Weights[0];
508  Weights[0] *= Weights[0] / 720.0;
509  Weights[1] = (361.0 / 192.0 - w * (59.0 / 8.0 + w
510  * (-185.0 / 16.0 + w * (25.0 / 3.0 + w * (-5.0 / 2.0 + w)
511  * (1.0 / 2.0 + w))))) / 120.0;
512  Weights[2] = (10543.0 / 960.0 + w * (-289.0 / 16.0 + w
513  * (79.0 / 16.0 + w * (43.0 / 6.0 + w * (-17.0 / 4.0 + w
514  * (-1.0 + w)))))) / 48.0;
515  w2 = w * w;
516  Weights[3] = (5887.0 / 320.0 - w2 * (231.0 / 16.0 - w2
517  * (21.0 / 4.0 - w2))) / 36.0;
518  Weights[4] = (10543.0 / 960.0 + w * (289.0 / 16.0 + w
519  * (79.0 / 16.0 + w * (-43.0 / 6.0 + w * (-17.0 / 4.0 + w
520  * (1.0 + w)))))) / 48.0;
521  Weights[6] = 1.0 / 2.0 + w;
522  Weights[6] *= Weights[6] * Weights[6];
523  Weights[6] *= Weights[6] / 720.0;
524  Weights[5] = 1.0 - Weights[0] - Weights[1] - Weights[2] - Weights[3]
525  - Weights[4] - Weights[6];
526  break;
527  case 7L:
528  /* x */
529  Weights.resize(8);
530  w = x - static_cast<double>(index_first[3]);
531  Weights[0] = 1.0 - w;
532  Weights[0] *= Weights[0];
533  Weights[0] *= Weights[0] * Weights[0];
534  Weights[0] *= (1.0 - w) / 5040.0;
535  w2 = w * w;
536  Weights[1] = (120.0 / 7.0 + w * (-56.0 + w * (72.0 + w
537  * (-40.0 + w2 * (12.0 + w * (-6.0 + w)))))) / 720.0;
538  Weights[2] = (397.0 / 7.0 - w * (245.0 / 3.0 + w * (-15.0 + w
539  * (-95.0 / 3.0 + w * (15.0 + w * (5.0 + w
540  * (-5.0 + w))))))) / 240.0;
541  Weights[3] = (2416.0 / 35.0 + w2 * (-48.0 + w2 * (16.0 + w2
542  * (-4.0 + w)))) / 144.0;
543  Weights[4] = (1191.0 / 35.0 - w * (-49.0 + w * (-9.0 + w
544  * (19.0 + w * (-3.0 + w) * (-3.0 + w2))))) / 144.0;
545  Weights[5] = (40.0 / 7.0 + w * (56.0 / 3.0 + w * (24.0 + w
546  * (40.0 / 3.0 + w2 * (-4.0 + w * (-2.0 + w)))))) / 240.0;
547  Weights[7] = w2;
548  Weights[7] *= Weights[7] * Weights[7];
549  Weights[7] *= w / 5040.0;
550  Weights[6] = 1.0 - Weights[0] - Weights[1] - Weights[2] - Weights[3]
551  - Weights[4] - Weights[5] - Weights[7];
552  break;
553  case 8L:
554  /* x */
555  Weights.resize(9);
556  w = x - static_cast<double>(index_first[4]);
557  Weights[0] = 1.0 / 2.0 - w;
558  Weights[0] *= Weights[0];
559  Weights[0] *= Weights[0];
560  Weights[0] *= Weights[0] / 40320.0;
561  w2 = w * w;
562  Weights[1] = (39.0 / 16.0 - w * (6.0 + w * (-9.0 / 2.0 + w2)))
563  * (21.0 / 16.0 + w * (-15.0 / 4.0 + w * (9.0 / 2.0 + w
564  * (-3.0 + w)))) / 5040.0;
565  Weights[2] = (82903.0 / 1792.0 + w * (-4177.0 / 32.0 + w
566  * (2275.0 / 16.0 + w * (-487.0 / 8.0 + w * (-85.0 / 8.0 + w
567  * (41.0 / 2.0 + w * (-5.0 + w * (-2.0 + w)))))))) / 1440.0;
568  Weights[3] = (310661.0 / 1792.0 - w * (14219.0 / 64.0 + w
569  * (-199.0 / 8.0 + w * (-1327.0 / 16.0 + w * (245.0 / 8.0 + w
570  * (53.0 / 4.0 + w * (-8.0 + w * (-1.0 + w)))))))) / 720.0;
571  Weights[4] = (2337507.0 / 8960.0 + w2 * (-2601.0 / 16.0 + w2
572  * (387.0 / 8.0 + w2 * (-9.0 + w2)))) / 576.0;
573  Weights[5] = (310661.0 / 1792.0 - w * (-14219.0 / 64.0 + w
574  * (-199.0 / 8.0 + w * (1327.0 / 16.0 + w * (245.0 / 8.0 + w
575  * (-53.0 / 4.0 + w * (-8.0 + w * (1.0 + w)))))))) / 720.0;
576  Weights[7] = (39.0 / 16.0 - w * (-6.0 + w * (-9.0 / 2.0 + w2)))
577  * (21.0 / 16.0 + w * (15.0 / 4.0 + w * (9.0 / 2.0 + w
578  * (3.0 + w)))) / 5040.0;
579  Weights[8] = 1.0 / 2.0 + w;
580  Weights[8] *= Weights[8];
581  Weights[8] *= Weights[8];
582  Weights[8] *= Weights[8] / 40320.0;
583  Weights[6] = 1.0 - Weights[0] - Weights[1] - Weights[2] - Weights[3]
584  - Weights[4] - Weights[5] - Weights[7] - Weights[8];
585  break;
586  case 9L:
587  /* x */
588  Weights.resize(10);
589  w = x - static_cast<double>(index_first[4]);
590  Weights[0] = 1.0 - w;
591  Weights[0] *= Weights[0];
592  Weights[0] *= Weights[0];
593  Weights[0] *= Weights[0] * (1.0 - w) / 362880.0;
594  Weights[1] = (502.0 / 9.0 + w * (-246.0 + w * (472.0 + w
595  * (-504.0 + w * (308.0 + w * (-84.0 + w * (-56.0 / 3.0 + w
596  * (24.0 + w * (-8.0 + w))))))))) / 40320.0;
597  Weights[2] = (3652.0 / 9.0 - w * (2023.0 / 2.0 + w * (-952.0 + w
598  * (938.0 / 3.0 + w * (112.0 + w * (-119.0 + w * (56.0 / 3.0 + w
599  * (14.0 + w * (-7.0 + w))))))))) / 10080.0;
600  Weights[3] = (44117.0 / 42.0 + w * (-2427.0 / 2.0 + w * (66.0 + w
601  * (434.0 + w * (-129.0 + w * (-69.0 + w * (34.0 + w * (6.0 + w
602  * (-6.0 + w))))))))) / 4320.0;
603  w2 = w * w;
604  Weights[4] = (78095.0 / 63.0 - w2 * (700.0 + w2 * (-190.0 + w2
605  * (100.0 / 3.0 + w2 * (-5.0 + w))))) / 2880.0;
606  Weights[5] = (44117.0 / 63.0 + w * (809.0 + w * (44.0 + w
607  * (-868.0 / 3.0 + w * (-86.0 + w * (46.0 + w * (68.0 / 3.0 + w
608  * (-4.0 + w * (-4.0 + w))))))))) / 2880.0;
609  Weights[6] = (3652.0 / 21.0 - w * (-867.0 / 2.0 + w * (-408.0 + w
610  * (-134.0 + w * (48.0 + w * (51.0 + w * (-4.0 + w) * (-1.0 + w)
611  * (2.0 + w))))))) / 4320.0;
612  Weights[7] = (251.0 / 18.0 + w * (123.0 / 2.0 + w * (118.0 + w
613  * (126.0 + w * (77.0 + w * (21.0 + w * (-14.0 / 3.0 + w
614  * (-6.0 + w * (-2.0 + w))))))))) / 10080.0;
615  Weights[9] = w2 * w2;
616  Weights[9] *= Weights[9] * w / 362880.0;
617  Weights[8] = 1.0 - Weights[0] - Weights[1] - Weights[2] - Weights[3]
618  - Weights[4] - Weights[5] - Weights[6] - Weights[7] - Weights[9];
619  break;
620  default:
621  std::cout<<"Invalid spline degree\n"<<std::endl;
622  }
623  }
624 
655 template <class RandomAccessIterator>
656 inline
657 double BSplineInterpolatedValue(RandomAccessIterator coef_first,
658  RandomAccessIterator coef_last,
659  double x,
660  std::size_t SplineDegree)
661 
662 {
663 
664  long coef_size = static_cast<long>(coef_last-coef_first);
665  long Width2 = 2L * coef_size - 2L;
666 
667  /* compute the interpolation indexes */
668  slip::Array<long> Index(SplineDegree+1);
669  slipalgo::BSplineComputeInterpolationIndexes(SplineDegree,
670  x,
671  Index.begin(),
672  Index.end());
673  /* compute the interpolation weights */
674  slip::Array<double> Weights;
675  slipalgo::BSplineComputeInterpolationWeights(SplineDegree,
676  x,
677  Index.begin(),
678  Index.end(),
679  Weights);
680 
681  long splinedegree = static_cast<long>(SplineDegree);
682  /* apply the mirror boundary conditions */
683  for (long k = 0L; k <= splinedegree; ++k)
684  {
685  Index[k] = (coef_size == 1L) ? (0L) : ((Index[k] < 0L) ?
686  (-Index[k] - Width2 * ((-Index[k]) / Width2))
687  : (Index[k] - Width2 * (Index[k] / Width2)));
688  if (coef_size <= Index[k])
689  {
690  Index[k] = Width2 - Index[k];
691  }
692  }
693 
694  /* perform interpolation */
695  double interpolated = 0.0;
696 
697  interpolated = 0.0;
698  for (long j = 0L; j <= splinedegree; ++j)
699  {
700  interpolated += Weights[j] * coef_first[Index[j]];
701  }
702 
703  return interpolated;
704 }
705 
706 
735 template <typename RandomAccessIterator2d1,
736  typename RandomAccessIterator2d2>
737 inline
738 void BSplineSamplesToCoefficients2d(RandomAccessIterator2d1 in_up,
739  RandomAccessIterator2d1 in_bot,
740  const std::size_t degree,
741  RandomAccessIterator2d2 out_up,
742  RandomAccessIterator2d2 out_bot,
743  const double tolerance =
744  std::numeric_limits<double>::epsilon())
745 {
746  assert((in_bot-in_up)[0] == (out_bot-out_up)[0]);
747  assert((in_bot-in_up)[1] == (out_bot-out_up)[1]);
748 
749  std::size_t rows = static_cast<std::size_t>((in_bot-in_up)[0]);
750  std::size_t cols = static_cast<std::size_t>((in_bot-in_up)[1]);
751 
752  slip::Array<double> poles;
753  slipalgo::BSplinePoles< slip::Array<double> >(degree,poles);
754 
755  /* convert the image samples into interpolation coefficients */
756  /* separable process, along x */
757  for (std::size_t y = 0; y < rows; ++y)
758  {
759  slipalgo::BSplineComputeInterpolationCoefficients(in_up.row_begin(y),
760  in_up.row_end(y),
761  out_up.row_begin(y),
762  poles,
763  tolerance);
764  }
765 
766  /* in-place separable process, along y */
767  slip::Array<double> tmp_col(rows);
768  for (std::size_t x = 0; x < cols; ++x)
769  {
770  std::copy(out_up.col_begin(x),out_up.col_end(x),tmp_col.begin());
771  slipalgo::BSplineComputeInterpolationCoefficients(tmp_col.begin(),tmp_col.end(),out_up.col_begin(x), poles,tolerance);
772 
773  }
774 }
775 
776 
777 
796 template <class Matrix1,
797  class Matrix2>
798 inline
799 void BSplineSamplesToCoefficients(const std::size_t degree,
800  const Matrix1 &input,
801  Matrix2 &output)
802 {
803 
804  slipalgo::BSplineSamplesToCoefficients2d(input.upper_left(),
805  input.bottom_right(),
806  degree,
807  output.upper_left(),
808  output.bottom_right());
809 }
810 
844 template <typename RandomAccessIterator3d1,
845  typename RandomAccessIterator3d2>
846 inline
847 void BSplineSamplesToCoefficients3d(RandomAccessIterator3d1 in_fup,
848  RandomAccessIterator3d1 in_bbot,
849  const std::size_t degree,
850  RandomAccessIterator3d2 out_fup,
851  RandomAccessIterator3d2 out_bbot,
852  const double tolerance =
853  std::numeric_limits<double>::epsilon())
854 {
855  assert((in_bbot-in_fup) == (out_bbot-out_fup));
856  std::size_t slices = static_cast<std::size_t>((in_bbot-in_fup)[0]);
857  std::size_t rows = static_cast<std::size_t>((in_bbot-in_fup)[1]);
858  std::size_t cols = static_cast<std::size_t>((in_bbot-in_fup)[2]);
859 
860  slip::Array<double> poles;
861  slipalgo::BSplinePoles< slip::Array<double> >(degree,poles);
862 
863  //slip::Array3d<double> tmp(slices,rows,cols);
864 
865  /* separable process, along z */
866  for (std::size_t x = 0; x < cols; ++x)
867  {
868  for (std::size_t y = 0; y < rows; ++y)
869  {
870  slipalgo::BSplineComputeInterpolationCoefficients(in_fup.slice_begin(y,x),in_fup.slice_end(y,x),out_fup.slice_begin(y,x), poles, tolerance);
871  }
872  }
873 
874  /* in-place separable process, along y */
875  slip::Array<double> tmp_row(rows);
876  for (std::size_t z = 0; z < slices; ++z)
877  {
878  for (std::size_t x = 0; x < cols; ++x)
879  {
880  slipalgo::BSplineComputeInterpolationCoefficients(out_fup.col_begin(z,x),out_fup.col_end(z,x),tmp_row.begin(), poles,tolerance);
881  std::copy(tmp_row.begin(),tmp_row.end(),out_fup.col_begin(z,x));
882  }
883  }
884 
885  /* in-place separable process, along x */
886  slip::Array<double> tmp_col(cols);
887  for (std::size_t z = 0; z < slices; ++z)
888  {
889  for (std::size_t y = 0; y < rows; ++y)
890  {
891  slipalgo::BSplineComputeInterpolationCoefficients(out_fup.row_begin(z,y),out_fup.row_end(z,y),tmp_col.begin(), poles, tolerance);
892  std::copy(tmp_col.begin(),tmp_col.end(),out_fup.row_begin(z,y));
893  }
894  }
895 }
896 
897 
911 template <typename RandomAccessIterator2d>
912 inline
913 double BSplineInterpolatedValue(RandomAccessIterator2d Bcoeff_up,
914  RandomAccessIterator2d Bcoeff_bot,
915  double x,
916  double y,
917  std::size_t SplineDegree)
918 
919 {
920  long rows = static_cast<long>((Bcoeff_bot - Bcoeff_up)[0]);
921  long cols = static_cast<long>((Bcoeff_bot - Bcoeff_up)[1]);
922 
923  long Width2 = 2L * cols - 2L;
924  long Height2 = 2L * rows - 2L;
925  long splinedegree = static_cast<long>(SplineDegree);
926 
927  /* compute the interpolation indexes */
928  slip::Array<long> xIndex(SplineDegree+1);
929  slipalgo::BSplineComputeInterpolationIndexes(SplineDegree,
930  x,
931  xIndex.begin(),
932  xIndex.end());
933  slip::Array<long> yIndex(SplineDegree+1);
934  slipalgo::BSplineComputeInterpolationIndexes(SplineDegree,
935  y,
936  yIndex.begin(),
937  yIndex.end());
938 
939 
940  /* compute the interpolation weights */
941  slip::Array<double> xWeights;
942  slipalgo::BSplineComputeInterpolationWeights(SplineDegree,
943  x,
944  xIndex.begin(),
945  xIndex.end(),
946  xWeights);
947  slip::Array<double> yWeights;
948  slipalgo::BSplineComputeInterpolationWeights(SplineDegree,
949  y,
950  yIndex.begin(),
951  yIndex.end(),
952  yWeights);
953 
954  /* apply the mirror boundary conditions */
955  for (long k = 0L; k <= splinedegree; k++)
956  {
957  xIndex[k] = (cols == 1L) ? (0L) : ((xIndex[k] < 0L) ?
958  (-xIndex[k] - Width2 * ((-xIndex[k]) / Width2))
959  : (xIndex[k] - Width2 * (xIndex[k] / Width2)));
960  if (cols <= xIndex[k])
961  {
962  xIndex[k] = Width2 - xIndex[k];
963  }
964  yIndex[k] = (rows == 1L) ? (0L) : ((yIndex[k] < 0L) ?
965  (-yIndex[k] - Height2 * ((-yIndex[k]) / Height2))
966  : (yIndex[k] - Height2 * (yIndex[k] / Height2)));
967  if (rows <= yIndex[k])
968  {
969  yIndex[k] = Height2 - yIndex[k];
970  }
971 
972  }
973 
974  /* perform interpolation */
975  double interpolated = 0.0;
976  slip::DPoint2d<int> d(0,0);
977  for (long j = 0L; j <= splinedegree; j++)
978  {
979  double w = 0.0;
980  d[0] = yIndex[j];
981  for (long i = 0L; i <= splinedegree; i++)
982  {
983  d[1] = xIndex[i];
984  w += xWeights[i] * Bcoeff_up[d];
985  }
986  interpolated += yWeights[j] * w;
987  }
988 
989  return(interpolated);
990 }
991 
1041 template <class Matrix1>
1042 inline
1043 double BSplineInterpolatedValue(Matrix1& Bcoeff,
1044  double x,
1045  double y,
1046  std::size_t SplineDegree)
1047 
1048 { /* begin InterpolatedValue */
1049 
1050 
1051  long Width2 = 2L * Bcoeff.cols() - 2L, Height2 = 2L * Bcoeff.rows() - 2L;
1052  long splinedegree = static_cast<long>(SplineDegree);
1053 
1054  /* compute the interpolation indexes */
1055  slip::Array<long> xIndex(SplineDegree+1);
1056  slipalgo::BSplineComputeInterpolationIndexes(SplineDegree,
1057  x,
1058  xIndex.begin(),
1059  xIndex.end());
1060  slip::Array<long> yIndex(SplineDegree+1);
1061  slipalgo::BSplineComputeInterpolationIndexes(SplineDegree,
1062  y,
1063  yIndex.begin(),
1064  yIndex.end());
1065 
1066 
1067  /* compute the interpolation weights */
1068  slip::Array<double> xWeights;
1069  slipalgo::BSplineComputeInterpolationWeights(SplineDegree,
1070  x,
1071  xIndex.begin(),
1072  xIndex.end(),
1073  xWeights);
1074  slip::Array<double> yWeights;
1075  slipalgo::BSplineComputeInterpolationWeights(SplineDegree,
1076  y,
1077  yIndex.begin(),
1078  yIndex.end(),
1079  yWeights);
1080 
1081  /* apply the mirror boundary conditions */
1082  for (long k = 0L; k <= splinedegree; k++)
1083  {
1084 
1085  xIndex[k] = (Bcoeff.cols() == 1L) ? (0L) : ((xIndex[k] < 0L) ?
1086  (-xIndex[k] - Width2 * ((-xIndex[k]) / Width2))
1087  : (xIndex[k] - Width2 * (xIndex[k] / Width2)));
1088  if (((long) Bcoeff.cols()) <= xIndex[k])
1089  {
1090  xIndex[k] = Width2 - xIndex[k];
1091  }
1092  yIndex[k] = (Bcoeff.rows() == 1L) ? (0L) : ((yIndex[k] < 0L) ?
1093  (-yIndex[k] - Height2 * ((-yIndex[k]) / Height2))
1094  : (yIndex[k] - Height2 * (yIndex[k] / Height2)));
1095  if (((long) Bcoeff.rows()) <= yIndex[k])
1096  {
1097  yIndex[k] = Height2 - yIndex[k];
1098  }
1099 
1100  }
1101 
1102  /* perform interpolation */
1103  double interpolated = 0.0;
1104  for (long j = 0L; j <= splinedegree; j++)
1105  {
1106  double w = 0.0;
1107  for (long i = 0L; i <= splinedegree; i++)
1108  {
1109  w += xWeights[i] * Bcoeff[yIndex[j]][xIndex[i]];
1110  }
1111  interpolated += yWeights[j] * w;
1112  }
1113 
1114  return(interpolated);
1115 }
1116 
1179 template <class Container3d>
1180 inline
1181 double BSplineInterpolatedValue(Container3d& Bcoeff,
1182  double x,
1183  double y,
1184  double z,
1185  std::size_t SplineDegree)
1186 
1187 {
1188 
1189  long Width2 = 2L * Bcoeff.cols() - 2L;
1190  long Height2 = 2L * Bcoeff.rows() - 2L;
1191  long Depth2 = 2L * Bcoeff.slices() - 2L;
1192 
1193  long splinedegree = static_cast<long>(SplineDegree);
1194 
1195  /* compute the interpolation indexes */
1196  slip::Array<long> xIndex(SplineDegree+1);
1197  slipalgo::BSplineComputeInterpolationIndexes(SplineDegree,
1198  x,
1199  xIndex.begin(),
1200  xIndex.end());
1201  slip::Array<long> yIndex(SplineDegree+1);
1202  slipalgo::BSplineComputeInterpolationIndexes(SplineDegree,
1203  y,
1204  yIndex.begin(),
1205  yIndex.end());
1206 
1207  slip::Array<long> zIndex(SplineDegree+1);
1208  slipalgo::BSplineComputeInterpolationIndexes(SplineDegree,
1209  z,
1210  zIndex.begin(),
1211  zIndex.end());
1212 
1213 
1214  /* compute the interpolation weights */
1215  slip::Array<double> xWeights;
1216  slipalgo::BSplineComputeInterpolationWeights(SplineDegree,
1217  x,
1218  xIndex.begin(),
1219  xIndex.end(),
1220  xWeights);
1221  slip::Array<double> yWeights;
1222  slipalgo::BSplineComputeInterpolationWeights(SplineDegree,
1223  y,
1224  yIndex.begin(),
1225  yIndex.end(),
1226  yWeights);
1227 
1228  slip::Array<double> zWeights;
1229  slipalgo::BSplineComputeInterpolationWeights(SplineDegree,
1230  z,
1231  zIndex.begin(),
1232  zIndex.end(),
1233  zWeights);
1234 
1235  /* apply the mirror boundary conditions */
1236  for (long k = 0L; k <= splinedegree; k++)
1237  {
1238 
1239  xIndex[k] = (Bcoeff.cols() == 1L) ? (0L) : ((xIndex[k] < 0L) ?
1240  (-xIndex[k] - Width2 * ((-xIndex[k]) / Width2))
1241  : (xIndex[k] - Width2 * (xIndex[k] / Width2)));
1242  if (((long) Bcoeff.cols()) <= xIndex[k])
1243  {
1244  xIndex[k] = Width2 - xIndex[k];
1245  }
1246  yIndex[k] = (Bcoeff.rows() == 1L) ? (0L) : ((yIndex[k] < 0L) ?
1247  (-yIndex[k] - Height2 * ((-yIndex[k]) / Height2))
1248  : (yIndex[k] - Height2 * (yIndex[k] / Height2)));
1249  if (((long) Bcoeff.rows()) <= yIndex[k])
1250  {
1251  yIndex[k] = Height2 - yIndex[k];
1252  }
1253  zIndex[k] = (Bcoeff.slices() == 1L) ? (0L) : ((zIndex[k] < 0L) ?
1254  (-zIndex[k] - Depth2 * ((-zIndex[k]) / Depth2))
1255  : (zIndex[k] - Depth2 * (zIndex[k] / Depth2)));
1256  if (((long) Bcoeff.slices()) <= zIndex[k])
1257  {
1258  zIndex[k] = Depth2 - zIndex[k];
1259  }
1260 
1261  }
1262 
1263  /* perform interpolation */
1264  double interpolated = 0.0;
1265  for (long j = 0L; j <= splinedegree; ++j)
1266  {
1267  double w = 0.0;
1268  for (long i = 0L; i <= splinedegree; ++i)
1269  {
1270  double w2 = 0.0;
1271  for (long k = 0L; k <= splinedegree; ++k)
1272  {
1273  w2 += zWeights[k] * Bcoeff[zIndex[k]][yIndex[j]][xIndex[i]];
1274  }
1275  w += xWeights[i] * w2;
1276  }
1277  interpolated += yWeights[j] * w;
1278  }
1279 
1280  return(interpolated);
1281 }
1282 
1296 template <typename RandomAccessIterator3d>
1297 inline
1298 double BSplineInterpolatedValue(RandomAccessIterator3d Bcoeff_fup,
1299  RandomAccessIterator3d Bcoeff_bbot,
1300  double x,
1301  double y,
1302  double z,
1303  std::size_t SplineDegree)
1304 
1305 {
1306  long slices = static_cast<long>((Bcoeff_bbot - Bcoeff_fup)[0]);
1307  long rows = static_cast<long>((Bcoeff_bbot - Bcoeff_fup)[1]);
1308  long cols = static_cast<long>((Bcoeff_bbot - Bcoeff_fup)[2]);
1309 
1310 
1311  long Width2 = 2L * cols - 2L;
1312  long Height2 = 2L * rows - 2L;
1313  long Depth2 = 2L * slices - 2L;
1314 
1315  long splinedegree = static_cast<long>(SplineDegree);
1316 
1317  /* compute the interpolation indexes */
1318  slip::Array<long> xIndex(SplineDegree+1);
1319  slipalgo::BSplineComputeInterpolationIndexes(SplineDegree,
1320  x,
1321  xIndex.begin(),
1322  xIndex.end());
1323  slip::Array<long> yIndex(SplineDegree+1);
1324  slipalgo::BSplineComputeInterpolationIndexes(SplineDegree,
1325  y,
1326  yIndex.begin(),
1327  yIndex.end());
1328 
1329  slip::Array<long> zIndex(SplineDegree+1);
1330  slipalgo::BSplineComputeInterpolationIndexes(SplineDegree,
1331  z,
1332  zIndex.begin(),
1333  zIndex.end());
1334 
1335 
1336  /* compute the interpolation weights */
1337  slip::Array<double> xWeights;
1338  slipalgo::BSplineComputeInterpolationWeights(SplineDegree,
1339  x,
1340  xIndex.begin(),
1341  xIndex.end(),
1342  xWeights);
1343  slip::Array<double> yWeights;
1344  slipalgo::BSplineComputeInterpolationWeights(SplineDegree,
1345  y,
1346  yIndex.begin(),
1347  yIndex.end(),
1348  yWeights);
1349 
1350  slip::Array<double> zWeights;
1351  slipalgo::BSplineComputeInterpolationWeights(SplineDegree,
1352  z,
1353  zIndex.begin(),
1354  zIndex.end(),
1355  zWeights);
1356 
1357  /* apply the mirror boundary conditions */
1358  for (long k = 0L; k <= splinedegree; k++)
1359  {
1360 
1361  xIndex[k] = (cols == 1L) ? (0L) : ((xIndex[k] < 0L) ?
1362  (-xIndex[k] - Width2 * ((-xIndex[k]) / Width2))
1363  : (xIndex[k] - Width2 * (xIndex[k] / Width2)));
1364  if (cols <= xIndex[k])
1365  {
1366  xIndex[k] = Width2 - xIndex[k];
1367  }
1368  yIndex[k] = (rows == 1L) ? (0L) : ((yIndex[k] < 0L) ?
1369  (-yIndex[k] - Height2 * ((-yIndex[k]) / Height2))
1370  : (yIndex[k] - Height2 * (yIndex[k] / Height2)));
1371  if (rows <= yIndex[k])
1372  {
1373  yIndex[k] = Height2 - yIndex[k];
1374  }
1375  zIndex[k] = (slices == 1L) ? (0L) : ((zIndex[k] < 0L) ?
1376  (-zIndex[k] - Depth2 * ((-zIndex[k]) / Depth2))
1377  : (zIndex[k] - Depth2 * (zIndex[k] / Depth2)));
1378  if (slices <= zIndex[k])
1379  {
1380  zIndex[k] = Depth2 - zIndex[k];
1381  }
1382 
1383  }
1384 
1385  /* perform interpolation */
1386  double interpolated = 0.0;
1387  slip::DPoint3d<int> d(0,0,0);
1388  for (long j = 0L; j <= splinedegree; ++j)
1389  {
1390  double w = 0.0;
1391  d[1] = yIndex[j];
1392  for (long i = 0L; i <= splinedegree; ++i)
1393  {
1394  double w2 = 0.0;
1395  d[2] = xIndex[i];
1396  for (long k = 0L; k <= splinedegree; ++k)
1397  {
1398  d[0] = zIndex[k];
1399  w2 += zWeights[k] * Bcoeff_fup[d];
1400  }
1401  w += xWeights[i] * w2;
1402  }
1403  interpolated += yWeights[j] * w;
1404  }
1405 
1406  return(interpolated);
1407 }
1408 
1409 
1434 template <typename RandomAccessIterator1,
1435  typename RandomAccessIterator2>
1436 void bspline_resampling_1d(RandomAccessIterator1 in_first,
1437  RandomAccessIterator1 in_last,
1438  const std::size_t spline_degree,
1439  RandomAccessIterator2 out_first,
1440  RandomAccessIterator2 out_last)
1441 {
1442  assert(static_cast<std::size_t>(out_last-out_first) >= static_cast<std::size_t>(in_last-in_first));
1443  std::size_t in_size = static_cast<std::size_t>(in_last-in_first);
1444  std::size_t out_size = static_cast<std::size_t>(out_last-out_first);
1445  slip::Array<double> Poles;
1446  slipalgo::BSplinePoles(spline_degree,Poles);
1447  slip::Array<double> Coef(in_size);
1448  slipalgo::BSplineComputeInterpolationCoefficients(in_first,in_last,
1449  Coef.begin(),
1450  Poles,
1451  std::numeric_limits<double>::epsilon());
1452 
1453  slip::Array<double> x(out_size);
1454  slip::iota(x.begin(),x.end(),0.0,1.0);
1455  slip::range_fun_interab<double,double> funab(0.0,double(out_size-1),0.0,double(in_size-1));
1456  slip::change_dynamic(x.begin(),x.end(),x.begin(),funab);
1457 
1458  slip::Array<double>::iterator it_x = x.begin();
1459  for(; out_first != out_last; ++out_first,++it_x)
1460  {
1461  *out_first =
1462  slipalgo::BSplineInterpolatedValue(Coef.begin(),Coef.end(),
1463  *it_x,spline_degree);
1464  }
1465 }
1466 
1467 
1496 template <typename RandomAccessIterator2d1,
1497  typename RandomAccessIterator2d2>
1498 void bspline_resampling_2d(RandomAccessIterator2d1 in_up,
1499  RandomAccessIterator2d1 in_bot,
1500  const std::size_t spline_degree,
1501  RandomAccessIterator2d2 out_up,
1502  RandomAccessIterator2d2 out_bot)
1503 {
1504  assert((in_bot-in_up)[0] <= (out_bot-out_up)[0]);
1505  assert((in_bot-in_up)[1] <= (out_bot-out_up)[1]);
1506 
1507  std::size_t in_rows = static_cast<std::size_t>((in_bot-in_up)[0]);
1508  std::size_t in_cols = static_cast<std::size_t>((in_bot-in_up)[1]);
1509 
1510  std::size_t out_rows = static_cast<std::size_t>((out_bot-out_up)[0]);
1511  std::size_t out_cols = static_cast<std::size_t>((out_bot-out_up)[1]);
1512 
1513  slip::Array<double> Poles;
1514  slipalgo::BSplinePoles(spline_degree,Poles);
1515  slip::Array2d<double> Coef(in_rows,in_cols);
1516  slipalgo::BSplineSamplesToCoefficients2d(in_up,in_bot,
1517  spline_degree,
1518  Coef.upper_left(),
1519  Coef.bottom_right());
1520 
1521  slip::Array<double> x(out_cols);
1522  slip::iota(x.begin(),x.end(),0.0,1.0);
1523  slip::range_fun_interab<double,double> funx(0.0,double(out_cols-1),0.0,double(in_cols-1));
1524  slip::change_dynamic(x.begin(),x.end(),x.begin(),funx);
1525 
1526  slip::Array<double> y(out_rows);
1527  slip::iota(y.begin(),y.end(),0.0,1.0);
1528  slip::range_fun_interab<double,double> funy(0.0,double(out_rows-1),0.0,double(in_rows-1));
1529  slip::change_dynamic(y.begin(),y.end(),y.begin(),funy);
1530 
1531 
1532  slip::DPoint2d<int> d(0,0);
1533  for(std::size_t i = 0; i < out_rows; ++i)
1534  {
1535  d[0] = i;
1536  for(std::size_t j = 0; j < out_cols; ++j)
1537  {
1538  d[1] = j;
1539  out_up[d] =
1540  slipalgo::BSplineInterpolatedValue(Coef.upper_left(),
1541  Coef.bottom_right(),
1542  x[j],
1543  y[i],
1544  spline_degree);
1545 
1546  }
1547  }
1548 
1549 }
1550 
1551 
1552 
1586 template <typename RandomAccessIterator3d1,
1587  typename RandomAccessIterator3d2>
1588 void bspline_resampling_3d(RandomAccessIterator3d1 in_fup,
1589  RandomAccessIterator3d1 in_bbot,
1590  const std::size_t spline_degree,
1591  RandomAccessIterator3d2 out_fup,
1592  RandomAccessIterator3d2 out_bbot)
1593 {
1594  // assert(static_cast<std::size_t>(in_last-in_first)
1595  std::size_t in_slices = static_cast<std::size_t>((in_bbot-in_fup)[0]);
1596  std::size_t in_rows = static_cast<std::size_t>((in_bbot-in_fup)[1]);
1597  std::size_t in_cols = static_cast<std::size_t>((in_bbot-in_fup)[2]);
1598 
1599  std::size_t out_slices = static_cast<std::size_t>((out_bbot-out_fup)[0]);
1600  std::size_t out_rows = static_cast<std::size_t>((out_bbot-out_fup)[1]);
1601  std::size_t out_cols = static_cast<std::size_t>((out_bbot-out_fup)[2]);
1602 
1603  slip::Array<double> Poles;
1604  slipalgo::BSplinePoles(spline_degree,Poles);
1605  slip::Array3d<double> Coef(in_slices,in_rows,in_cols);
1606  slipalgo::BSplineSamplesToCoefficients3d(in_fup,in_bbot,
1607  spline_degree,
1608  Coef.front_upper_left(),Coef.back_bottom_right());
1609 
1610  slip::Array<double> x(out_cols);
1611  slip::iota(x.begin(),x.end(),0.0,1.0);
1612  slip::range_fun_interab<double,double> funx(0.0,double(out_cols-1),0.0,double(in_cols-1));
1613  slip::change_dynamic(x.begin(),x.end(),x.begin(),funx);
1614 
1615  slip::Array<double> y(out_rows);
1616  slip::iota(y.begin(),y.end(),0.0,1.0);
1617  slip::range_fun_interab<double,double> funy(0.0,double(out_rows-1),0.0,double(in_rows-1));
1618  slip::change_dynamic(y.begin(),y.end(),y.begin(),funy);
1619 
1620  slip::Array<double> z(out_slices);
1621  slip::iota(z.begin(),z.end(),0.0,1.0);
1622  slip::range_fun_interab<double,double> funz(0.0,double(out_slices-1),0.0,double(in_slices-1));
1623  slip::change_dynamic(z.begin(),z.end(),z.begin(),funz);
1624 
1625 
1626  slip::DPoint3d<int> d(0,0,0);
1627  for(std::size_t k = 0; k < out_slices; ++k)
1628  {
1629  d[0] = k;
1630  for(std::size_t i = 0; i < out_rows; ++i)
1631  {
1632  d[1] = i;
1633  for(std::size_t j = 0; j < out_cols; ++j)
1634  {
1635  d[2] = j;
1636  out_fup[d] =
1637  slipalgo::BSplineInterpolatedValue(Coef.front_upper_left(),
1638  Coef.back_bottom_right(),
1639  x[j],
1640  y[i],
1641  z[k],
1642  spline_degree);
1643 
1644  }
1645  }
1646  }
1647 }
1648 
1649 
1675 template <class GrayScaleImage>
1676 inline
1677 void BSplineGeometricTransformation (const GrayScaleImage &InputImage,
1678  GrayScaleImage& OutputImage,
1679  const int SplineDegree=3,
1680  const double OriginX=0.0,
1681  const double OriginY=0.0,
1682  const double Angle=0.0,
1683  const double ShiftX=0.0,
1684  const double ShiftY=0.0,
1685  const bool Masking=true)
1686 
1687 {
1688  std::size_t Width = OutputImage.rows();
1689  std::size_t Height = OutputImage.cols();
1690  double width = static_cast<double>(Width) - 0.5;
1691  double height = static_cast<double>(Height) - 0.5;
1692 
1693  /* convert between a representation based on image samples */
1694  /* and a representation based on image B-spline coefficients */
1695  slip::Array2d<double> Coefficients(InputImage.rows(), InputImage.cols(),0.0);
1696  slipalgo::BSplineSamplesToCoefficients2d(InputImage.upper_left(),
1697  InputImage.bottom_right(),
1698  SplineDegree,
1699  Coefficients.upper_left(),
1700  Coefficients.bottom_right());
1701  /* prepare the geometry */
1702  double angle = Angle * slip::constants<double>::pi() / 180.0;
1703  double a11 = std::cos(angle);
1704  double a12 = -std::sin(angle);
1705  double a21 = -a12;
1706  double a22 = a11;
1707  double x0 = a11 * (ShiftX + OriginX) + a12 * (ShiftY + OriginY);
1708  double y0 = a21 * (ShiftX + OriginX) + a22 * (ShiftY + OriginY);
1709  double shiftX = OriginX - x0;
1710  double shiftY = OriginY - y0;
1711 
1712  /* visit all pixels of the output image and assign their value */
1713  // p = OutputImage;
1714 
1715  for (std::size_t y = 0; y < Height; ++y)
1716  {
1717  x0 = a12 * (double)y + shiftX;
1718  y0 = a22 * (double)y + shiftY;
1719  for (std::size_t x = 0; x < Width; ++x)
1720  {
1721  double x1 = x0 + a11 * static_cast<double>(x);
1722  double y1 = y0 + a21 * static_cast<double>(x);
1723 
1724  if (Masking)
1725  {
1726  if ( (x1 <= -0.5) || (width <= x1)
1727  || (y1 <= -0.5) || (height <= y1))
1728  {
1729  OutputImage[y][x] = 0.0F;
1730  }
1731  else
1732  {
1733  OutputImage[y][x] =
1734  slipalgo::BSplineInterpolatedValue(Coefficients,
1735  x1,
1736  y1,
1737  SplineDegree);
1738  }
1739  }
1740  else
1741  {
1742  OutputImage[y][x] =
1743  slipalgo::BSplineInterpolatedValue(Coefficients,
1744  x1,
1745  y1,
1746  SplineDegree);
1747  }
1748  }
1749  }
1750 }
1751 
1752  /* @} */
1753 
1754 }//::slipalgo
1755 #endif //SLIP_BSPLINE_INTERPOLATION_HPP
Array2d.hpp
Provides a class to manipulate 2d dynamic and generic arrays.
slipalgo::bspline_resampling_3d
void bspline_resampling_3d(RandomAccessIterator3d1 in_fup, RandomAccessIterator3d1 in_bbot, const std::size_t spline_degree, RandomAccessIterator3d2 out_fup, RandomAccessIterator3d2 out_bbot)
Resample a 3d range using bspline interpolation.
Definition: bspline_interpolation.hpp:1588
slip::Array2d::upper_left
iterator2d upper_left()
Returns a read/write iterator2d that points to the first element of the Array2d. It points to the upp...
Definition: Array2d.hpp:2744
DPoint2d.hpp
Provides a class to modelize the difference of slip::Point2d.
slipalgo::BSplineInitialCausalCoefficient
double BSplineInitialCausalCoefficient(InputIterator infirst, InputIterator inlast, double z, double Tolerance)
BSpline computation of causal coefficients based on http://bigwww.epfl.ch/thevenaz/interpolation/.
Definition: bspline_interpolation.hpp:185
slip::Array3d::front_upper_left
iterator3d front_upper_left()
Returns a read/write iterator3d that points to the first element of the Array3d. It points to the fro...
Definition: Array3d.hpp:4882
slip::DPoint3d
Difference of Point3D class, specialization of DPoint<CoordType,DIM> with DIM = 3.
Definition: DenseVector3dField3d.hpp:122
slipalgo::BSplineInterpolatedValue
double BSplineInterpolatedValue(RandomAccessIterator coef_first, RandomAccessIterator coef_last, double x, std::size_t SplineDegree)
BSpline interpolation from the interpolation coefficients. based on http://bigwww.epfl.ch/thevenaz/interpolation/.
Definition: bspline_interpolation.hpp:657
slip::change_dynamic
void change_dynamic(InputIterator first, InputIterator last, OutputIterator result, UnaryOperation dynamic_fun)
Changes the dynamic of a container.
Definition: dynamic.hpp:357
slipalgo::BSplineInitialAntiCausalCoefficient
double BSplineInitialAntiCausalCoefficient(InputIterator outlast, const double z)
BSpline computation of anticausal coefficients based on http://bigwww.epfl.ch/thevenaz/interpolation/...
Definition: bspline_interpolation.hpp:246
slip::Array3d
This is a three-dimensional dynamic and generic container. This container statisfies the Bidirectionn...
Definition: Array3d.hpp:109
slip::epsilon
slip::lin_alg_traits< T >::value_type epsilon()
Returns the epsilon value of a real or a complex.
Definition: linear_algebra_traits.hpp:211
slip::Array::begin
iterator begin()
Returns a read/write iterator that points to the first element in the Array. Iteration is done in ord...
Definition: Array.hpp:1077
slip::Array2d::bottom_right
iterator2d bottom_right()
Returns a read/write iterator2d that points to the past the end element of the Array2d. It points to past the end element of the bottom right element of the Array2d.
Definition: Array2d.hpp:2759
slipalgo::BSplineComputeInterpolationCoefficients
void BSplineComputeInterpolationCoefficients(InputIterator infirst, InputIterator inlast, OutputIterator outfirst, InputIterator2 Poles_first, InputIterator2 Poles_last, double Tolerance=std::numeric_limits< double >::epsilon())
BSpline computation of coefficients based on http://bigwww.epfl.ch/thevenaz/interpolation/.
Definition: bspline_interpolation.hpp:285
DPoint3d.hpp
Provides a class to modelize the difference of slip::Point3d.
slipalgo::bspline_resampling_1d
void bspline_resampling_1d(RandomAccessIterator1 in_first, RandomAccessIterator1 in_last, const std::size_t spline_degree, RandomAccessIterator2 out_first, RandomAccessIterator2 out_last)
Resample a range using bspline interpolation.
Definition: bspline_interpolation.hpp:1436
slip::iota
void iota(ForwardIterator first, ForwardIterator last, T value, T step=T(1))
Iota assigns sequential increasing values based on a predefined step to a range. That is...
Definition: arithmetic_op.hpp:2087
slip::constants::pi
static _Tp pi()
Constant .
Definition: macros.hpp:402
slip::sin
HyperVolume< T > sin(const HyperVolume< T > &M)
Definition: HyperVolume.hpp:5101
Array.hpp
Provides a class to manipulate 1d dynamic and generic arrays.
slip::cos
HyperVolume< T > cos(const HyperVolume< T > &M)
Definition: HyperVolume.hpp:5083
slipalgo::bspline_resampling_2d
void bspline_resampling_2d(RandomAccessIterator2d1 in_up, RandomAccessIterator2d1 in_bot, const std::size_t spline_degree, RandomAccessIterator2d2 out_up, RandomAccessIterator2d2 out_bot)
Resample a 2d range using bspline interpolation.
Definition: bspline_interpolation.hpp:1498
slip::Array::iterator
pointer iterator
Definition: Array.hpp:158
slipalgo::BSplineSamplesToCoefficients2d
void BSplineSamplesToCoefficients2d(RandomAccessIterator2d1 in_up, RandomAccessIterator2d1 in_bot, const std::size_t degree, RandomAccessIterator2d2 out_up, RandomAccessIterator2d2 out_bot, const double tolerance=std::numeric_limits< double >::epsilon())
Convert a 2d samples range to a 2d bspline coefficients range.
Definition: bspline_interpolation.hpp:738
slip::Array::end
iterator end()
Returns a read/write iterator that points one past the last element in the Array. Iteration is done i...
Definition: Array.hpp:1084
slip::copy
void copy(_II first, _II last, _OI output_first)
Copy algorithm optimized for slip iterators.
Definition: copy_ext.hpp:177
slip::DPoint2d
Difference of Point2D class, specialization of DPoint<CoordType,DIM> with DIM = 2.
Definition: Array2d.hpp:129
macros.hpp
Provides some mathematical functors and constants.
slip::sqrt
HyperVolume< T > sqrt(const HyperVolume< T > &M)
Definition: HyperVolume.hpp:5074
slipalgo::BSplineComputeInterpolationWeights
void BSplineComputeInterpolationWeights(const std::size_t spline_degree, double x, RandomAccessIterator index_first, RandomAccessIterator index_last, Array &Weights)
Computes the bspline interpolations weights. based on http://bigwww.epfl.ch/thevenaz/interpolation/.
Definition: bspline_interpolation.hpp:432
slipalgo::BSplineGeometricTransformation
void BSplineGeometricTransformation(const GrayScaleImage &InputImage, GrayScaleImage &OutputImage, const int SplineDegree=3, const double OriginX=0.0, const double OriginY=0.0, const double Angle=0.0, const double ShiftX=0.0, const double ShiftY=0.0, const bool Masking=true)
BSpline computation of a geometric transformation. based on http://bigwww.epfl.ch/thevenaz/interpolat...
Definition: bspline_interpolation.hpp:1677
arithmetic_op.hpp
Provides some algorithms to computes arithmetical operations on ranges.
slip::range_fun_interab
Functor object uses to change the range of container according to an affine transformation.
Definition: dynamic.hpp:129
slip::Array3d::back_bottom_right
iterator3d back_bottom_right()
Returns a read/write iterator3d that points to the past the end element of the Array3d. It points to past the end element of the back bottom right element of the Array3d.
Definition: Array3d.hpp:4897
slip::log
HyperVolume< T > log(const HyperVolume< T > &M)
Definition: HyperVolume.hpp:5146
dynamic.hpp
Provides algorithms and functors to change the dynamic of ranges.
Array3d.hpp
Provides a class to manipulate 3d dynamic and generic arrays.
slipalgo::BSplineSamplesToCoefficients3d
void BSplineSamplesToCoefficients3d(RandomAccessIterator3d1 in_fup, RandomAccessIterator3d1 in_bbot, const std::size_t degree, RandomAccessIterator3d2 out_fup, RandomAccessIterator3d2 out_bbot, const double tolerance=std::numeric_limits< double >::epsilon())
Convert a 3d samples range to a 3d bspline coefficients range.
Definition: bspline_interpolation.hpp:847
slipalgo::BSplineComputeInterpolationIndexes
void BSplineComputeInterpolationIndexes(const std::size_t spline_degree, double x, RandomAccessIterator index_first, RandomAccessIterator index_last)
Computes the bspline interpolations indices. based on http://bigwww.epfl.ch/thevenaz/interpolation/.
Definition: bspline_interpolation.hpp:392
slip::Array
This is a linear (one-dimensional) dynamic template container. This container statisfies the RandomAc...
Definition: Array.hpp:94
slipalgo::BSplinePoles
void BSplinePoles(const std::size_t degree, Array &Poles)
BSpline Poles computations based on http://bigwww.epfl.ch/thevenaz/interpolation/.
Definition: bspline_interpolation.hpp:110
slipalgo::BSplineSamplesToCoefficients
void BSplineSamplesToCoefficients(const std::size_t degree, const Matrix1 &input, Matrix2 &output)
BSpline computation of coefficients for 2D container based on http://bigwww.epfl.ch/thevenaz/interpol...
Definition: bspline_interpolation.hpp:799
slip::Array2d
This is a two-dimensional dynamic and generic container. This container statisfies the Bidirectionnal...
Definition: Array2d.hpp:135
Generated on Tue Sep 30 2014 10:05:56 for SLIP by   doxygen 1.8.6