GNU Octave  4.0.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
oct-convn.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2010-2015 VZLU Prague
4 
5 This file is part of Octave.
6 
7 Octave is free software; you can redistribute it and/or modify it
8 under the terms of the GNU General Public License as published by the
9 Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
11 
12 Octave is distributed in the hope that it will be useful, but WITHOUT
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15 for more details.
16 
17 You should have received a copy of the GNU General Public License
18 along with Octave; see the file COPYING. If not, see
19 <http://www.gnu.org/licenses/>.
20 
21 */
22 
23 #ifdef HAVE_CONFIG_H
24 #include <config.h>
25 #endif
26 
27 #include <iostream>
28 #include <algorithm>
29 
30 #include "f77-fcn.h"
31 
32 #include "oct-convn.h"
33 #include "oct-locbuf.h"
34 
35 // 2d convolution with a matrix kernel.
36 template <class T, class R>
37 static void
38 convolve_2d (const T *a, octave_idx_type ma, octave_idx_type na,
39  const R *b, octave_idx_type mb, octave_idx_type nb,
40  T *c, bool inner);
41 
42 // Forward instances to our Fortran implementations.
43 #define FORWARD_IMPL(T,R,f,F) \
44 extern "C" \
45 F77_RET_T \
46 F77_FUNC (f##conv2o, F##CONV2O) (const octave_idx_type&, \
47  const octave_idx_type&, \
48  const T*, const octave_idx_type&, \
49  const octave_idx_type&, const R*, T *); \
50 \
51 extern "C" \
52 F77_RET_T \
53 F77_FUNC (f##conv2i, F##CONV2I) (const octave_idx_type&, \
54  const octave_idx_type&, \
55  const T*, const octave_idx_type&, \
56  const octave_idx_type&, const R*, T *); \
57 \
58 template <> void \
59 convolve_2d<T, R> (const T *a, octave_idx_type ma, octave_idx_type na, \
60  const R *b, octave_idx_type mb, octave_idx_type nb, \
61  T *c, bool inner) \
62 { \
63  if (inner) \
64  F77_XFCN (f##conv2i, F##CONV2I, (ma, na, a, mb, nb, b, c)); \
65  else \
66  F77_XFCN (f##conv2o, F##CONV2O, (ma, na, a, mb, nb, b, c)); \
67 }
68 
69 FORWARD_IMPL (double, double, d, D)
70 FORWARD_IMPL (float, float, s, S)
71 FORWARD_IMPL (Complex, Complex, z, Z)
72 FORWARD_IMPL (FloatComplex, FloatComplex, c, C)
73 FORWARD_IMPL (Complex, double, zd, ZD)
74 FORWARD_IMPL (FloatComplex, float, cs, CS)
75 
76 template <class T, class R>
77 void convolve_nd (const T *a, const dim_vector& ad, const dim_vector& acd,
78  const R *b, const dim_vector& bd, const dim_vector& bcd,
79  T *c, const dim_vector& ccd, int nd, bool inner)
80 {
81  if (nd == 2)
82  convolve_2d<T, R> (a, ad(0), ad(1), b, bd(0), bd(1), c, inner);
83  else
84  {
85  octave_idx_type ma = acd(nd-2);
86  octave_idx_type na = ad(nd-1);
87  octave_idx_type mb = bcd(nd-2);
88  octave_idx_type nb = bd(nd-1);
89  octave_idx_type ldc = ccd(nd-2);
90  if (inner)
91  {
92  for (octave_idx_type ja = 0; ja < na - nb + 1; ja++)
93  for (octave_idx_type jb = 0; jb < nb; jb++)
94  convolve_nd<T, R> (a + ma*(ja+jb), ad, acd,
95  b + mb*(nb-jb-1), bd, bcd,
96  c + ldc*ja, ccd, nd-1, inner);
97  }
98  else
99  {
100  for (octave_idx_type ja = 0; ja < na; ja++)
101  for (octave_idx_type jb = 0; jb < nb; jb++)
102  convolve_nd<T, R> (a + ma*ja, ad, acd, b + mb*jb, bd, bcd,
103  c + ldc*(ja+jb), ccd, nd-1, inner);
104  }
105  }
106 }
107 
108 // Arbitrary convolutor.
109 // The 2nd array is assumed to be the smaller one.
110 template <class T, class R>
111 static MArray<T>
112 convolve (const MArray<T>& a, const MArray<R>& b,
113  convn_type ct)
114 {
115  if (a.is_empty () || b.is_empty ())
116  return MArray<T> ();
117 
118  int nd = std::max (a.ndims (), b.ndims ());
119  const dim_vector adims = a.dims ().redim (nd);
120  const dim_vector bdims = b.dims ().redim (nd);
121  dim_vector cdims = dim_vector::alloc (nd);
122 
123  for (int i = 0; i < nd; i++)
124  {
125  if (ct == convn_valid)
126  cdims(i) = std::max (adims(i) - bdims(i) + 1,
127  static_cast<octave_idx_type> (0));
128  else
129  cdims(i) = std::max (adims(i) + bdims(i) - 1,
130  static_cast<octave_idx_type> (0));
131  }
132 
133  MArray<T> c (cdims, T ());
134 
135  convolve_nd<T, R> (a.fortran_vec (), adims, adims.cumulative (),
136  b.fortran_vec (), bdims, bdims.cumulative (),
137  c.fortran_vec (), cdims.cumulative (),
138  nd, ct == convn_valid);
139 
140  if (ct == convn_same)
141  {
142  // Pick the relevant part.
143  Array<idx_vector> sidx (dim_vector (nd, 1));
144 
145  for (int i = 0; i < nd; i++)
146  sidx(i) = idx_vector::make_range (bdims(i)/2, 1, adims(i));
147  c = c.index (sidx);
148  }
149 
150  return c;
151 }
152 
153 #define CONV_DEFS(TPREF, RPREF) \
154 TPREF ## NDArray \
155 convn (const TPREF ## NDArray& a, const RPREF ## NDArray& b, convn_type ct) \
156 { \
157  return convolve (a, b, ct); \
158 } \
159 TPREF ## Matrix \
160 convn (const TPREF ## Matrix& a, const RPREF ## Matrix& b, convn_type ct) \
161 { \
162  return convolve (a, b, ct); \
163 } \
164 TPREF ## Matrix \
165 convn (const TPREF ## Matrix& a, const RPREF ## ColumnVector& c, \
166  const RPREF ## RowVector& r, convn_type ct) \
167 { \
168  return convolve (a, c * r, ct); \
169 }
170 
173 CONV_DEFS (Complex, Complex)
174 CONV_DEFS (Float, Float)
176 CONV_DEFS (FloatComplex, FloatComplex)
bool is_empty(void) const
Definition: Array.h:472
#define C(a, b)
Definition: Faddeeva.cc:255
static MArray< T > convolve(const MArray< T > &a, const MArray< R > &b, convn_type ct)
Definition: oct-convn.cc:112
convn_type
Definition: oct-convn.h:47
int ndims(void) const
Definition: Array.h:487
dim_vector cumulative(void) const
Return cumulative dimensions.
Definition: dim-vector.h:489
#define FORWARD_IMPL(T, R, f, F)
Definition: oct-convn.cc:43
Definition: MArray.h:36
F77_RET_T const octave_idx_type const octave_idx_type const octave_idx_type double const octave_idx_type double const octave_idx_type double const octave_idx_type double * Z
Definition: qz.cc:114
F77_RET_T const double const double double * d
const dim_vector & dims(void) const
Return a const-reference so that dims ()(i) works efficiently.
Definition: Array.h:337
void convolve_nd(const T *a, const dim_vector &ad, const dim_vector &acd, const R *b, const dim_vector &bd, const dim_vector &bcd, T *c, const dim_vector &ccd, int nd, bool inner)
Definition: oct-convn.cc:77
#define CONV_DEFS(TPREF, RPREF)
Definition: oct-convn.cc:153
static void convolve_2d(const T *a, octave_idx_type ma, octave_idx_type na, const R *b, octave_idx_type mb, octave_idx_type nb, T *c, bool inner)
dim_vector redim(int n) const
Definition: dim-vector.cc:266
static dim_vector alloc(int n)
Definition: dim-vector.h:256
Handles the reference counting for all the derived classes.
Definition: Array.h:45
charNDArray max(char d, const charNDArray &m)
Definition: chNDArray.cc:233
std::complex< float > FloatComplex
Definition: oct-cmplx.h:30
std::complex< double > Complex
Definition: oct-cmplx.h:29
const T * fortran_vec(void) const
Definition: Array.h:481
static idx_vector make_range(octave_idx_type start, octave_idx_type step, octave_idx_type len)
Definition: idx-vector.h:476
Array< T > index(const idx_vector &i) const
Indexing without resizing.
Definition: Array.cc:716