SLIP  1.4
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
convolution.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 
68 
75 #ifndef SLIP_CONVOLUTION_HPP
76 #define SLIP_CONVOLUTION_HPP
77 
78 
79 #include <algorithm>
80 #include <numeric>
81 #include <cassert>
82 #include <iterator>
83 #include <vector>
84 #include "arithmetic_op.hpp"
85 #include "border_treatment.hpp"
86 #include "Array.hpp"
87 #include "Array3d.hpp"
88 #include "Box2d.hpp"
89 #include "Vector.hpp"
90 #include "FFT.hpp"
91 #include "apply.hpp"
92 
93 namespace slip
94 {
95 
97  /* @{ */
132  template <typename SrcIter, typename KernelIter, typename ResIter>
133  inline
134  void valid_convolution(SrcIter first,
135  SrcIter last,
136  KernelIter kernel_first,
137  KernelIter kernel_last,
138  std::size_t ksize_left,
139  std::size_t ksize_right,
140  ResIter result)
141  {
142  assert(static_cast<std::size_t>(kernel_last - kernel_first) == (ksize_left + ksize_right + 1));
143  SrcIter first_offset = first + ksize_left;
144  SrcIter last_offset = last - ksize_right;
145  std::size_t K = ksize_left + ksize_right + 1;
146  for(;first_offset < last_offset; ++first_offset, ++result)
147  {
148  SrcIter first_offset_tmp = first_offset - ksize_left;
149  SrcIter last_offset_tmp = first_offset_tmp + K;
150  std::reverse_iterator<KernelIter> kernel_last_tmp =
151  std::reverse_iterator<KernelIter>(kernel_last);
152 
153  *result = std::inner_product(first_offset_tmp,last_offset_tmp,
154  kernel_last_tmp,
155  typename std::iterator_traits<ResIter>::value_type());
156 
157  }
158  }
159 
160 
190  template <typename SrcIter, typename KernelIter, typename ResIter>
191  inline
192  void full_convolution(SrcIter first,
193  SrcIter last,
194  KernelIter kernel_first,
195  KernelIter kernel_last,
196  ResIter result)
197  {
198 
199  std::size_t K = static_cast<std::size_t>(kernel_last - kernel_first - 1);
200 
201  //valid convolution
202  slip::valid_convolution(first,last,
203  kernel_first,kernel_last,
204  K,0,
205  result + K);
206 
207  //zero-padded treatment
208  //left
209  std::reverse_iterator<KernelIter> rkernel_last_tmp =
210  std::reverse_iterator<KernelIter>(kernel_last);
211 
212  for(size_t i = 0; i < K; ++i)
213  {
214  *(result + i ) = std::inner_product(first,first + i + 1,rkernel_last_tmp + (K - i),typename std::iterator_traits<ResIter>::value_type());
215  }
216 
217  //right
218  SrcIter first_tmp = last - K;
219  ResIter result_tmp = result + (last - first);
220  rkernel_last_tmp = std::reverse_iterator<KernelIter>(kernel_last);
221 
222  for(size_t i = 0; i < K; ++i)
223  {
224  *(result_tmp + i) = std::inner_product(first_tmp + i,last,rkernel_last_tmp,typename std::iterator_traits<ResIter>::value_type());
225  }
226  }
227 
297  template <typename SrcIter, typename KernelIter, typename ResIter>
298  inline
299  void same_convolution(SrcIter first,
300  SrcIter last,
301  KernelIter kernel_first,
302  KernelIter kernel_last,
303  std::size_t ksize_left,
304  std::size_t ksize_right,
305  ResIter result,
306  slip::BORDER_TREATMENT bt = slip::BORDER_TREATMENT_ZERO_PADDED)
307  {
308  //valid convolution
309  slip::valid_convolution(first,last,
310  kernel_first,kernel_last,
311  ksize_left,ksize_right,
312  result + ksize_right);
313 
314  //border_treatments
315  if(bt == slip::BORDER_TREATMENT_ZERO_PADDED)
316  {
317  std::size_t K = ksize_left + ksize_right;
318  //left
319  for(size_t i = 0; i < ksize_right; ++i)
320  {
321  *(result + (ksize_right - i - 1)) = std::inner_product(first,first + (K - i),std::reverse_iterator<KernelIter>(kernel_last - i - 1),typename std::iterator_traits<ResIter>::value_type());
322 
323  }
324  //right
325  SrcIter first_tmp = last - K;
326  ResIter result_tmp = result + ((last - first) - ksize_left);
327  for(size_t i = 0; i < ksize_left; ++i)
328  {
329  *(result_tmp + i) = std::inner_product(first_tmp + i,last,std::reverse_iterator<KernelIter>(kernel_last),typename std::iterator_traits<ResIter>::value_type());
330 
331  }
332 
333 
334  }
335  else if(bt == slip::BORDER_TREATMENT_NEUMANN)
336  {
337  std::size_t K = ksize_left + ksize_right;
338  typedef typename std::iterator_traits<ResIter>::value_type res_value_type;
339  typedef typename std::iterator_traits<SrcIter>::value_type src_value_type;
340  typedef typename std::iterator_traits<KernelIter>::value_type k_value_type;
341 
342 
343  //left
344  //compute the left accumulate sum of the elements of the reverse
345  //kernel and store each partial result in accu_vec
346  // k_value_type accum = k_value_type(0);
347  std::reverse_iterator<KernelIter> rkernel =
348  std::reverse_iterator<KernelIter>(kernel_last);
349 
350  std::vector<k_value_type> accu_vec(ksize_right);
351  std::partial_sum(rkernel,rkernel + ksize_right,accu_vec.begin());
352 
353  res_value_type left_bvalue = src_value_type(*first);
354  for(size_t i = 0; i < ksize_right; ++i)
355  {
356 
357  *(result + (ksize_right - i - 1)) =
358  std::inner_product(first,
359  first + (K - i),
360  std::reverse_iterator<KernelIter>(kernel_last - i - 1),
361  res_value_type(accu_vec[i]) * left_bvalue);
362  }
363 
364  //right
365  //compute the right accumulate sum of the elements of the reverse
366  //kernel and store each partial result in accu_vec2
367  //accum = k_value_type(0);
368  KernelIter rkernel2 = kernel_last - K - 1;
369 
370  std::vector<k_value_type> accu_vec2(ksize_left);
371  std::partial_sum(rkernel2,rkernel2 + ksize_left,accu_vec2.begin());
372  res_value_type right_bvalue = src_value_type(*(last - 1));
373  SrcIter first_tmp = last - K;
374  ResIter result_tmp = result + ((last - first) - ksize_left);
375  for(size_t i = 0; i < ksize_left; ++i)
376  {
377  *(result_tmp + i) =
378  std::inner_product(first_tmp + i,
379  last,
380  std::reverse_iterator<KernelIter>(kernel_last),
381  res_value_type(accu_vec2[i]) * right_bvalue);
382  }
383 
384  }
385  else if(bt == slip::BORDER_TREATMENT_AVOID)
386  {
387  //left
388  for(size_t i = 0; i < ksize_right; ++i)
389  {
390  *(result + (ksize_right - i - 1)) =
391  typename std::iterator_traits<ResIter>::value_type();
392  }
393  //right
394  ResIter result_tmp = result + ((last - first) - ksize_left);
395  for(size_t i = 0; i < ksize_left; ++i)
396  {
397  *(result_tmp + i) =
398  typename std::iterator_traits<ResIter>::value_type();
399  }
400  }
401 
402  else if(bt == slip::BORDER_TREATMENT_CIRCULAR)
403  {
404  typedef typename std::iterator_traits<SrcIter>::value_type src_value_type;
405 
406  std::size_t kernel_size = ksize_right + ksize_left + 1;
407  std::size_t tmp_size = ksize_right + ksize_right+ksize_left;
408  slip::Array<src_value_type> tmp_data(tmp_size);
409 
410  // Copy data :
411  for(size_t i = 0; i < tmp_size-ksize_right; ++i)
412  {
413  tmp_data[i+ksize_right]=*(first+i);
414  }
415  for(size_t i = 0; i < ksize_right; ++i)
416  {
417  tmp_data[ksize_right-i-1]=*(last-i-1);
418  }
419 
420  //left
421  for(size_t i = 0; i < ksize_right; ++i)
422  {
423  *(result + i)=std::inner_product(tmp_data.begin() + i,tmp_data.begin() + i + kernel_size,std::reverse_iterator<KernelIter>(kernel_last),typename std::iterator_traits<ResIter>::value_type());
424  }
425 
426  std::size_t tmp_size2 = ksize_left + ksize_left+ksize_right;
427  slip::Array<src_value_type> tmp_data2(tmp_size2);
428  // Copy data :
429  for(size_t i = 0; i < tmp_size2-ksize_left; ++i)
430  {
431  tmp_data2[ksize_right+ksize_left-i-1]=*(last-i-1);
432  }
433  for(size_t i = 0; i < ksize_left; ++i)
434  {
435  tmp_data2[ksize_right+ksize_left+i]=*(first+i);
436  }
437  //right
438  for(size_t i = 0; i < ksize_left; ++i)
439  {
440  *(result + ((last - first) - ksize_left) + i)= std::inner_product(tmp_data2.begin() + i,tmp_data2.begin() + i + kernel_size,std::reverse_iterator<KernelIter>(kernel_last),typename std::iterator_traits<ResIter>::value_type());
441  }
442  }
443  else
444  {
445  std::size_t K = ksize_left + ksize_right;
446  //left
447  for(size_t i = 0; i < ksize_right; ++i)
448  {
449  *(result + (ksize_right - i - 1)) = std::inner_product(first,first + (K - i),std::reverse_iterator<KernelIter>(kernel_last - i - 1),typename std::iterator_traits<ResIter>::value_type());
450 
451  }
452  //right
453  SrcIter first_tmp = last - K;
454  ResIter result_tmp = result + ((last - first) - ksize_left);
455  for(size_t i = 0; i < ksize_left; ++i)
456  {
457  *(result_tmp + i) = std::inner_product(first_tmp + i,last,std::reverse_iterator<KernelIter>(kernel_last),typename std::iterator_traits<ResIter>::value_type());
458 
459  }
460 
461  }
462  }
541 template <typename SrcIter, typename KernelIter, typename ResIter>
542  inline
543  void same_convolution(SrcIter first,
544  SrcIter last,
545  KernelIter kernel_first,
546  KernelIter kernel_last,
547  ResIter result,
548  slip::BORDER_TREATMENT bt= slip::BORDER_TREATMENT_ZERO_PADDED)
549  {
550 
551  std::size_t ksize_left = (kernel_last - kernel_first)/ 2;
552  std::size_t ksize_right = (kernel_last - kernel_first - 1) - ksize_left;
553 
554  slip::same_convolution(first,
555  last,
556  kernel_first,
557  kernel_last,
558  ksize_left,
559  ksize_right,
560  result,
561  bt);
562  }
563 
600 template <typename Vector1, typename Vector2, typename Vector3>
601 
602 inline void same_convolution(const Vector1& Signal,
603  const Vector2& Kernel,
604  Vector3& Result,
605  slip::BORDER_TREATMENT bt= slip::BORDER_TREATMENT_ZERO_PADDED)
606 
607 
608  {
609 
610  slip::same_convolution(Signal.begin(),Signal.end(),Kernel.begin(),Kernel.end(),Result.begin(),bt);
611 
612  }
613 
648  template<typename RandomAccessIterator1,
649  typename RandomAccessIterator2,
650  typename RandomAccessIterator3>
651  inline
652  void fft_circular_convolution(RandomAccessIterator1 first , RandomAccessIterator1 last,
653  RandomAccessIterator2 first2, RandomAccessIterator2 last2,
654  RandomAccessIterator3 conv_first, RandomAccessIterator3 conv_last )
655  {
656  assert((last - first) == (last2 - first2));
657  assert((last - first) == (conv_last - conv_first));
658  assert((last - first) == (last2 - first2));
659 
660  typedef typename std::iterator_traits<RandomAccessIterator1>::value_type value_type;
661  typedef typename slip::lin_alg_traits<value_type>::value_type T;
662 
663 
664  std::size_t N = static_cast<std::size_t>(last-first);
665  slip::Array<std::complex<T> > fft1(N);
666  slip::Array<std::complex<T> > fft2(N);
667  slip::real_fft(first,last,fft1.begin());
668  slip::real_fft(first2,last2,fft2.begin());
669  slip::multiplies(fft1.begin(),fft1.end(),
670  fft2.begin(),
671  fft1.begin());
672  slip::ifft(fft1.begin(),fft1.end(),fft2.begin());
673  std::transform(fft2.begin(),fft2.end(),conv_first,
674  slip::un_real<std::complex<T>,T>());
675 
676 
677  }
678 
679 
708  template<typename RandomAccessIterator1, typename RandomAccessIterator2, typename RandomAccessIterator3>
709  inline
710  void fft_convolution(RandomAccessIterator1 first , RandomAccessIterator1 last,
711  RandomAccessIterator2 first2, RandomAccessIterator2 last2,
712  RandomAccessIterator3 conv_first, RandomAccessIterator3 conv_last )
713  {
714 
715  assert((last - first) != 0);
716  assert((last2 - first2) != 0);
717  assert((conv_last-conv_first) >= (last - first) + (last2 - first2) -1);
718  typedef typename std::iterator_traits<RandomAccessIterator1>::difference_type _Difference1;
719  _Difference1 size1 = last - first;
720  typedef typename std::iterator_traits<RandomAccessIterator2>::difference_type _Difference2;
721  _Difference2 size2 = last2 - first2;
722 
723  _Difference1 size = size1 + _Difference1(size2) - 1;
724 
725 
726  typedef typename std::iterator_traits<RandomAccessIterator1>::value_type value_type;
727  typedef typename slip::lin_alg_traits<value_type>::value_type Real;
728 
729  //zero padding of the input data
730  slip::Vector<value_type> data1(size,value_type());
731  std::copy(first,last,data1.begin());
732  slip::Vector<value_type> data2(size,value_type());
733  std::copy(first2,last2,data2.begin());
734 
735  //compute the two FFT
736  slip::Vector<std::complex<value_type> > fft1(size,std::complex<value_type>());
737  slip::Vector<std::complex<value_type> > fft2(size,std::complex<value_type>());
738  slip::real_fft(data1.begin(), data1.end(), fft1.begin());
739  slip::real_fft(data2.begin(), data2.end(), fft2.begin());
740  //multiplies the two FFT
741  slip::multiplies(fft1.begin(),fft1.end(),
742  fft2.begin(),
743  fft1.begin());
744  //invert the product of the two FFT
745  slip::ifft(fft1.begin(),fft1.end(),fft2.begin());
746  //get the real part of the inverse FFT
747  std::transform(fft2.begin(),fft2.end(),conv_first,slip::un_real<std::complex<Real>,Real>());
748  }
749 
772  template <typename Vector1, typename Vector2, typename Vector3>
773  inline void fft_convolution(const Vector1& Signal,
774  const Vector2& Kernel,
775  Vector3& Result)
776 
777 
778  {
779  slip::fft_convolution(Signal.begin(),Signal.end(),
780  Kernel.begin(),Kernel.end(),
781  Result.begin(),Result.end());
782  }
783 
784 
785  /* @} */
786 
787 
789  /* @{ */
837  template <typename RandomAccessIterator2d1,
838  typename KernelIter1,
839  typename KernelIter2,
840  typename RandomAccessIterator2d2>
841 inline
842 void separable_convolution2d(RandomAccessIterator2d1 I_up,
843  RandomAccessIterator2d1 I_bot,
844  KernelIter1 k1_first,
845  KernelIter1 k1_last,
846  std::size_t ksize_left1,
847  std::size_t ksize_right1,
848  slip::BORDER_TREATMENT bt1,
849  KernelIter2 k2_first,
850  KernelIter2 k2_last,
851  std::size_t ksize_left2,
852  std::size_t ksize_right2,
853  slip::BORDER_TREATMENT bt2,
854  RandomAccessIterator2d2 R_up)
855 {
856 
857  typename RandomAccessIterator2d1::difference_type size2d =
858  I_bot - I_up;
859  typedef typename RandomAccessIterator2d1::size_type size_type;
860  typedef typename RandomAccessIterator2d1::value_type value_type;
861  size_type nb_rows = size2d[0];
862  size_type nb_cols = size2d[1];
863 
864 
865  for(size_type i = 0; i < nb_rows; ++i)
866  {
867  slip::same_convolution(I_up.row_begin(i),I_up.row_end(i),
868  k1_first,k1_last,
869  ksize_left1,
870  ksize_right1,
871  R_up.row_begin(i),
872  bt1);
873  }
874  slip::Array<value_type> col_tmp(nb_rows);
875  for(size_type j = 0; j < nb_cols; ++j)
876  {
877  std::copy(R_up.col_begin(j),R_up.col_end(j),col_tmp.begin());
878  slip::same_convolution(col_tmp.begin(),col_tmp.end(),
879  k2_first,k2_last,
880  ksize_left2,
881  ksize_right2,
882  R_up.col_begin(j),
883  bt2);
884  }
885 
886 }
887 
888 
889 
926 template <typename RandomAccessIterator2d1,
927  typename KernelIter1,
928  typename KernelIter2,
929  typename RandomAccessIterator2d2>
930 inline
931 void separable_convolution2d(RandomAccessIterator2d1 I_up,
932  RandomAccessIterator2d1 I_bot,
933  KernelIter1 k1_first,
934  KernelIter1 k1_last,
935  KernelIter2 k2_first,
936  KernelIter2 k2_last,
937  RandomAccessIterator2d2 R_up,
938  slip::BORDER_TREATMENT bt = slip::BORDER_TREATMENT_ZERO_PADDED)
939 {
940  std::size_t size1 = (k1_last-k1_first)/2;
941  std::size_t size1_bis = (k1_last-k1_first) - 1 - size1;
942  std::size_t size2 = (k2_last-k2_first)/2;
943  std::size_t size2_bis = (k2_last-k2_first) - 1 - size2;
944 
945  slip::separable_convolution2d(I_up,I_bot,
946  k1_first,
947  k1_last,
948  size1,
949  size1_bis,
950  bt,
951  k2_first,
952  k2_last,
953  size2,
954  size2_bis,
955  bt,
956  R_up);
957 }
958 
991  template <typename Container1,
992  typename Kernel1,
993  typename Kernel2,
994  typename Container2>
995  inline void separable_convolution2d(const Container1& Signal,
996  const Kernel1& kernel1,
997  const Kernel2& kernel2,
998  Container2& Result,
999  slip::BORDER_TREATMENT bt = slip::BORDER_TREATMENT_ZERO_PADDED)
1000 
1001 
1002  {
1003  slip::separable_convolution2d(Signal.upper_left(),
1004  Signal.bottom_right(),
1005  kernel1.begin(),kernel1.end(),
1006  kernel2.begin(),kernel2.end(),
1007  Result.upper_left());
1008  }
1009 
1010 
1055 template <typename RandomAccessIterator2d1,
1056  typename KernelIter1,
1057  typename KernelIter2>
1058 inline
1059 void separable_convolution2d(RandomAccessIterator2d1 I_up,
1060  RandomAccessIterator2d1 I_bot,
1061  KernelIter1 k1_first,
1062  KernelIter1 k1_last,
1063  std::size_t ksize_left1,
1064  std::size_t ksize_right1,
1065  slip::BORDER_TREATMENT bt1,
1066  KernelIter2 k2_first,
1067  KernelIter2 k2_last,
1068  std::size_t ksize_left2,
1069  std::size_t ksize_right2,
1070  slip::BORDER_TREATMENT bt2)
1071 {
1072 
1073  typename RandomAccessIterator2d1::difference_type size2d =
1074  I_bot - I_up;
1075  typedef typename RandomAccessIterator2d1::size_type size_type;
1076  typedef typename RandomAccessIterator2d1::value_type value_type;
1077  const size_type nb_rows = size2d[0];
1078  const size_type nb_cols = size2d[1];
1079 
1080  //slip::Array2d<value_type> Atmp(nb_rows,nb_cols);
1081  slip::Array<value_type> row_tmp(nb_cols);
1082  for(size_type i = 0; i < nb_rows; ++i)
1083  {
1084  std::copy(I_up.row_begin(i),I_up.row_end(i),row_tmp.begin());
1085  slip::same_convolution(row_tmp.begin(),row_tmp.end(),
1086  k1_first,k1_last,
1087  ksize_left1,
1088  ksize_right1,
1089  I_up.row_begin(i),
1090  bt1);
1091  }
1092  slip::Array<value_type> col_tmp(nb_rows);
1093  for(size_type j = 0; j < nb_cols; ++j)
1094  {
1095  std::copy(I_up.col_begin(j),I_up.col_end(j),col_tmp.begin());
1096  slip::same_convolution(col_tmp.begin(),col_tmp.end(),
1097  k2_first,k2_last,
1098  ksize_left2,
1099  ksize_right2,
1100  I_up.col_begin(j),
1101  bt2);
1102  }
1103 
1104 }
1105 
1106 
1137 template <typename RandomAccessIterator2d1,
1138  typename KernelIter1,
1139  typename KernelIter2,
1140  typename RandomAccessIterator2d2>
1141 inline
1142 void separable_fft_convolution2d(RandomAccessIterator2d1 I_up,
1143  RandomAccessIterator2d1 I_bot,
1144  KernelIter1 k1_first,
1145  KernelIter1 k1_last,
1146  KernelIter2 k2_first,
1147  KernelIter2 k2_last,
1148  RandomAccessIterator2d2 R_up)
1149 {
1150  typename RandomAccessIterator2d1::difference_type size2d =
1151  I_bot - I_up;
1152  typedef typename RandomAccessIterator2d1::size_type size_type;
1153  typedef typename RandomAccessIterator2d1::value_type value_type;
1154  size_type nb_rows = size2d[0];
1155  size_type nb_cols = size2d[1];
1156  size_type size1 = (k1_last - k1_first);
1157  size_type size2 = (k2_last - k2_first);
1158  size_type size2o2 = size2/2;
1159  slip::Array2d<value_type> Atmp(nb_rows+size2-1,nb_cols+size1-1);
1160  slip::Array<value_type> Vtmp(Atmp.cols());
1161 
1162  for(size_t i = 0; i < nb_rows; ++i)
1163  {
1164  slip::fft_convolution(I_up.row_begin(i),I_up.row_end(i),
1165  k1_first,k1_last,
1166  Atmp.row_begin(i),Atmp.row_end(i));
1167  }
1168  for(size_t j = 0; j < nb_cols; ++j)
1169  {
1170  slip::fft_convolution(Atmp.col_begin(j)+size2o2,Atmp.col_end(j)-size2o2,
1171  k2_first,k2_last,
1172  Vtmp.begin(),Vtmp.end());
1173  std::copy(Vtmp.begin()+size2o2,Vtmp.end()-size2o2,R_up.col_begin(j));
1174  }
1175 }
1202  template <typename Container1,
1203  typename Kernel1,
1204  typename Kernel2,
1205  typename Container2>
1206  inline void separable_fft_convolution2d(const Container1& Signal,
1207  const Kernel1& kernel1,
1208  const Kernel2& kernel2,
1209  Container2& Result)
1210 
1211 
1212  {
1213  slip::separable_fft_convolution2d(Signal.upper_left(),
1214  Signal.bottom_right(),
1215  kernel1.begin(),kernel1.end(),
1216  kernel2.begin(),kernel2.end(),
1217  Result.upper_left());
1218  }
1219 
1220 
1269  template <typename SrcIter2d,
1270  typename KernelIter2d,
1271  typename ResIter2d>
1272  inline
1273  void valid_convolution2d(SrcIter2d S_up,
1274  SrcIter2d S_bot,
1275  KernelIter2d kernel_up,
1276  KernelIter2d kernel_bot,
1277  std::size_t ksize_left,
1278  std::size_t ksize_right,
1279  std::size_t ksize_up,
1280  std::size_t ksize_bot,
1281  ResIter2d R_up,
1282  ResIter2d R_bot)
1283  {
1284 
1285  assert((R_bot - R_up)[0] == ((S_bot - S_up)[0] - (kernel_bot - kernel_up)[0] + 1));
1286  assert((R_bot - R_up)[1] == ((S_bot - S_up)[1] - (kernel_bot - kernel_up)[1] + 1));
1287  assert(static_cast<int>((kernel_bot - kernel_up)[1]) == static_cast<int>(ksize_left + ksize_right + 1));
1288  assert(static_cast<int>((kernel_bot - kernel_up)[0]) == static_cast<int>(ksize_up + ksize_bot + 1));
1289  std::size_t k_rows = static_cast<std::size_t>((kernel_bot - kernel_up)[0]);
1290  std::size_t k_cols = static_cast<std::size_t>((kernel_bot - kernel_up)[1]);
1291  std::size_t r_rows = static_cast<std::size_t>((R_bot - R_up)[0]);
1292  std::size_t r_cols = static_cast<std::size_t>((R_bot - R_up)[1]);
1293 
1294 
1295  typedef typename SrcIter2d::value_type value_type;
1296  typedef typename KernelIter2d::value_type k_value_type;
1297 
1298  //temporary array used to compute the reverse kernel
1299  slip::DPoint2d<int> dp(1,0);
1300  std::reverse_iterator<KernelIter2d> k_rupper_left = std::reverse_iterator<KernelIter2d>(kernel_bot-dp);
1301  std::reverse_iterator<KernelIter2d> k_rbottom_right = std::reverse_iterator<KernelIter2d>(kernel_up);
1302 
1303  slip::Array2d<k_value_type> reverse_kernel(k_rows,k_cols);
1304  std::copy(k_rupper_left,k_rbottom_right,reverse_kernel.upper_left());
1305 
1306  // apply the convolution on a point where the kernel does
1307  // fit in the signal
1308  for(std::size_t i = 0, ii = ksize_bot; i < r_rows; ++i, ++ii)
1309  {
1310  for(std::size_t j = 0, jj = ksize_right; j < r_cols; ++j, ++jj)
1311  {
1312  for(std::size_t k = 0, kk = (ii - ksize_bot); k < k_rows; ++k, ++kk)
1313  {
1314  *(R_up.row_begin(i)+j)+=std::inner_product(S_up.row_begin(kk)+j,
1315  S_up.row_begin(kk)+j+k_cols,
1316  reverse_kernel.row_begin(k),value_type(0));
1317  }
1318  }
1319  }
1320  }
1321 
1322 
1368  template <typename SrcIter2d,
1369  typename KernelIter2d,
1370  typename ResIter2d>
1371  inline
1372  void same_convolution2d(SrcIter2d S_up,
1373  SrcIter2d S_bot,
1374  KernelIter2d kernel_up,
1375  KernelIter2d kernel_bot,
1376  int ksize_left,
1377  int ksize_right,
1378  int ksize_up,
1379  int ksize_bot,
1380  ResIter2d R_up,
1381  ResIter2d R_bot,
1382  slip::BORDER_TREATMENT bt = slip::BORDER_TREATMENT_ZERO_PADDED)
1383  {
1384 
1385  assert((R_bot - R_up)[0] == ((S_bot - S_up)[0]));
1386  assert((R_bot - R_up)[1] == ((S_bot - S_up)[1]));
1387  assert(static_cast<int>((kernel_bot - kernel_up)[1]) == static_cast<int>(ksize_left + ksize_right + 1));
1388  assert(static_cast<int>((kernel_bot - kernel_up)[0]) == static_cast<int>(ksize_up + ksize_bot + 1));
1389  int s_rows = static_cast<int>((S_bot - S_up)[0]);
1390  int s_cols = static_cast<int>((S_bot - S_up)[1]);
1391  int k_rows = static_cast<int>((kernel_bot - kernel_up)[0]);
1392  int k_cols = static_cast<int>((kernel_bot - kernel_up)[1]);
1393  int r_rows = static_cast<int>((R_bot - R_up)[0]);
1394  int r_cols = static_cast<int>((R_bot - R_up)[1]);
1395 
1396 
1397  typedef typename SrcIter2d::value_type value_type;
1398  typedef typename KernelIter2d::value_type k_value_type;
1399 
1400  //temporary array used to compute the reverse kernel
1401  slip::DPoint2d<int> dp(1,0);
1402  std::reverse_iterator<KernelIter2d> k_rupper_left = std::reverse_iterator<KernelIter2d>(kernel_bot-dp);
1403  std::reverse_iterator<KernelIter2d> k_rbottom_right = std::reverse_iterator<KernelIter2d>(kernel_up);
1404 
1405  slip::Array2d<k_value_type> reverse_kernel(k_rows,k_cols);
1406  std::copy(k_rupper_left,k_rbottom_right,reverse_kernel.upper_left());
1407 
1408  // apply the convolution on a point where the kernel does
1409  // fit in the signal
1410  for(int ii = ksize_bot; ii < (r_rows - ksize_up); ++ii)
1411  {
1412  for(int j = 0, jj = ksize_right; jj < (r_cols - ksize_left); ++j, ++jj)
1413  {
1414  value_type tmp = value_type();
1415  for(int k = 0, kk = (ii - ksize_bot); k < k_rows; ++k, ++kk)
1416  {
1417  tmp += std::inner_product(S_up.row_begin(kk)+j,
1418  S_up.row_begin(kk)+j+k_cols,
1419  reverse_kernel.row_begin(k),
1420  value_type());
1421  }
1422  *(R_up.row_begin(ii)+jj) = tmp;
1423  }
1424  }
1425 
1426 
1427  //border treatment
1428  switch(bt)
1429  {
1430  case slip::BORDER_TREATMENT_AVOID: break;
1431  case slip::BORDER_TREATMENT_NEUMANN:
1432  {
1433  //up
1434  for(int ii = 0; ii < ksize_bot; ++ii)
1435  {
1436  for(int j = 0,jj = ksize_right; jj < (r_cols - ksize_left); ++j, ++jj)
1437  {
1438  value_type tmp = value_type();
1439  for(int k = 0; k < (ksize_bot-ii); ++k)
1440  {
1441  tmp += std::inner_product(S_up.row_begin(0)+j,
1442  S_up.row_begin(0)+j+k_cols,
1443  reverse_kernel.row_begin(k),
1444  value_type());
1445  }
1446 
1447  for(int k = (ksize_bot-ii), kk = 0; k < k_rows; ++k, ++kk)
1448  {
1449 
1450  tmp += std::inner_product(S_up.row_begin(kk)+j,
1451  S_up.row_begin(kk)+j+k_cols,
1452  reverse_kernel.row_begin(k),
1453  value_type());
1454  }
1455 
1456  *(R_up.row_begin(ii)+jj) = tmp;
1457 
1458  }
1459  }
1460 
1461  //bottom
1462  for(int i = 0, ii = (r_rows - ksize_up); ii < r_rows; ++i, ++ii)
1463  {
1464  for(int j = 0, jj = ksize_right; jj < (r_cols - ksize_left); ++j, ++jj)
1465  {
1466  value_type tmp = value_type();
1467  for(int kk = (ii - ksize_bot), k = 0; kk < r_rows; ++k, ++kk)
1468  {
1469  tmp += std::inner_product(S_up.row_begin(kk)+j,
1470  S_up.row_begin(kk)+j+k_cols,
1471  reverse_kernel.row_begin(k),
1472  value_type());
1473  }
1474  for(int k = k_rows-1; k > (k_rows - 2 - i); --k)
1475 
1476  {
1477  tmp += std::inner_product(S_up.row_begin(r_rows-1)+j,
1478  S_up.row_begin(r_rows-1)+j+k_cols,
1479  reverse_kernel.row_begin(k),
1480  value_type());
1481 
1482  }
1483  *(R_up.row_begin(ii)+jj) = tmp;
1484 
1485  }
1486  }
1487  //right
1488  for(int ii = ksize_bot; ii < (r_rows - ksize_up); ++ii)
1489  {
1490  for(int j = (r_cols - ksize_left - ksize_right), jj = (r_cols - ksize_left),cpt=0; jj < r_cols; ++j, ++jj,++cpt)
1491  {
1492  value_type tmp = value_type();
1493  for(int k = 0, kk = (ii - ksize_bot); k < k_rows; ++k, ++kk)
1494  {
1495  tmp += std::inner_product(S_up.row_begin(kk)+j,
1496  S_up.row_end(kk),
1497  reverse_kernel.row_begin(k),
1498  value_type());
1499 
1500  }
1501  for(int c = (k_cols - 1); c >= (k_cols - 1-cpt); --c)
1502  {
1503  tmp += std::inner_product(reverse_kernel.col_begin(c),
1504  reverse_kernel.col_end(c),
1505  S_up.col_begin(s_cols-1)+(ii-ksize_bot),
1506  value_type());
1507  }
1508  *(R_up.row_begin(ii)+jj) = tmp;
1509  }
1510  }
1511  //left
1512  for(int ii = ksize_bot; ii < (r_rows - ksize_up); ++ii)
1513  {
1514  for(int j = ksize_right, jj = 0; jj < ksize_right; --j, ++jj)
1515  {
1516  value_type tmp = value_type();
1517  for(int k = 0, kk = (ii - ksize_bot); k < k_rows; ++k, ++kk)
1518  {
1519  tmp += std::inner_product(reverse_kernel.row_begin(k) + j,
1520  reverse_kernel.row_end(k),
1521  S_up.row_begin(kk),
1522  value_type());
1523  }
1524  for(int c = (j-1); c >= 0; --c)
1525  {
1526  tmp += std::inner_product(reverse_kernel.col_begin(c),
1527  reverse_kernel.col_end(c),
1528  S_up.col_begin(0)+(ii-ksize_bot),
1529  value_type());
1530  }
1531  *(R_up.row_begin(ii)+jj) = tmp;
1532  }
1533  }
1534  //upper_left
1535  //temporary array to extend S
1536  slip::Array2d<value_type> Stmp(k_rows+ksize_bot-1,k_cols+ksize_right-1);
1537  int stmp_rows = static_cast<int>(Stmp.rows());
1538  int stmp_cols = static_cast<int>(Stmp.cols());
1539  for(int i = ksize_bot, ii = 0 ; i < stmp_rows; ++i, ++ii)
1540  {
1541  std::copy(S_up.row_begin(ii),
1542  S_up.row_begin(ii) + (k_cols - 1),
1543  Stmp.row_begin(i) + ksize_right);
1544  }
1545  for(int i = 0; i < ksize_bot; ++i)
1546  {
1547  std::copy(S_up.row_begin(0),
1548  S_up.row_begin(0) + (k_cols - 1),
1549  Stmp.row_begin(i) + ksize_right);
1550  }
1551  for(int j = 0; j < ksize_right; ++j)
1552  {
1553  std::copy(Stmp.col_begin(ksize_right),
1554  Stmp.col_end(ksize_right),
1555  Stmp.col_begin(j));
1556  }
1557  //compute the convolution
1558  for(int ii = 0; ii < ksize_bot; ++ii)
1559  {
1560  for(int j = 0,jj = 0; jj < ksize_right; ++jj, ++j)
1561  {
1562  value_type tmp = value_type();
1563  for(int k = 0, kk = ii; k < k_rows; ++k, ++kk)
1564  {
1565  tmp += std::inner_product(reverse_kernel.row_begin(k),
1566  reverse_kernel.row_end(k),
1567  Stmp.row_begin(kk)+j,
1568  value_type());
1569  }
1570 
1571  *(R_up.row_begin(ii)+jj) = tmp;
1572  }
1573  }
1574 
1575  //bottom_right
1576  //temporary array to extend S
1577  Stmp.resize(ksize_up+k_rows-1,ksize_left+k_cols-1);
1578  stmp_rows = static_cast<int>(Stmp.rows());
1579  stmp_cols = static_cast<int>(Stmp.cols());
1580  int offsetj = (s_cols - k_cols + 1);
1581 
1582  for(int i = 0, ii = ((r_rows - 1) - (ksize_bot + ksize_up - 1)) ; i < (ksize_up + 1); ++i, ++ii)
1583  {
1584  std::copy(S_up.row_begin(ii) + offsetj,
1585  S_up.row_end(ii),
1586  Stmp.row_begin(i));
1587  }
1588  for(int i = (ksize_up + 1); i < stmp_rows; ++i)
1589  {
1590  std::copy(S_up.row_begin(s_rows - 1) + offsetj,
1591  S_up.row_end(s_rows - 1),
1592  Stmp.row_begin(i));
1593  }
1594 
1595  for(int j = stmp_cols-ksize_left; j < stmp_cols ; ++j)
1596  {
1597  std::copy(Stmp.col_begin(stmp_cols-1-ksize_left),
1598  Stmp.col_end(stmp_cols-1-ksize_left),
1599  Stmp.col_begin(j));
1600  }
1601 
1602  //computes the convolution
1603  for(int ii = (r_rows - ksize_up), i = 0; ii < r_rows; ++ii,++i)
1604  {
1605  for(int j = 0, jj = (r_cols - ksize_left); jj < r_cols; ++j, ++jj)
1606  {
1607  value_type tmp = value_type();
1608  for(int kk = i, k = 0; k < k_rows; ++k, ++kk)
1609  {
1610  tmp += std::inner_product(reverse_kernel.row_begin(k),
1611  reverse_kernel.row_end(k),
1612  Stmp.row_begin(kk)+j,
1613  value_type());
1614  }
1615  *(R_up.row_begin(ii)+jj) = tmp;
1616  }
1617  }
1618 
1619  //upper_right
1620  //temporary array to extend S
1621  Stmp.resize(k_rows + ksize_bot - 1,k_cols + ksize_left - 1);
1622  stmp_rows = static_cast<int>(Stmp.rows());
1623  stmp_cols = static_cast<int>(Stmp.cols());
1624  for(int i = ksize_bot, ii = 0 ; i < stmp_rows; ++i, ++ii)
1625  {
1626  std::copy(S_up.row_begin(ii)+ offsetj,
1627  S_up.row_end(ii),
1628  Stmp.row_begin(i));
1629  }
1630 
1631  for(int i = 0; i < ksize_bot; ++i)
1632  {
1633  std::copy(S_up.row_begin(0) + offsetj,
1634  S_up.row_end(0),
1635  Stmp.row_begin(i));
1636  }
1637  for(int j = (k_cols - 1); j < stmp_cols; ++j)
1638  {
1639  std::copy(Stmp.col_begin(k_cols - 2),
1640  Stmp.col_end(k_cols - 2),
1641  Stmp.col_begin(j));
1642  }
1643 
1644 
1645  //computes the convolution
1646  for(int ii = 0; ii < ksize_bot; ++ii)
1647  {
1648  for(int j = 0, jj = (r_cols - ksize_left); jj < r_cols;++j,++jj)
1649  {
1650  value_type tmp = value_type();
1651  for(int k = 0, kk = ii; k < k_rows; ++k, ++kk)
1652  {
1653  tmp += std::inner_product(reverse_kernel.row_begin(k),
1654  reverse_kernel.row_end(k),
1655  Stmp.row_begin(kk)+j,
1656  value_type());
1657 
1658  }
1659  *(R_up.row_begin(ii)+jj) = tmp;
1660  }
1661  }
1662 
1663  //bottom_left
1664  //temporary array to extend S
1665  Stmp.resize(k_rows + ksize_up - 1,k_cols + ksize_right - 1);
1666  stmp_rows = static_cast<int>(Stmp.rows());
1667  stmp_cols = static_cast<int>(Stmp.cols());
1668  offsetj = (s_cols - 1 - ksize_left);
1669  for(int i = 0, ii = ((r_rows - 1) - (ksize_bot + ksize_up - 1)) ; i < (ksize_up + 1); ++i, ++ii)
1670  {
1671  std::copy(S_up.row_begin(ii),
1672  S_up.row_begin(ii) + (k_cols - 1),
1673  Stmp.row_begin(i) + ksize_right);
1674  }
1675  for(int i = (ksize_up + 1); i < stmp_rows; ++i)
1676  {
1677  std::copy(S_up.row_begin(r_rows - 1),
1678  S_up.row_begin(r_rows - 1)+ (k_cols - 1),
1679  Stmp.row_begin(i)+ ksize_right);
1680  }
1681  for(int j = 0; j < ksize_right; ++j)
1682  {
1683  std::copy(Stmp.col_begin(ksize_right),
1684  Stmp.col_end(ksize_right),
1685  Stmp.col_begin(j));
1686  }
1687  //computes the convolution
1688  for(int ii = (r_rows - ksize_up), i = 0; ii < r_rows; ++ii, ++i)
1689  {
1690  for(int j = 0, jj = 0; jj < ksize_right; ++j, ++jj)
1691  {
1692  value_type tmp = value_type();
1693  for(int kk = i, k = 0; k < k_rows; ++k, ++kk)
1694  {
1695  tmp += std::inner_product(reverse_kernel.row_begin(k),
1696  reverse_kernel.row_end(k),
1697  Stmp.row_begin(kk)+j,
1698  value_type());
1699 
1700 
1701  }
1702  *(R_up.row_begin(ii)+jj) = tmp;
1703  }
1704  }
1705 
1706  }
1707  break;
1708  case slip::BORDER_TREATMENT_ZERO_PADDED:
1709  //up
1710  for(int ii = 0; ii < ksize_bot; ++ii)
1711  {
1712  for(int j = 0, jj = ksize_right; jj < (r_cols - ksize_left); ++j, ++jj)
1713  {
1714  value_type tmp = value_type();
1715  for(int k = (ksize_bot - ii), kk = 0; k < k_rows; ++k, ++kk)
1716  {
1717  tmp += std::inner_product(S_up.row_begin(kk)+j,
1718  S_up.row_begin(kk)+j+k_cols,
1719  reverse_kernel.row_begin(k),
1720  value_type());
1721  }
1722 
1723  *(R_up.row_begin(ii)+jj) = tmp;
1724 
1725  }
1726  }
1727  //bottom
1728  for(int ii = (r_rows - ksize_up); ii < r_rows; ++ii)
1729  {
1730  for(int j = 0, jj = ksize_right; jj < (r_cols - ksize_left); ++j, ++jj)
1731  {
1732 
1733  value_type tmp = value_type();
1734  for(int kk = (ii - ksize_bot), k = 0; kk < r_rows; ++k, ++kk)
1735  {
1736  tmp += std::inner_product(S_up.row_begin(kk)+j,
1737  S_up.row_begin(kk)+j+k_cols,
1738  reverse_kernel.row_begin(k),
1739  value_type());
1740  }
1741 
1742  *(R_up.row_begin(ii)+jj) = tmp;
1743 
1744  }
1745  }
1746 
1747  //right
1748  for(int ii = ksize_bot; ii < (r_rows - ksize_up); ++ii)
1749  {
1750  for(int j = (r_cols - ksize_left - ksize_right), jj = (r_cols - ksize_left); jj < r_cols; ++j, ++jj)
1751  {
1752  value_type tmp = value_type();
1753  for(int k = 0, kk = (ii - ksize_bot); k < k_rows; ++k, ++kk)
1754  {
1755  tmp += std::inner_product(S_up.row_begin(kk)+j,
1756  S_up.row_end(kk),
1757  reverse_kernel.row_begin(k),
1758  value_type());
1759 
1760  }
1761  *(R_up.row_begin(ii)+jj) = tmp;
1762  }
1763  }
1764 
1765  //left
1766  for(int ii = ksize_bot; ii < (r_rows - ksize_up); ++ii)
1767  {
1768  for(int j = ksize_right, jj = 0; jj < ksize_right; --j, ++jj)
1769  {
1770  value_type tmp = value_type();
1771  for(int k = 0, kk = (ii - ksize_bot); k < k_rows; ++k, ++kk)
1772  {
1773  tmp += std::inner_product(reverse_kernel.row_begin(k) + j,
1774  reverse_kernel.row_end(k),
1775  S_up.row_begin(kk),
1776  value_type());
1777  }
1778  *(R_up.row_begin(ii)+jj) = tmp;
1779  }
1780  }
1781  //upper_left
1782  for(int ii = 0; ii < ksize_bot; ++ii)
1783  {
1784  for(int j = ksize_right, jj = 0; jj < ksize_right; --j, ++jj)
1785  {
1786  value_type tmp = value_type();
1787  for(int k = (ksize_bot - ii), kk = 0; k < k_rows; ++k, ++kk)
1788  {
1789  tmp += std::inner_product(reverse_kernel.row_begin(k) + j,
1790  reverse_kernel.row_end(k),
1791  S_up.row_begin(kk),
1792  value_type());
1793  }
1794  *(R_up.row_begin(ii)+jj) = tmp;
1795  }
1796  }
1797 
1798  //upper_right
1799  for(int ii = 0; ii < ksize_bot; ++ii)
1800  {
1801  for(int j = (r_cols - ksize_left - ksize_right), jj = (r_cols - ksize_left); jj < r_cols; ++j, ++jj)
1802  {
1803  //std::cout<<"ii = "<<ii<<" jj = "<<jj<<std::endl;
1804  value_type tmp = value_type();
1805  for(int k = (ksize_bot - ii), kk = 0; k < k_rows; ++k, ++kk)
1806  {
1807  tmp += std::inner_product(S_up.row_begin(kk)+j,
1808  S_up.row_end(kk),
1809  reverse_kernel.row_begin(k),
1810  value_type());
1811  }
1812  *(R_up.row_begin(ii)+jj) = tmp;
1813  }
1814  }
1815 
1816  //bottom_right
1817  for(int ii = (r_rows - ksize_up); ii < r_rows; ++ii)
1818  {
1819  for(int j = (r_cols - ksize_left - ksize_right), jj = (r_cols - ksize_left); jj < r_cols; ++j, ++jj)
1820  {
1821 
1822  value_type tmp = value_type();
1823  for(int kk = (ii - ksize_bot), k = 0; kk < r_rows; ++k, ++kk)
1824  {
1825  tmp += std::inner_product(S_up.row_begin(kk)+j,
1826  S_up.row_end(kk),
1827  reverse_kernel.row_begin(k),
1828  value_type());
1829  }
1830  *(R_up.row_begin(ii)+jj) = tmp;
1831  }
1832  }
1833  // bottom_left
1834  for(int ii = (r_rows - ksize_up); ii < r_rows; ++ii)
1835  {
1836  for(int j = ksize_right, jj = 0; jj < ksize_right; --j, ++jj)
1837  {
1838  value_type tmp = value_type();
1839  for(int kk = (ii - ksize_bot), k = 0; kk < r_rows; ++k, ++kk)
1840  {
1841  tmp += std::inner_product(reverse_kernel.row_begin(k) + j,
1842  reverse_kernel.row_end(k),
1843  S_up.row_begin(kk),
1844  value_type());
1845  }
1846  *(R_up.row_begin(ii)+jj) = tmp;
1847  }
1848  }
1849  break;
1850  case slip::BORDER_TREATMENT_CIRCULAR:
1851  {
1852 
1853  //up
1854  for(unsigned int ii = 0; ii < static_cast<unsigned int>(ksize_bot); ++ii)
1855  {
1856  for(unsigned int j = 0,jj = static_cast<unsigned int>(ksize_right); jj < static_cast<unsigned int>(r_cols - ksize_left); ++j, ++jj)
1857  {
1858  value_type tmp = value_type();
1859  unsigned int jjj = j%s_cols;
1860  for(unsigned int k = 0, kk = (ii - static_cast<unsigned int>(ksize_bot)); k < static_cast<unsigned int>(k_rows); ++k, ++kk)
1861  {
1862  unsigned int kkk = kk%s_rows;
1863  tmp += std::inner_product(S_up.row_begin(kkk)+jjj,
1864  S_up.row_begin(kkk)+jjj+k_cols,
1865  reverse_kernel.row_begin(k),
1866  value_type());
1867  }
1868  *(R_up.row_begin(ii)+jj) = tmp;
1869  }
1870  }
1871  //bottom
1872  for(unsigned int ii = static_cast<unsigned int>(r_rows - ksize_up); ii < static_cast<unsigned int>(r_rows); ++ii)
1873  {
1874  for(unsigned int j = 0, jj = static_cast<unsigned int>(ksize_right); jj < static_cast<unsigned int>(r_cols - ksize_left); ++j, ++jj)
1875  {
1876  value_type tmp = value_type();
1877  unsigned int jjj = j%s_cols;
1878  for(unsigned int k = 0, kk = (ii - static_cast<unsigned int>(ksize_bot)); k < static_cast<unsigned int>(k_rows); ++k, ++kk)
1879  {
1880  unsigned int kkk = kk%s_rows;
1881  tmp += std::inner_product(S_up.row_begin(kkk)+jjj,
1882  S_up.row_begin(kkk)+jjj+k_cols,
1883  reverse_kernel.row_begin(k),
1884  value_type());
1885  }
1886  *(R_up.row_begin(ii)+jj) = tmp;
1887  }
1888  }
1889  //right
1890 for(unsigned int ii = static_cast<unsigned int>(ksize_bot); ii < static_cast<unsigned int>(r_rows - ksize_up); ++ii)
1891  {
1892  for(unsigned int j = static_cast<unsigned int>(r_cols - ksize_left - ksize_right), jj = static_cast<unsigned int>(r_cols - ksize_left); jj < static_cast<unsigned int>(r_cols); ++j, ++jj)
1893  {
1894  value_type tmp = value_type();
1895  slip::DPoint2d<int> d;
1896  for(unsigned int k = 0, kk = (ii - static_cast<unsigned int>(ksize_bot)); k < static_cast<unsigned int>(k_rows); ++k, ++kk)
1897  {
1898  d[0] = static_cast<int>(kk);
1899  for(unsigned int l = 0, ll = (jj-static_cast<unsigned int>(ksize_right)); l < static_cast<unsigned int>(k_cols); ++l,++ll)
1900  {
1901  d[1] = static_cast<int>(ll%s_cols);
1902  tmp += S_up[d] * reverse_kernel[k][l];
1903  }
1904  }
1905  *(R_up.row_begin(ii)+jj) = tmp;
1906  }
1907  }
1908  //left
1909 for(unsigned int ii = static_cast<unsigned int>(ksize_bot); ii < static_cast<unsigned int>(r_rows - ksize_up); ++ii)
1910  {
1911  for(unsigned int j = static_cast<unsigned int>(ksize_right), jj = 0; jj < static_cast<unsigned int>(ksize_right); --j, ++jj)
1912  {
1913  value_type tmp = value_type();
1914  slip::DPoint2d<int> d;
1915  for(unsigned int k = 0, kk = (ii - static_cast<unsigned int>(ksize_bot)); k < static_cast<unsigned int>(k_rows); ++k, ++kk)
1916  {
1917  d[0] = static_cast<int>(kk);
1918  for(unsigned int l = 0, ll = (s_cols - ksize_right+jj); l < static_cast<unsigned int>(k_cols); ++l,++ll)
1919  {
1920  d[1] = static_cast<int>(ll%s_cols);
1921  tmp += S_up[d] * reverse_kernel[k][l];
1922  }
1923  }
1924  *(R_up.row_begin(ii)+jj) = tmp;
1925  }
1926  }
1927  //upper_left
1928 for(int ii = 0; ii < ksize_bot; ++ii)
1929  {
1930  for(int j = ksize_right, jj = 0; jj < ksize_right; --j, ++jj)
1931  {
1932  value_type tmp = value_type();
1933  slip::DPoint2d<int> d;
1934  for(unsigned int k = 0, kk = static_cast<unsigned int>(s_rows-ksize_bot+ii); k < static_cast<unsigned int>(k_rows); ++k, ++kk)
1935  {
1936  d[0] = static_cast<int>(kk%s_rows);
1937  for(unsigned int l = 0, ll = static_cast<unsigned int>(s_cols - ksize_right+jj); l < static_cast<unsigned int>(k_cols); ++l,++ll)
1938  {
1939  d[1] = static_cast<int>(ll%s_cols);
1940  tmp += S_up[d] * reverse_kernel[k][l];
1941  }
1942  }
1943  *(R_up.row_begin(ii)+jj) = tmp;
1944  }
1945  }
1946 
1947  //upper_right
1948  for(int ii = 0; ii < ksize_bot; ++ii)
1949  {
1950  for(int j = (r_cols - ksize_left - ksize_right), jj = (r_cols - ksize_left); jj < r_cols; ++j, ++jj)
1951  {
1952  value_type tmp = value_type();
1953  slip::DPoint2d<int> d;
1954  for(int k = 0, kk = (s_rows-ksize_bot+ii); k < k_rows; ++k, ++kk)
1955  {
1956  d[0] = static_cast<int>(kk%s_rows);
1957  for(unsigned int l = 0, ll = static_cast<unsigned int>(jj-ksize_right); l < static_cast<unsigned int>(k_cols); ++l,++ll)
1958  {
1959  d[1] = static_cast<int>(ll%s_cols);
1960  tmp += S_up[d] * reverse_kernel[k][l];
1961 
1962  }
1963  }
1964  *(R_up.row_begin(ii)+jj) = tmp;
1965  }
1966  }
1967  //bottom_right
1968  for(int ii = (r_rows - ksize_up); ii < r_rows; ++ii)
1969  {
1970  for(int j = (r_cols - ksize_left - ksize_right), jj = (r_cols - ksize_left); jj < r_cols; ++j, ++jj)
1971  {
1972 
1973  value_type tmp = value_type();
1974  slip::DPoint2d<int> d;
1975  for(unsigned int kk = static_cast<unsigned int>(s_rows-ksize_bot+ii), k = 0; k < static_cast<unsigned int>(k_rows); ++k, ++kk)
1976  {
1977  d[0] = static_cast<int>(kk%s_rows);
1978  for(unsigned int l = 0, ll = static_cast<unsigned int>(jj-ksize_right); l < static_cast<unsigned int>(k_cols); ++l,++ll)
1979  {
1980  d[1] = static_cast<int>(ll%s_cols);
1981  tmp += S_up[d] * reverse_kernel[k][l];
1982 
1983  }
1984  }
1985  *(R_up.row_begin(ii)+jj) = tmp;
1986  }
1987  }
1988  // bottom_left
1989  for(int ii = (r_rows - ksize_up); ii < r_rows; ++ii)
1990  {
1991  for(int j = ksize_right, jj = 0; jj < ksize_right; --j, ++jj)
1992  {
1993  value_type tmp = value_type();
1994  slip::DPoint2d<int> d;
1995  for(unsigned int kk = static_cast<unsigned int>(s_rows-ksize_bot+ii), k = 0; k < static_cast<unsigned int>(k_rows); ++k, ++kk)
1996  {
1997  d[0] = static_cast<int>(kk%s_rows);
1998  for(unsigned int l = 0, ll = static_cast<unsigned int>(s_cols - ksize_right+jj); l < static_cast<unsigned int>(k_cols); ++l,++ll)
1999  {
2000  d[1] = static_cast<int>(ll%s_cols);
2001  tmp += S_up[d] * reverse_kernel[k][l];
2002 
2003  }
2004  }
2005  *(R_up.row_begin(ii)+jj) = tmp;
2006  }
2007  }
2008  }
2009  break;
2010  default: std::cout<<"default"<<std::endl;
2011  };
2012  }
2013 
2014 
2052  template <typename SrcIter2d,
2053  typename KernelIter2d,
2054  typename ResIter2d>
2055  inline
2056  void same_convolution2d(SrcIter2d S_up,
2057  SrcIter2d S_bot,
2058  KernelIter2d kernel_up,
2059  KernelIter2d kernel_bot,
2060  ResIter2d R_up,
2061  ResIter2d R_bot,
2062  slip::BORDER_TREATMENT bt = slip::BORDER_TREATMENT_ZERO_PADDED)
2063  {
2064  int k_rows = static_cast<int>((kernel_bot - kernel_up)[0] - 1);
2065  int k_cols = static_cast<int>((kernel_bot - kernel_up)[1] - 1);
2066  int k_rows_2 = k_rows / 2;
2067  int k_cols_2 = k_cols / 2;
2068  slip::same_convolution2d(S_up,S_bot,
2069  kernel_up,kernel_bot,
2070  k_cols - k_cols_2,k_cols_2,
2071  k_rows - k_rows_2,k_rows_2,
2072  R_up,R_bot,
2073  bt);
2074  }
2075 
2116  template <typename Container2d1,
2117  typename Kernel2d,
2118  typename Container2d2>
2119  inline
2120  void same_convolution2d(const Container2d1& Signal,
2121  const Kernel2d& Kernel,
2122  Container2d2& Result,
2123  slip::BORDER_TREATMENT bt)
2124  {
2125  slip::same_convolution2d(Signal.upper_left(),Signal.bottom_right(),
2126  Kernel.upper_left(),Kernel.bottom_right(),
2127  Result.upper_left(),Result.bottom_right(),
2128  bt);
2129  }
2161  template <typename SrcIter2d,
2162  typename KernelIter2d,
2163  typename ResIter2d>
2164  inline
2165  void full_convolution2d(SrcIter2d S_up,
2166  SrcIter2d S_bot,
2167  KernelIter2d kernel_up,
2168  KernelIter2d kernel_bot,
2169  ResIter2d R_up,
2170  ResIter2d R_bot)
2171  {
2172  assert( (R_bot - R_up)[0] == ((S_bot - S_up)[0] + (kernel_bot - kernel_up)[0] - 1));
2173  assert( (R_bot - R_up)[1] == ((S_bot - S_up)[1] + (kernel_bot - kernel_up)[1] - 1));
2174 
2175  int s_rows = static_cast<int>((S_bot - S_up)[0]);
2176  int s_cols = static_cast<int>((S_bot - S_up)[1]);
2177  int k_rows = static_cast<int>((kernel_bot - kernel_up)[0]);
2178  int k_cols = static_cast<int>((kernel_bot - kernel_up)[1]);
2179 
2180  typedef typename SrcIter2d::value_type value_type;
2181 
2182  //temporary array used to compute the reverse kernel
2183  slip::DPoint2d<int> dp(1,0);
2184  std::reverse_iterator<KernelIter2d> k_rupper_left = std::reverse_iterator<KernelIter2d>(kernel_bot-dp);
2185  std::reverse_iterator<KernelIter2d> k_rbottom_right = std::reverse_iterator<KernelIter2d>(kernel_up);
2186 
2187  slip::Array2d<typename KernelIter2d::value_type>
2188  reverse_kernel(k_rows,k_cols);
2189  std::copy(k_rupper_left,k_rbottom_right,reverse_kernel.upper_left());
2190 
2191 
2192  slip::DPoint2d<int> d(0,0);
2193  //up
2194  int imax = (k_rows - 1);
2195  int jmax = (s_cols - k_cols + 1);
2196  for(int i = 0, ii = 0; i < imax; ++i,++ii)
2197  {
2198  d[0] = ii;
2199  for(int j = 0, jj =(k_cols-1); j < jmax; ++j, ++jj)
2200  {
2201  d[1] = jj;
2202  value_type tmp = value_type();
2203  for(int k = (k_rows-1-i),kk = 0; k < k_rows; ++k,++kk)
2204  {
2205  tmp += std::inner_product(reverse_kernel.row_begin(k),
2206  reverse_kernel.row_end(k),
2207  S_up.row_begin(kk)+j,
2208  value_type(0));
2209  }
2210  R_up[d] = tmp;
2211  }
2212 
2213  }
2214 
2215 
2216  //center
2217  imax = (s_rows - k_rows+1);
2218  jmax = (s_cols - k_cols + 1);
2219  for(int i = 0, ii = (k_rows -1); i < imax; ++i,++ii)
2220  {
2221  d[0] = ii;
2222  for(int j = 0, jj =(k_cols-1); j < jmax; ++j, ++jj)
2223  {
2224  d[1] = jj;
2225  value_type tmp = value_type();
2226  for(int k = 0; k < k_rows; ++k)
2227  {
2228  tmp += std::inner_product(reverse_kernel.row_begin(k),
2229  reverse_kernel.row_end(k),
2230  S_up.row_begin(i+k)+j,
2231  value_type(0));
2232  }
2233  R_up[d] = tmp;
2234  }
2235 
2236  }
2237 
2238  //bottom
2239  jmax = (s_cols - k_cols + 1);
2240  for(int i = (s_rows - k_rows+1), ii = s_rows; i < s_rows; ++i,++ii)
2241  {
2242  d[0] = ii;
2243  for(int j = 0, jj =(k_cols-1); j < jmax; ++j, ++jj)
2244  {
2245  d[1] = jj;
2246  value_type tmp = value_type();
2247  for(int k = 0, kk = i; k < (s_rows-i); ++k,++kk)
2248  {
2249  tmp += std::inner_product(reverse_kernel.row_begin(k),
2250  reverse_kernel.row_end(k),
2251  S_up.row_begin(kk)+j,
2252  value_type(0));
2253  }
2254  R_up[d] = tmp;
2255  }
2256 
2257  }
2258 
2259  //left
2260  imax = (s_rows - k_rows + 1);
2261  jmax = k_cols - 1;
2262  for(int i = 0, ii = (k_rows -1); i < imax; ++i,++ii)
2263  {
2264  d[0] = ii;
2265  for(int j = 0, jj =0; j < jmax; ++j, ++jj)
2266  {
2267  d[1] = jj;
2268  value_type tmp = value_type();
2269  for(int k = 0; k < k_rows; ++k)
2270  {
2271  tmp += std::inner_product(reverse_kernel.row_begin(k)
2272  + (jmax - j),
2273  reverse_kernel.row_end(k),
2274  S_up.row_begin(i+k),
2275  value_type(0));
2276  }
2277  R_up[d] = tmp;
2278  }
2279 
2280  }
2281 
2282  //right
2283  imax = (s_rows - k_rows + 1);
2284  int jtmp = s_cols - k_cols;
2285  for(int i = 0, ii = (k_rows -1); i < imax; ++i,++ii)
2286  {
2287  d[0] = ii;
2288  for(int j = (s_cols - k_cols + 1), jj =s_cols; j < s_cols;++j, ++jj)
2289  {
2290  d[1] = jj;
2291  value_type tmp = value_type();
2292  for(int k = 0; k < k_rows; ++k)
2293  {
2294  tmp += std::inner_product(reverse_kernel.row_begin(k),
2295  reverse_kernel.row_end(k) - (j - jtmp),
2296  S_up.row_begin(i+k)+j,
2297  value_type(0));
2298  }
2299  R_up[d] = tmp;
2300  }
2301 
2302  }
2303 
2304  //upper_left
2305  imax = k_rows - 1;
2306  jmax = k_cols - 1;
2307  for(int i = 0, ii = 0; i < imax; ++i,++ii)
2308  {
2309  d[0] = ii;
2310  for(int j = 0, jj = 0; j < jmax; ++j, ++jj)
2311  {
2312  d[1] = jj;
2313  value_type tmp = value_type();
2314  for(int k = (k_rows-1-i),kk = 0; k < k_rows; ++k,++kk)
2315  {
2316  tmp += std::inner_product(reverse_kernel.row_begin(k)
2317  + (jmax - j),
2318  reverse_kernel.row_end(k),
2319  S_up.row_begin(kk),
2320  value_type(0));
2321 
2322  }
2323  R_up[d] = tmp;
2324  }
2325 
2326  }
2327 
2328  //upper_right
2329  imax = k_rows - 1;
2330  for(int i = 0, ii = 0; i < imax; ++i,++ii)
2331  {
2332  d[0] = ii;
2333  for(int j = (s_cols - k_cols + 1), jj = s_cols; j < s_cols; ++j, ++jj)
2334  {
2335  d[1] = jj;
2336  value_type tmp = value_type();
2337  for(int k = (k_rows-1-i),kk = 0; k < k_rows; ++k,++kk)
2338  {
2339  tmp += std::inner_product(reverse_kernel.row_begin(k),
2340  reverse_kernel.row_end(k)
2341  - (j - jtmp),
2342  S_up.row_begin(kk)+j,
2343  value_type(0));
2344  }
2345  R_up[d] = tmp;
2346  }
2347 
2348  }
2349 
2350  //bottom_right
2351  imax = (s_rows - k_rows + 1);
2352  jmax = (s_cols - k_cols + 1);
2353  for(int i = (s_rows - k_rows + 1), ii = s_rows; i < s_rows; ++i,++ii)
2354  {
2355  d[0] = ii;
2356  for(int j = (s_cols - k_cols + 1), jj = s_cols; j < s_cols; ++j, ++jj)
2357  {
2358  d[1] = jj;
2359  value_type tmp = value_type();
2360  for(int k = 0, kk = i; k < (s_rows-i); ++k,++kk)
2361  {
2362  tmp += std::inner_product(reverse_kernel.row_begin(k),
2363  reverse_kernel.row_end(k)- (j - jtmp),
2364  S_up.row_begin(kk)+j,
2365  value_type(0));
2366  }
2367  R_up[d] = tmp;
2368  }
2369 
2370  }
2371 
2372  //bottom_left
2373  imax = (s_rows - k_rows+1);
2374  jmax = k_cols - 1;
2375  for(int i = (s_rows - k_rows + 1), ii = s_rows; i < s_rows; ++i,++ii)
2376  {
2377  d[0] = ii;
2378  for(int j = 0, jj =0; j < jmax; ++j, ++jj)
2379  {
2380  d[1] = jj;
2381  value_type tmp = value_type();
2382  for(int k = 0, kk = i; k < (s_rows-i); ++k,++kk)
2383  {
2384  tmp += std::inner_product(reverse_kernel.row_begin(k)
2385  + (jmax - j),
2386  reverse_kernel.row_end(k),
2387  S_up.row_begin(kk),
2388  value_type(0));
2389  }
2390  R_up[d] = tmp;
2391  }
2392 
2393  }
2394 
2395 
2396  }
2397 
2430  template<typename InputIterator2d1, typename InputIterator2d2, typename OutputIterator2d>
2431  inline
2432  void fft_convolution2d(InputIterator2d1 in1_upper_left,
2433  InputIterator2d1 in1_bottom_right,
2434  InputIterator2d2 in2_upper_left,
2435  InputIterator2d2 in2_bottom_right,
2436  OutputIterator2d out_upper_left,
2437  OutputIterator2d out_bottom_right)
2438  {
2439  assert(((in1_bottom_right - in1_upper_left)[0]+(in2_bottom_right - in2_upper_left)[0] - 1) == (out_bottom_right - out_upper_left)[0]);
2440  assert(((in1_bottom_right - in1_upper_left)[1]+(in2_bottom_right - in2_upper_left)[1] - 1) == (out_bottom_right - out_upper_left)[1]);
2441 
2442  typename InputIterator2d1::difference_type size2din1 = in1_bottom_right - in1_upper_left;
2443  typename InputIterator2d2::difference_type size2din2 = in2_bottom_right - in2_upper_left;
2444  typename OutputIterator2d::difference_type size2dout = out_bottom_right - out_upper_left;
2445 
2446  typedef typename std::iterator_traits<InputIterator2d1>::value_type value_type;
2447  typedef typename slip::lin_alg_traits<value_type>::value_type Real;
2448 
2449  std::size_t rows = size2din1[0] + size2din2[0] - 1;
2450  std::size_t cols = size2din1[1] + size2din2[1] - 1;
2451  //zero-padd Signal and Kernel data
2452  slip::Array2d<Real> Signal_zero_padded(rows,cols,Real());
2453  slip::Array2d<Real> Kernel_zero_padded(rows,cols,Real());
2454  slip::Box2d<int> signal_box(0,0,
2455  static_cast<int>(size2din1[0]-1),
2456  static_cast<int>(size2din1[1]-1));
2457  slip::Box2d<int> kernel_box(0,0,
2458  static_cast<int>(size2din2[0]-1),
2459  static_cast<int>(size2din2[1]-1));
2460  std::copy(in1_upper_left,in1_bottom_right,Signal_zero_padded.upper_left(signal_box));
2461 
2462  std::copy(in2_upper_left,in2_bottom_right,
2463  Kernel_zero_padded.upper_left(kernel_box));
2464 
2465  slip::Array2d<std::complex<Real> > fft1(rows,cols,std::complex<Real>());
2466  slip::Array2d<std::complex<Real> > fft2(rows,cols,std::complex<Real>());
2467 
2468  //computes the two fft
2469  slip::real_fft2d(Signal_zero_padded,fft1);
2470  slip::real_fft2d(Kernel_zero_padded,fft2);
2471  //computes fft1*fft2
2472  slip::multiplies(fft1.begin(),fft1.end(),fft2.begin(),fft1.begin());
2473  //ifft2d
2474  slip::ifft2d(fft1.upper_left(), fft1.bottom_right(), fft2.upper_left());
2475 
2476  //the result should be real
2477  std::transform(fft2.begin(),fft2.end(),out_upper_left,slip::un_real<std::complex<Real>,Real>());
2478 
2479  }
2480 
2512 template <typename Container2d1,
2513  typename Kernel2d,
2514  typename Container2d2>
2515 inline void fft_convolution2d(const Container2d1& Signal,
2516  const Kernel2d& Kernel,
2517  Container2d2& Result)
2518 
2519 
2520  {
2521  slip::fft_convolution2d(Signal.upper_left(),Signal.bottom_right(),
2522  Kernel.upper_left(),Kernel.bottom_right(),
2523  Result.upper_left(),Result.bottom_right());
2524  }
2525 
2526  /* @} */
2527 
2528 
2530  /* @{ */
2531 
2590 template <typename RandomAccessIterator3d1,
2591  typename KernelIter1,
2592  typename KernelIter2,
2593  typename KernelIter3,
2594  typename RandomAccessIterator3d2>
2595 inline
2596 void separable_convolution3d(RandomAccessIterator3d1 I_up,
2597  RandomAccessIterator3d1 I_bot,
2598  KernelIter1 k1_first,
2599  KernelIter1 k1_last,
2600  std::size_t ksize_left1,
2601  std::size_t ksize_right1,
2602  slip::BORDER_TREATMENT bt1,
2603  KernelIter2 k2_first,
2604  KernelIter2 k2_last,
2605  std::size_t ksize_left2,
2606  std::size_t ksize_right2,
2607  slip::BORDER_TREATMENT bt2,
2608  KernelIter3 k3_first,
2609  KernelIter3 k3_last,
2610  std::size_t ksize_left3,
2611  std::size_t ksize_right3,
2612  slip::BORDER_TREATMENT bt3,
2613  RandomAccessIterator3d2 R_up)
2614 {
2615 
2616  typename RandomAccessIterator3d1::difference_type size3d =I_bot - I_up;
2617  typedef typename RandomAccessIterator3d1::size_type size_type;
2618  typedef typename RandomAccessIterator3d1::value_type value_type;
2619 
2620  size_type nb_slice = size3d[0];
2621  size_type nb_rows = size3d[1];
2622  size_type nb_cols = size3d[2];
2623 
2624 
2625  slip::Array<value_type> col_tmp(nb_cols);
2626  for(size_type k = 0; k < nb_slice; ++k)
2627  {
2628  for(size_type i = 0; i < nb_rows; ++i)
2629  {
2630  slip::same_convolution(I_up.row_begin(k,i),
2631  I_up.row_end(k,i),
2632  k1_first,
2633  k1_last,
2634  ksize_left1,
2635  ksize_right1,
2636  col_tmp.begin(),
2637  bt1);
2638  std::copy(col_tmp.begin(),col_tmp.end(),
2639  R_up.row_begin(k,i));
2640  }
2641  }
2642  slip::Array<value_type> row_tmp(nb_rows);
2643  for(size_type k = 0; k < nb_slice; ++k)
2644  {
2645  for(size_type j = 0; j < nb_cols; ++j)
2646  {
2647  slip::same_convolution(R_up.col_begin(k,j),
2648  R_up.col_end(k,j),
2649  k2_first,
2650  k2_last,
2651  ksize_left2,
2652  ksize_right2,
2653  row_tmp.begin(),
2654  bt2);
2655  std::copy(row_tmp.begin(),row_tmp.end(),
2656  R_up.col_begin(k,j));
2657  }
2658  }
2659 
2660  slip::Array<value_type> slice_tmp(nb_slice);
2661  for(size_type i = 0; i < nb_rows; ++i)
2662  {
2663  for(size_type j = 0; j < nb_cols; ++j)
2664  {
2665  slip::same_convolution(R_up.slice_begin(i,j),
2666  R_up.slice_end(i,j),
2667  k3_first,
2668  k3_last,
2669  ksize_left3,
2670  ksize_right3,
2671  slice_tmp.begin(),
2672  bt3);
2673  std::copy(slice_tmp.begin(),slice_tmp.end(),
2674  R_up.slice_begin(i,j));
2675  }
2676  }
2677 
2678 
2679 }
2680 
2681 
2722 template <typename RandomAccessIterator3d1,
2723  typename KernelIter1,
2724  typename KernelIter2,
2725  typename KernelIter3,
2726  typename RandomAccessIterator3d2>
2727 inline
2728 void separable_convolution3d(RandomAccessIterator3d1 I_up,
2729  RandomAccessIterator3d1 I_bot,
2730  KernelIter1 k1_first,
2731  KernelIter1 k1_last,
2732  KernelIter2 k2_first,
2733  KernelIter2 k2_last,
2734  KernelIter3 k3_first,
2735  KernelIter3 k3_last,
2736  RandomAccessIterator3d2 R_up,
2737  slip::BORDER_TREATMENT bt = slip::BORDER_TREATMENT_ZERO_PADDED)
2738 {
2739 
2740  std::size_t size1 = (k1_last-k1_first)/2;
2741  std::size_t size1_bis = (k1_last-k1_first) - 1 - size1;
2742  std::size_t size2 = (k2_last-k2_first)/2;
2743  std::size_t size2_bis = (k2_last-k2_first) - 1 - size2;
2744  std::size_t size3 = (k3_last-k3_first)/2;
2745  std::size_t size3_bis = (k3_last-k3_first) - 1 - size3;
2746 
2747 
2748  slip::separable_convolution3d(I_up,I_bot,
2749  k1_first,
2750  k1_last,
2751  size1,
2752  size1_bis,
2753  bt,
2754  k2_first,
2755  k2_last,
2756  size2,
2757  size2_bis,
2758  bt,
2759  k3_first,
2760  k3_last,
2761  size3,
2762  size3_bis,
2763  bt,
2764  R_up);
2765 
2766 }
2767 
2768 
2803  template <typename Container1,
2804  typename Kernel1,
2805  typename Kernel2,
2806  typename Kernel3,
2807  typename Container2>
2808  inline void separable_convolution3d(const Container1& Signal,
2809  const Kernel1& kernel1,
2810  const Kernel2& kernel2,
2811  const Kernel3& kernel3,
2812  Container2& Result,
2813  slip::BORDER_TREATMENT bt = slip::BORDER_TREATMENT_ZERO_PADDED)
2814 
2815 
2816  {
2817  slip::separable_convolution3d(Signal.front_upper_left(),
2818  Signal.back_bottom_right(),
2819  kernel1.begin(),kernel1.end(),
2820  kernel2.begin(),kernel2.end(),
2821  kernel3.begin(),kernel3.end(),
2822  Result.front_upper_left(),
2823  bt);
2824  }
2825 
2826 
2882 template <typename RandomAccessIterator3d1,
2883  typename KernelIter1,
2884  typename KernelIter2,
2885  typename KernelIter3>
2886 inline
2887 void separable_convolution3d(RandomAccessIterator3d1 I_up,
2888  RandomAccessIterator3d1 I_bot,
2889  KernelIter1 k1_first,
2890  KernelIter1 k1_last,
2891  std::size_t ksize_left1,
2892  std::size_t ksize_right1,
2893  slip::BORDER_TREATMENT bt1,
2894  KernelIter2 k2_first,
2895  KernelIter2 k2_last,
2896  std::size_t ksize_left2,
2897  std::size_t ksize_right2,
2898  slip::BORDER_TREATMENT bt2,
2899  KernelIter3 k3_first,
2900  KernelIter3 k3_last,
2901  std::size_t ksize_left3,
2902  std::size_t ksize_right3,
2903  slip::BORDER_TREATMENT bt3)
2904 {
2905 
2906  typename RandomAccessIterator3d1::difference_type size3d =I_bot - I_up;
2907  typedef typename RandomAccessIterator3d1::size_type size_type;
2908  typedef typename RandomAccessIterator3d1::value_type value_type;
2909 
2910  const size_type nb_slice = size3d[0];
2911  const size_type nb_rows = size3d[1];
2912  const size_type nb_cols = size3d[2];
2913 
2914 
2915  slip::Array<value_type> col_tmp(nb_cols);
2916  for(size_type k = 0; k < nb_slice; ++k)
2917  {
2918  for(size_type i = 0; i < nb_rows; ++i)
2919  {
2920  slip::same_convolution(I_up.row_begin(k,i),
2921  I_up.row_end(k,i),
2922  k1_first,
2923  k1_last,
2924  ksize_left1,
2925  ksize_right1,
2926  col_tmp.begin(),
2927  bt1);
2928  std::copy(col_tmp.begin(),col_tmp.end(),
2929  I_up.row_begin(k,i));
2930  }
2931  }
2932  slip::Array<value_type> row_tmp(nb_rows);
2933  for(size_type k = 0; k < nb_slice; ++k)
2934  {
2935  for(size_type j = 0; j < nb_cols; ++j)
2936  {
2937  slip::same_convolution(I_up.col_begin(k,j),
2938  I_up.col_end(k,j),
2939  k2_first,
2940  k2_last,
2941  ksize_left2,
2942  ksize_right2,
2943  row_tmp.begin(),
2944  bt2);
2945  std::copy(row_tmp.begin(),row_tmp.end(),
2946  I_up.col_begin(k,j));
2947  }
2948  }
2949 
2950  slip::Array<value_type> slice_tmp(nb_slice);
2951  for(size_type i = 0; i < nb_rows; ++i)
2952  {
2953  for(size_type j = 0; j < nb_cols; ++j)
2954  {
2955  slip::same_convolution(I_up.slice_begin(i,j),
2956  I_up.slice_end(i,j),
2957  k3_first,
2958  k3_last,
2959  ksize_left3,
2960  ksize_right3,
2961  slice_tmp.begin(),
2962  bt3);
2963  std::copy(slice_tmp.begin(),slice_tmp.end(),
2964  I_up.slice_begin(i,j));
2965  }
2966  }
2967 
2968 
2969 }
2970 
3007 template <typename RandomAccessIterator3d1,
3008  typename KernelIter1,
3009  typename KernelIter2,
3010  typename KernelIter3,
3011  typename RandomAccessIterator3d2>
3012 inline
3013 void separable_fft_convolution3d(RandomAccessIterator3d1 I_up,
3014  RandomAccessIterator3d1 I_bot,
3015  KernelIter1 k1_first,
3016  KernelIter1 k1_last,
3017  KernelIter2 k2_first,
3018  KernelIter2 k2_last,
3019  KernelIter3 k3_first,
3020  KernelIter3 k3_last,
3021  RandomAccessIterator3d2 R_up)
3022 {
3023  typename RandomAccessIterator3d1::difference_type size3d =
3024  I_bot - I_up;
3025  typedef typename RandomAccessIterator3d1::size_type size_type;
3026  typedef typename RandomAccessIterator3d1::value_type value_type;
3027  size_type nb_slices = size3d[0];
3028  size_type nb_rows = size3d[1];
3029  size_type nb_cols = size3d[2];
3030 
3031  size_type size1 = (k1_last - k1_first);
3032  size_type size2 = (k2_last - k2_first);
3033  size_type size3 = (k3_last - k3_first);
3034 
3035  size_type size2o2 = size2/2;
3036  size_type size1o2 = size1/2;
3037 
3038  slip::Array3d<value_type> Atmp1(nb_slices + size3 - 1,
3039  nb_rows + size1 - 1,
3040  nb_cols + size2 - 1);
3041  slip::Array3d<value_type> Atmp2(nb_slices + size3 - 1,
3042  nb_rows + size1 - 1,
3043  nb_cols + size2 - 1);
3044  slip::Array<value_type> Vtmp(Atmp2.slices());
3045 
3046 
3047  for(size_t k = 0; k < nb_slices; ++k)
3048  {
3049  for(size_t i = 0; i < nb_rows; ++i)
3050  {
3051  slip::fft_convolution(I_up.row_begin(k,i),
3052  I_up.row_end(k,i),
3053  k1_first,k1_last,
3054  Atmp1.row_begin(k,i),
3055  Atmp1.row_end(k,i));
3056  }
3057  }
3058  for(size_t k = 0; k < nb_slices; ++k)
3059  {
3060  for(size_t j = 0; j < nb_cols; ++j)
3061  {
3062  slip::fft_convolution(Atmp1.col_begin(k,j) + size2o2,
3063  Atmp1.col_end(k,j) - size2o2,
3064  k2_first,k2_last,
3065  Atmp2.col_begin(k,j),
3066  Atmp2.col_end(k,j));
3067  }
3068  }
3069  Atmp1.resize(0,0,0);
3070  for(size_t i = 0; i < nb_rows; ++i)
3071  {
3072  for(size_t j = 0; j < nb_cols; ++j)
3073  {
3074  slip::fft_convolution(Atmp2.slice_begin(i,j) + size1o2,
3075  Atmp2.slice_end(i,j) - size1o2,
3076  k3_first,k3_last,
3077  Vtmp.begin(),Vtmp.end());
3078  std::copy(Vtmp.begin()+size2o2,Vtmp.end()-size2o2,R_up.col_begin(i,j));
3079  }
3080  }
3081 
3082 
3083 }
3115  template <typename Container1,
3116  typename Kernel1,
3117  typename Kernel2,
3118  typename Kernel3,
3119  typename Container2>
3120  inline void separable_fft_convolution3d(const Container1& Signal,
3121  const Kernel1& kernel1,
3122  const Kernel2& kernel2,
3123  const Kernel3& kernel3,
3124  Container2& Result)
3125 
3126 
3127  {
3128  slip::separable_fft_convolution3d(Signal.front_upper_left(),
3129  Signal.back_bottom_right(),
3130  kernel1.begin(),kernel1.end(),
3131  kernel2.begin(),kernel2.end(),
3132  kernel3.begin(),kernel3.end(),
3133  Result.front_upper_left());
3134  }
3135 
3136 
3171 template<typename InputIterator3d1,
3172  typename InputIterator3d2,
3173  typename OutputIterator3d>
3174 inline
3175 void
3176 fft_convolution3d(InputIterator3d1 in1_front_upper_left,
3177  InputIterator3d1 in1_back_bottom_right,
3178  InputIterator3d2 in2_front_upper_left,
3179  InputIterator3d2 in2_back_bottom_right,
3180  OutputIterator3d out_front_upper_left,
3181  OutputIterator3d out_back_bottom_right)
3182  {
3183  assert(((in1_back_bottom_right - in1_front_upper_left)[0]+(in2_back_bottom_right - in2_front_upper_left)[0] - 1) == (out_back_bottom_right - out_front_upper_left)[0]);
3184  assert(((in1_back_bottom_right - in1_front_upper_left)[1]+(in2_back_bottom_right - in2_front_upper_left)[1] - 1) == (out_back_bottom_right - out_front_upper_left)[1]);
3185 assert(((in1_back_bottom_right - in1_front_upper_left)[2]+(in2_back_bottom_right - in2_front_upper_left)[2] - 1) == (out_back_bottom_right - out_front_upper_left)[2]);
3186 
3187  typename InputIterator3d1::difference_type size3din1 = in1_back_bottom_right - in1_front_upper_left;
3188  typename InputIterator3d2::difference_type size3din2 = in2_back_bottom_right - in2_front_upper_left;
3189  typename OutputIterator3d::difference_type size3dout = out_back_bottom_right - out_front_upper_left;
3190 
3191  typedef typename slip::lin_alg_traits<typename std::iterator_traits<InputIterator3d1>::value_type>::value_type Real;
3192 
3193  std::size_t slices = size3din1[0] + size3din2[0] - 1;
3194  std::size_t rows = size3din1[1] + size3din2[1] - 1;
3195  std::size_t cols = size3din1[2] + size3din2[2] - 1;
3196 
3197  //zero-padd Signal and Kernel data
3198  slip::Array3d<Real> Signal_zero_padded(slices,rows,cols,Real());
3199  slip::Array3d<Real> Kernel_zero_padded(slices,rows,cols,Real());
3200  slip::Box3d<int> signal_box(0,0,0,
3201  static_cast<int>(size3din1[0]-1),
3202  static_cast<int>(size3din1[1]-1),
3203  static_cast<int>(size3din1[2]-1));
3204  slip::Box3d<int> kernel_box(0,0,0,
3205  static_cast<int>(size3din2[0]-1),
3206  static_cast<int>(size3din2[1]-1),
3207  static_cast<int>(size3din2[2]-1));
3208  std::copy(in1_front_upper_left,
3209  in1_back_bottom_right,
3210  Signal_zero_padded.front_upper_left(signal_box));
3211  std::copy(in2_front_upper_left,in2_back_bottom_right,
3212  Kernel_zero_padded.front_upper_left(kernel_box));
3213 
3214 
3215  slip::Array3d<std::complex<Real> > fft1(slices,rows,cols,std::complex<Real>());
3216  slip::Array3d<std::complex<Real> > fft2(slices,rows,cols,std::complex<Real>());
3217  //computes the two fft
3218  slip::real_fft3d(Signal_zero_padded.front_upper_left(),
3219  Signal_zero_padded.back_bottom_right(),
3220  fft1.front_upper_left());
3221  slip::real_fft3d(Kernel_zero_padded.front_upper_left(),
3222  Kernel_zero_padded.back_bottom_right(),
3223  fft2.front_upper_left());
3224  //computes fft1*fft2
3225  slip::multiplies(fft1.begin(),fft1.end(),fft2.begin(),fft1.begin());
3226 
3227  //ifft3d
3228  slip::ifft3d(fft1.front_upper_left(), fft1.back_bottom_right(),
3229  fft2.front_upper_left());
3230  //the result should be real
3231  std::transform(fft2.front_upper_left(),fft2.back_bottom_right(),
3232  out_front_upper_left,
3233  slip::un_real<std::complex<Real>,Real>());
3234  }
3235 
3264  template <typename Container3d1,
3265  typename Kernel3d,
3266  typename Container3d2>
3267  inline void fft_convolution3d(const Container3d1& Signal,
3268  const Kernel3d& Kernel,
3269  Container3d2& Result)
3270 
3271 
3272  {
3273  slip::fft_convolution3d(Signal.front_upper_left(),
3274  Signal.back_bottom_right(),
3275  Kernel.front_upper_left(),
3276  Kernel.back_bottom_right(),
3277  Result.front_upper_left(),
3278  Result.back_bottom_right());
3279  }
3280 
3281 
3302  template<typename InputIterator3d1,
3303  typename InputIterator3d2,
3304  typename OutputIterator3d>
3305  inline
3306  void
3307  fft_convolution3d_same(InputIterator3d1 in1_front_upper_left,
3308  InputIterator3d1 in1_back_bottom_right,
3309  InputIterator3d2 in2_front_upper_left,
3310  InputIterator3d2 in2_back_bottom_right,
3311  OutputIterator3d out_front_upper_left,
3312  OutputIterator3d out_back_bottom_right)
3313  {
3314  assert( (in1_back_bottom_right - in1_front_upper_left)
3315  == (out_back_bottom_right - out_front_upper_left));
3316 
3317  typename InputIterator3d1::difference_type size3din1 = in1_back_bottom_right - in1_front_upper_left;
3318  typename InputIterator3d2::difference_type size3din2 = in2_back_bottom_right - in2_front_upper_left;
3319  typename OutputIterator3d::difference_type size3dout = out_back_bottom_right - out_front_upper_left;
3320 
3321  typedef typename slip::lin_alg_traits<typename std::iterator_traits<InputIterator3d1>::value_type>::value_type Real;
3322 
3323  std::size_t slices = size3din1[0] + size3din2[0] - 1;
3324  std::size_t rows = size3din1[1] + size3din2[1] - 1;
3325  std::size_t cols = size3din1[2] + size3din2[2] - 1;
3326 
3327  //zero-padd Signal and Kernel data
3328  slip::Array3d<Real> Signal_zero_padded(slices,rows,cols,Real());
3329  slip::Array3d<Real> Kernel_zero_padded(slices,rows,cols,Real());
3330  slip::Box3d<int> signal_box(0,0,0,
3331  static_cast<int>(size3din1[0]-1),
3332  static_cast<int>(size3din1[1]-1),
3333  static_cast<int>(size3din1[2]-1));
3334  slip::Box3d<int> kernel_box(0,0,0,
3335  static_cast<int>(size3din2[0]-1),
3336  static_cast<int>(size3din2[1]-1),
3337  static_cast<int>(size3din2[2]-1));
3338  std::copy(in1_front_upper_left,
3339  in1_back_bottom_right,
3340  Signal_zero_padded.front_upper_left(signal_box));
3341  std::copy(in2_front_upper_left,in2_back_bottom_right,
3342  Kernel_zero_padded.front_upper_left(kernel_box));
3343 
3344 
3345  slip::Array3d<std::complex<Real> > fft1(slices,rows,cols,std::complex<Real>());
3346  slip::Array3d<std::complex<Real> > fft2(slices,rows,cols,std::complex<Real>());
3347  //computes the two fft
3348  slip::real_fft3d(Signal_zero_padded.front_upper_left(),
3349  Signal_zero_padded.back_bottom_right(),
3350  fft1.front_upper_left());
3351  slip::real_fft3d(Kernel_zero_padded.front_upper_left(),
3352  Kernel_zero_padded.back_bottom_right(),
3353  fft2.front_upper_left());
3354  //computes fft1*fft2
3355  slip::multiplies(fft1.begin(),fft1.end(),fft2.begin(),fft1.begin());
3356 
3357  //ifft3d
3358  slip::ifft3d(fft1.front_upper_left(), fft1.back_bottom_right(),
3359  fft2.front_upper_left());
3360  //fftshift3d
3361  //slip::fftshift3d(fft2.front_upper_left(),fft2.back_bottom_right());
3362  //the result should be real
3363  slip::Box3d<int> out_box(size3din2[0]/2,size3din2[1]/2,size3din2[2]/2,
3364  static_cast<int>(size3din2[0]/2 + size3din1[0] - 1),
3365  static_cast<int>(size3din2[1]/2 + size3din1[1] - 1),
3366  static_cast<int>(size3din2[2]/2 + size3din1[2] - 1));
3367 
3368  std::transform(fft2.front_upper_left(out_box),
3369  fft2.back_bottom_right(out_box),
3370  out_front_upper_left,slip::un_real<std::complex<Real>,Real>());
3371 
3372  }
3373 
3430  template <typename SrcIter3d,
3431  typename KernelIter3d,
3432  typename ResIter3d>
3433  inline
3434  void valid_convolution3d(SrcIter3d S_up,
3435  SrcIter3d S_bot,
3436  KernelIter3d kernel_up,
3437  KernelIter3d kernel_bot,
3438  std::size_t ksize_front,
3439  std::size_t ksize_back,
3440  std::size_t ksize_left,
3441  std::size_t ksize_right,
3442  std::size_t ksize_up,
3443  std::size_t ksize_bot,
3444  ResIter3d R_up,
3445  ResIter3d R_bot)
3446  {
3447 
3448  assert((R_bot - R_up)[0] == ((S_bot - S_up)[0] - (kernel_bot - kernel_up)[0] + 1));
3449  assert((R_bot - R_up)[1] == ((S_bot - S_up)[1] - (kernel_bot - kernel_up)[1] + 1));
3450  assert((R_bot - R_up)[2] == ((S_bot - S_up)[2] - (kernel_bot - kernel_up)[2] + 1));
3451 
3452  assert(static_cast<int>((kernel_bot - kernel_up)[0]) == static_cast<int>(ksize_front + ksize_back + 1));
3453  assert(static_cast<int>((kernel_bot - kernel_up)[1]) == static_cast<int>(ksize_up + ksize_bot + 1));
3454  assert(static_cast<int>((kernel_bot - kernel_up)[2]) == static_cast<int>(ksize_left + ksize_right + 1));
3455 
3456 
3457  std::size_t k_slices = static_cast<std::size_t>((kernel_bot - kernel_up)[0]);
3458  std::size_t k_rows = static_cast<std::size_t>((kernel_bot - kernel_up)[1]);
3459  std::size_t k_cols = static_cast<std::size_t>((kernel_bot - kernel_up)[2]);
3460  std::size_t r_slices = static_cast<std::size_t>((R_bot - R_up)[0]);
3461  std::size_t r_rows = static_cast<std::size_t>((R_bot - R_up)[1]);
3462  std::size_t r_cols = static_cast<std::size_t>((R_bot - R_up)[2]);
3463 
3464 
3465  typedef typename SrcIter3d::value_type value_type;
3466  typedef typename KernelIter3d::value_type k_value_type;
3467 
3468  //temporary array used to compute the reverse kernel
3469  slip::DPoint3d<int> dp(1,1,0);
3470  std::reverse_iterator<KernelIter3d> k_rupper_left = std::reverse_iterator<KernelIter3d>(kernel_bot-dp);
3471  std::reverse_iterator<KernelIter3d> k_rbottom_right = std::reverse_iterator<KernelIter3d>(kernel_up);
3472 
3473  slip::Array3d<k_value_type> reverse_kernel(k_slices,k_rows,k_cols);
3474  std::copy(k_rupper_left,k_rbottom_right,reverse_kernel.front_upper_left());
3475 
3476  // apply the convolution on a point where the kernel does
3477  // fit in the signal
3478  for(std::size_t k = 0, kk = ksize_back; k < r_slices; ++k, ++kk)
3479  {
3480  for(std::size_t i = 0, ii = ksize_bot; i < r_rows; ++i, ++ii)
3481  {
3482  for(std::size_t j = 0, jj = ksize_right; j < r_cols; ++j, ++jj)
3483  {
3484  for(std::size_t s = 0, ss = (kk - ksize_back); s < k_slices; ++s, ++ss)
3485  {
3486  for(std::size_t l = 0, ll = (ii - ksize_bot); l < k_rows; ++l, ++ll)
3487  {
3488  *(R_up.row_begin(k,i)+j)+=
3489  std::inner_product(S_up.row_begin(ss,ll)+j,
3490  S_up.row_begin(ss,ll)+j+k_cols,
3491  reverse_kernel.row_begin(s,l),
3492  value_type(0));
3493  }
3494  }
3495  }
3496  }
3497  }
3498  }
3499 
3500  /* @} */
3501 }//slip::
3502 
3503 
3504 
3505 #endif //SLIP_CONVOLUTION_HPP
slip::valid_convolution3d
void valid_convolution3d(SrcIter3d S_up, SrcIter3d S_bot, KernelIter3d kernel_up, KernelIter3d kernel_bot, std::size_t ksize_front, std::size_t ksize_back, std::size_t ksize_left, std::size_t ksize_right, std::size_t ksize_up, std::size_t ksize_bot, ResIter3d R_up, ResIter3d R_bot)
Computes the valid convolution of 3d signal by a 3d-kernel.
Definition: convolution.hpp:3434
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
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::fft_circular_convolution
void fft_circular_convolution(RandomAccessIterator1 first, RandomAccessIterator1 last, RandomAccessIterator2 first2, RandomAccessIterator2 last2, RandomAccessIterator3 conv_first, RandomAccessIterator3 conv_last)
Computes the FFT circular convolution.
Definition: convolution.hpp:652
slip::separable_convolution2d
void separable_convolution2d(RandomAccessIterator2d1 I_up, RandomAccessIterator2d1 I_bot, KernelIter1 k1_first, KernelIter1 k1_last, std::size_t ksize_left1, std::size_t ksize_right1, slip::BORDER_TREATMENT bt1, KernelIter2 k2_first, KernelIter2 k2_last, std::size_t ksize_left2, std::size_t ksize_right2, slip::BORDER_TREATMENT bt2, RandomAccessIterator2d2 R_up)
Computes the separable convolution of a 2d signal by two asymmetric 1d-kernels support.
Definition: convolution.hpp:842
slip::lin_alg_traits::value_type
_Tp value_type
Definition: linear_algebra_traits.hpp:188
slip::DPoint3d
Difference of Point3D class, specialization of DPoint<CoordType,DIM> with DIM = 3.
Definition: DenseVector3dField3d.hpp:122
slip::Array3d::slice_end
slice_iterator slice_end(const size_type row, const size_type col)
Returns a read/write iterator that points to the one past the end element of the line (row...
Box2d.hpp
Provides a class to manipulate 2d box.
slip::Vector::end
iterator end()
Returns a read/write iterator that points one past the last element in the Vector. Iteration is done in ordinary element order.
Definition: Vector.hpp:1481
slip::Array3d
This is a three-dimensional dynamic and generic container. This container statisfies the Bidirectionn...
Definition: Array3d.hpp:109
slip::Array3d::begin
const_iterator begin() const
Returns a read-only (constant) iterator that points to the first element in the Array3d. Iteration is done in ordinary element order.
Definition: Array3d.hpp:3644
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::un_real
Real functor. Return the real part of x.
Definition: macros.hpp:98
slip::Array3d::slices
size_type slices() const
Returns the number of slices (first dimension size) in the Array3d.
Definition: Array3d.hpp:5196
slip::fft_convolution3d
void fft_convolution3d(InputIterator3d1 in1_front_upper_left, InputIterator3d1 in1_back_bottom_right, InputIterator3d2 in2_front_upper_left, InputIterator3d2 in2_back_bottom_right, OutputIterator3d out_front_upper_left, OutputIterator3d out_back_bottom_right)
Computes the convolution of two 3D sequences using fft3d.
Definition: convolution.hpp:3176
slip::fft_convolution3d_same
void fft_convolution3d_same(InputIterator3d1 in1_front_upper_left, InputIterator3d1 in1_back_bottom_right, InputIterator3d2 in2_front_upper_left, InputIterator3d2 in2_back_bottom_right, OutputIterator3d out_front_upper_left, OutputIterator3d out_back_bottom_right)
Computes the convolution between two 3D sequences using fft3d.
Definition: convolution.hpp:3307
slip::Array3d::row_begin
row_iterator row_begin(const size_type slice, const size_type row)
Returns a read/write iterator that points to the first element of the row row of the slice slice in t...
slip::Array3d::row_end
row_iterator row_end(const size_type slice, const size_type row)
Returns a read/write iterator that points to the past-the-end element of the row row of the slice sli...
slip::real_fft2d
void real_fft2d(InputBidirectionalIterator2d in_upper_left, InputBidirectionalIterator2d in_bottom_right, OutputBidirectionalIterator2d out_upper_left)
Computes the real fft2d of a container with 2d iterators.
Definition: FFT.hpp:2553
slip::fft_convolution2d
void fft_convolution2d(InputIterator2d1 in1_upper_left, InputIterator2d1 in1_bottom_right, InputIterator2d2 in2_upper_left, InputIterator2d2 in2_bottom_right, OutputIterator2d out_upper_left, OutputIterator2d out_bottom_right)
Computes the standard crosscorrelation between two 2D sequences using fft2d.
Definition: convolution.hpp:2432
Array.hpp
Provides a class to manipulate 1d dynamic and generic arrays.
slip::ifft3d
void ifft3d(InputBidirectionalIterator3d in_front_upper_left, InputBidirectionalIterator3d in_back_bottom_right, OutputBidirectionalIterator3d out_front_upper_left)
Computes the ifft3d of a container with 3d iterators.
Definition: FFT.hpp:2852
slip::same_convolution2d
void same_convolution2d(SrcIter2d S_up, SrcIter2d S_bot, KernelIter2d kernel_up, KernelIter2d kernel_bot, int ksize_left, int ksize_right, int ksize_up, int ksize_bot, ResIter2d R_up, ResIter2d R_bot, slip::BORDER_TREATMENT bt=slip::BORDER_TREATMENT_ZERO_PADDED)
Computes the same convolution of 2d signal by a 2d-kernel.
Definition: convolution.hpp:1372
border_treatment.hpp
Provides some algorithms to handle the border conditions of ranges.
slip::ifft
void ifft(InputIter in_begin, InputIter in_end, OutputIter out_begin)
Computes the inverse fft of a container.
Definition: FFT.hpp:2507
slip::valid_convolution2d
void valid_convolution2d(SrcIter2d S_up, SrcIter2d S_bot, KernelIter2d kernel_up, KernelIter2d kernel_bot, std::size_t ksize_left, std::size_t ksize_right, std::size_t ksize_up, std::size_t ksize_bot, ResIter2d R_up, ResIter2d R_bot)
Computes the valid convolution of 2d signal by a 2d-kernel.
Definition: convolution.hpp:1273
slip::Array2d::resize
void resize(const size_type d1, const size_type d2, const T &val=T())
Resizes a Array2d.
Definition: Array2d.hpp:2323
slip::separable_fft_convolution3d
void separable_fft_convolution3d(RandomAccessIterator3d1 I_up, RandomAccessIterator3d1 I_bot, KernelIter1 k1_first, KernelIter1 k1_last, KernelIter2 k2_first, KernelIter2 k2_last, KernelIter3 k3_first, KernelIter3 k3_last, RandomAccessIterator3d2 R_up)
Computes the separable 3d convolution of signal using the fft algorithm.
Definition: convolution.hpp:3013
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::separable_fft_convolution2d
void separable_fft_convolution2d(RandomAccessIterator2d1 I_up, RandomAccessIterator2d1 I_bot, KernelIter1 k1_first, KernelIter1 k1_last, KernelIter2 k2_first, KernelIter2 k2_last, RandomAccessIterator2d2 R_up)
Computes the separable 2d convolution of signal using the fft algorithm.
Definition: convolution.hpp:1142
slip::multiplies
void multiplies(InputIterator1 __first1, InputIterator1 __last1, InputIterator2 __first2, OutputIterator __result)
Computes the pointwise product of two ranges.
Definition: arithmetic_op.hpp:508
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
slip::Vector
Numerical Vector class. This container statisfies the RandomAccessContainer concepts of the Standard ...
Definition: Vector.hpp:104
slip::Array2d::rows
size_type rows() const
Returns the number of rows (first dimension size) in the Array2d.
Definition: Array2d.hpp:3139
slip::Array2d::col_begin
col_iterator col_begin(const size_type col)
Returns a read/write iterator that points to the first element of the column column in the Array2d...
slip::Array3d::col_end
col_iterator col_end(const size_type slice, const size_type col)
Returns a read/write iterator that points to the past-the-end element of the column column of the sli...
slip::Array2d::end
iterator end()
Returns a read/write iterator that points one past the last element in the Array2d. Iteration is done in ordinary element order.
Definition: Array2d.hpp:2352
slip::fft_convolution
void fft_convolution(RandomAccessIterator1 first, RandomAccessIterator1 last, RandomAccessIterator2 first2, RandomAccessIterator2 last2, RandomAccessIterator3 conv_first, RandomAccessIterator3 conv_last)
Computes the FFT convolution.
Definition: convolution.hpp:710
FFT.hpp
Provides some FFT algorithms.
slip::Array3d::col_begin
col_iterator col_begin(const size_type slice, const size_type col)
Returns a read/write iterator that points to the first element of the column column of the slice slic...
arithmetic_op.hpp
Provides some algorithms to computes arithmetical operations on ranges.
slip::Array2d::row_begin
row_iterator row_begin(const size_type row)
Returns a read/write iterator that points to the first element of the row row in the Array2d...
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::same_convolution
void same_convolution(SrcIter first, SrcIter last, KernelIter kernel_first, KernelIter kernel_last, std::size_t ksize_left, std::size_t ksize_right, ResIter result, slip::BORDER_TREATMENT bt=slip::BORDER_TREATMENT_ZERO_PADDED)
Computes the same convolution of signal by a 1d-kernel.
Definition: convolution.hpp:299
slip::valid_convolution
void valid_convolution(SrcIter first, SrcIter last, KernelIter kernel_first, KernelIter kernel_last, std::size_t ksize_left, std::size_t ksize_right, ResIter result)
Computes the valid convolution of signal by a 1d-kernel.
Definition: convolution.hpp:134
slip::separable_convolution3d
void separable_convolution3d(RandomAccessIterator3d1 I_up, RandomAccessIterator3d1 I_bot, KernelIter1 k1_first, KernelIter1 k1_last, std::size_t ksize_left1, std::size_t ksize_right1, slip::BORDER_TREATMENT bt1, KernelIter2 k2_first, KernelIter2 k2_last, std::size_t ksize_left2, std::size_t ksize_right2, slip::BORDER_TREATMENT bt2, KernelIter3 k3_first, KernelIter3 k3_last, std::size_t ksize_left3, std::size_t ksize_right3, slip::BORDER_TREATMENT bt3, RandomAccessIterator3d2 R_up)
Computes the separable convolution of a 3d signal by 3 asymmetric 1d-kernels.
Definition: convolution.hpp:2596
slip::Array2d::cols
size_type cols() const
Returns the number of columns (second dimension size) in the Array2d.
Definition: Array2d.hpp:3155
slip::inner_product
std::iterator_traits< RandomAccessIterator1 >::value_type inner_product(RandomAccessIterator1 first1, RandomAccessIterator1 last1, RandomAccessIterator2 first2, typename std::iterator_traits< RandomAccessIterator1 >::value_type init=typename std::iterator_traits< RandomAccessIterator1 >::value_type())
Computes the inner_product of two ranges X and Y: .
Definition: linear_algebra.hpp:323
slip::full_convolution2d
void full_convolution2d(SrcIter2d S_up, SrcIter2d S_bot, KernelIter2d kernel_up, KernelIter2d kernel_bot, ResIter2d R_up, ResIter2d R_bot)
Computes the full convolution of 2d signal by a 2d-kernel.
Definition: convolution.hpp:2165
slip::Array3d::slice_begin
slice_iterator slice_begin(const size_type row, const size_type col)
Returns a read/write iterator that points to the first element of the line (row,col) threw the slices...
slip::ifft2d
void ifft2d(InputBidirectionalIterator2d in_upper_left, InputBidirectionalIterator2d in_bottom_right, OutputBidirectionalIterator2d out_upper_left)
Computes the ifft2d of a container with 2d iterators.
Definition: FFT.hpp:2629
Array3d.hpp
Provides a class to manipulate 3d dynamic and generic arrays.
slip::Array2d::col_end
col_iterator col_end(const size_type col)
Returns a read/write iterator that points one past the end element of the column column in the Array2...
slip::Array3d::resize
void resize(std::size_t d1, std::size_t d2, std::size_t d3, const T &val=T())
Resizes a Array3d.
Definition: Array3d.hpp:3606
slip::Array2d::row_end
row_iterator row_end(const size_type row)
Returns a read/write iterator that points one past the end element of the row row in the Array2d...
slip::real_fft
void real_fft(InputIter in_begin, InputIter in_end, OutputIter out_begin)
Computes the real fft of a container.
Definition: FFT.hpp:2430
apply.hpp
Provides some algorithms to apply C like functions to ranges.
slip::Signal
Numerical Signal class. This container statisfies the RandomAccessContainer concepts of the Standard ...
Definition: Signal.hpp:105
slip::full_convolution
void full_convolution(SrcIter first, SrcIter last, KernelIter kernel_first, KernelIter kernel_last, ResIter result)
Computes the full convolution of signal by a 1d-kernel.
Definition: convolution.hpp:192
slip::Array2d::begin
iterator begin()
Returns a read/write iterator that points to the first element in the Array2d. Iteration is done in o...
Definition: Array2d.hpp:2345
slip::BORDER_TREATMENT
BORDER_TREATMENT
Choose between different border treatment modes for convolution algorithms.
Definition: border_treatment.hpp:117
slip::Array
This is a linear (one-dimensional) dynamic template container. This container statisfies the RandomAc...
Definition: Array.hpp:94
slip::real_fft3d
void real_fft3d(InputBidirectionalIterator3d in_front_upper_left, InputBidirectionalIterator3d in_back_bottom_right, OutputBidirectionalIterator3d out_front_upper_left)
Computes the real fft3d of a container with 3d iterators.
Definition: FFT.hpp:2804
slip::lin_alg_traits
Definition: linear_algebra_traits.hpp:186
Vector.hpp
Provides a class to manipulate numerical vectors.
slip::Array2d
This is a two-dimensional dynamic and generic container. This container statisfies the Bidirectionnal...
Definition: Array2d.hpp:135
slip::Vector::begin
iterator begin()
Returns a read/write iterator that points to the first element in the Vector. Iteration is done in or...
Definition: Vector.hpp:1474
Generated on Tue Sep 30 2014 10:05:56 for SLIP by   doxygen 1.8.6