GNU Octave  4.2.1
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-2017 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 #if defined (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 <typename T, typename R>
37 static void
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_CXX, R_CXX, T, R, T_CAST, T_CONST_CAST, \
44  R_CONST_CAST, f, F) \
45  extern "C" \
46  F77_RET_T \
47  F77_FUNC (f##conv2o, F##CONV2O) (const F77_INT&, \
48  const F77_INT&, \
49  const T*, const F77_INT&, \
50  const F77_INT&, const R*, T *); \
51  \
52  extern "C" \
53  F77_RET_T \
54  F77_FUNC (f##conv2i, F##CONV2I) (const F77_INT&, \
55  const F77_INT&, \
56  const T*, const F77_INT&, \
57  const F77_INT&, const R*, T *); \
58  \
59  template <> void \
60  convolve_2d<T_CXX, R_CXX> (const T_CXX *a, F77_INT ma, F77_INT na, \
61  const R_CXX *b, F77_INT mb, F77_INT nb, \
62  T_CXX *c, bool inner) \
63  { \
64  if (inner) \
65  F77_XFCN (f##conv2i, F##CONV2I, (ma, na, T_CONST_CAST (a), \
66  mb, nb, R_CONST_CAST (b), \
67  T_CAST (c))); \
68  else \
69  F77_XFCN (f##conv2o, F##CONV2O, (ma, na, T_CONST_CAST (a), \
70  mb, nb, R_CONST_CAST (b), \
71  T_CAST (c))); \
72  }
73 
74 FORWARD_IMPL (double, double, F77_DBLE, F77_DBLE, , , , d, D)
75 FORWARD_IMPL (float, float, F77_REAL, F77_REAL, , , , s, S)
76 
77 FORWARD_IMPL (std::complex<double>, std::complex<double>,
78  F77_DBLE_CMPLX, F77_DBLE_CMPLX, F77_DBLE_CMPLX_ARG,
79  F77_CONST_DBLE_CMPLX_ARG, F77_CONST_DBLE_CMPLX_ARG, z, Z)
80 FORWARD_IMPL (std::complex<float>, std::complex<float>,
81  F77_CMPLX, F77_CMPLX, F77_CMPLX_ARG,
82  F77_CONST_CMPLX_ARG, F77_CONST_CMPLX_ARG, c, C)
83 
84 FORWARD_IMPL (std::complex<double>, double,
85  F77_DBLE_CMPLX, F77_DBLE, F77_DBLE_CMPLX_ARG,
86  F77_CONST_DBLE_CMPLX_ARG, , zd, ZD)
87 FORWARD_IMPL (std::complex<float>, float, F77_CMPLX, F77_REAL, F77_CMPLX_ARG,
88  F77_CONST_CMPLX_ARG, , cs, CS)
89 
90 template <typename T, typename R>
91 void convolve_nd (const T *a, const dim_vector& ad, const dim_vector& acd,
92  const R *b, const dim_vector& bd, const dim_vector& bcd,
93  T *c, const dim_vector& ccd, int nd, bool inner)
94 {
95  if (nd == 2)
96  convolve_2d<T, R> (a, ad(0), ad(1), b, bd(0), bd(1), c, inner);
97  else
98  {
99  octave_idx_type ma = acd(nd-2);
100  octave_idx_type na = ad(nd-1);
101  octave_idx_type mb = bcd(nd-2);
102  octave_idx_type nb = bd(nd-1);
103  octave_idx_type ldc = ccd(nd-2);
104  if (inner)
105  {
106  for (octave_idx_type ja = 0; ja < na - nb + 1; ja++)
107  for (octave_idx_type jb = 0; jb < nb; jb++)
108  convolve_nd<T, R> (a + ma*(ja+jb), ad, acd,
109  b + mb*(nb-jb-1), bd, bcd,
110  c + ldc*ja, ccd, nd-1, inner);
111  }
112  else
113  {
114  for (octave_idx_type ja = 0; ja < na; ja++)
115  for (octave_idx_type jb = 0; jb < nb; jb++)
116  convolve_nd<T, R> (a + ma*ja, ad, acd, b + mb*jb, bd, bcd,
117  c + ldc*(ja+jb), ccd, nd-1, inner);
118  }
119  }
120 }
121 
122 // Arbitrary convolutor.
123 // The 2nd array is assumed to be the smaller one.
124 template <typename T, typename R>
125 static MArray<T>
126 convolve (const MArray<T>& a, const MArray<R>& b,
127  convn_type ct)
128 {
129  if (a.is_empty () || b.is_empty ())
130  return MArray<T> ();
131 
132  int nd = std::max (a.ndims (), b.ndims ());
133  const dim_vector adims = a.dims ().redim (nd);
134  const dim_vector bdims = b.dims ().redim (nd);
135  dim_vector cdims = dim_vector::alloc (nd);
136 
137  for (int i = 0; i < nd; i++)
138  {
139  if (ct == convn_valid)
140  cdims(i) = std::max (adims(i) - bdims(i) + 1,
141  static_cast<octave_idx_type> (0));
142  else
143  cdims(i) = std::max (adims(i) + bdims(i) - 1,
144  static_cast<octave_idx_type> (0));
145  }
146 
147  MArray<T> c (cdims, T ());
148 
149  convolve_nd<T, R> (a.fortran_vec (), adims, adims.cumulative (),
150  b.fortran_vec (), bdims, bdims.cumulative (),
151  c.fortran_vec (), cdims.cumulative (),
152  nd, ct == convn_valid);
153 
154  if (ct == convn_same)
155  {
156  // Pick the relevant part.
157  Array<idx_vector> sidx (dim_vector (nd, 1));
158 
159  for (int i = 0; i < nd; i++)
160  sidx(i) = idx_vector::make_range (bdims(i)/2, 1, adims(i));
161  c = c.index (sidx);
162  }
163 
164  return c;
165 }
166 
167 #define CONV_DEFS(TPREF, RPREF) \
168  TPREF ## NDArray \
169  convn (const TPREF ## NDArray& a, const RPREF ## NDArray& b, \
170  convn_type ct) \
171  { \
172  return convolve (a, b, ct); \
173  } \
174  TPREF ## Matrix \
175  convn (const TPREF ## Matrix& a, const RPREF ## Matrix& b, \
176  convn_type ct) \
177  { \
178  return convolve (a, b, ct); \
179  } \
180  TPREF ## Matrix \
181  convn (const TPREF ## Matrix& a, const RPREF ## ColumnVector& c, \
182  const RPREF ## RowVector& r, convn_type ct) \
183  { \
184  return convolve (a, c * r, ct); \
185  }
186 
189 CONV_DEFS (Complex, Complex)
190 CONV_DEFS (Float, Float)
192 CONV_DEFS (FloatComplex, FloatComplex)
bool is_empty(void) const
Definition: Array.h:575
#define C(a, b)
Definition: Faddeeva.cc:246
static MArray< T > convolve(const MArray< T > &a, const MArray< R > &b, convn_type ct)
Definition: oct-convn.cc:126
convn_type
Definition: oct-convn.h:49
int ndims(void) const
Definition: Array.h:590
dim_vector cumulative(void) const
Return cumulative dimensions.
Definition: dim-vector.h:509
#define F77_DBLE_CMPLX_ARG(x)
Definition: f77-fcn.h:345
STL namespace.
Template for N-dimensional array classes with like-type math operators.
Definition: MArray.h:32
s
Definition: file-io.cc:2682
F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T F77_REAL F77_REAL &F77_RET_T F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE * d
calling an anonymous function involves an overhead quite comparable to the overhead of an m file function Passing a handle to a built in function is because the interpreter is not involved in the internal loop For a
Definition: cellfun.cc:398
#define F77_REAL
Definition: f77-fcn.h:332
const dim_vector & dims(void) const
Return a const-reference so that dims ()(i) works efficiently.
Definition: Array.h:439
#define FORWARD_IMPL(T_CXX, R_CXX, T, R, T_CAST, T_CONST_CAST,R_CONST_CAST, f, F)
Definition: oct-convn.cc:43
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:91
#define F77_CMPLX
Definition: f77-fcn.h:334
#define CONV_DEFS(TPREF, RPREF)
Definition: oct-convn.cc:167
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:275
F77_RET_T const F77_INT const F77_INT const F77_INT F77_DBLE const F77_INT F77_DBLE const F77_INT F77_DBLE const F77_INT F77_DBLE * Z
the sparsity preserving column transformation such that that defines the pivoting threshold can be given in which case it defines the c
Definition: lu.cc:138
static dim_vector alloc(int n)
Definition: dim-vector.h:270
N Dimensional Array with copy-on-write semantics.
Definition: Array.h:126
charNDArray max(char d, const charNDArray &m)
Definition: chNDArray.cc:228
#define F77_CMPLX_ARG(x)
Definition: f77-fcn.h:339
#define F77_CONST_DBLE_CMPLX_ARG(x)
Definition: f77-fcn.h:348
=val(i)}if ode{val(i)}occurs in table i
Definition: lookup.cc:239
b
Definition: cellfun.cc:398
#define F77_DBLE_CMPLX
Definition: f77-fcn.h:333
#define F77_CONST_CMPLX_ARG(x)
Definition: f77-fcn.h:342
std::complex< float > FloatComplex
Definition: oct-cmplx.h:32
std::complex< double > Complex
Definition: oct-cmplx.h:31
const T * fortran_vec(void) const
Definition: Array.h:584
Vector representing the dimensions (size) of an Array.
Definition: dim-vector.h:87
static idx_vector make_range(octave_idx_type start, octave_idx_type step, octave_idx_type len)
Definition: idx-vector.h:466
#define F77_DBLE
Definition: f77-fcn.h:331
Array< T > index(const idx_vector &i) const
Indexing without resizing.
Definition: Array.cc:718