Rev Author Line No. Line
178 kaklik 1 #ifndef __MULTI_ARRAY_OP_H
2 #define __MULTI_ARRAY_OP_H
3  
4 #include <boost/array.hpp>
5 #include <boost/multi_array.hpp>
6 #include <functional>
7 #include "functions.h"
8  
9 namespace mimas {
10 /** @defgroup arrayOp Operators for boost::multi_array and mimas::image
11 Mimas has support for doing element-wise operations and convolutions on
12 \c boost::multi_array as well as mimas::image.
13  
14 Implementations for element-wise operations on multi-dimensional arrays
15 are provided. \c boost::multi_array is used to represent the
16 multidimensional arrays.
17  
18 The supported operations are:
19 \li Associative array-scalar operators: multiplication, plus and minus.
20 \li Non-associative array-scalar operators: division, threshold.
21 \li Array-array operators: multiplication, division, plus and minus.
22  
23 The following example demonstrates, how versatile the array-operators
24 provided by mimas are:
25 \include arrayop/main.cc
26  
27  
28 @author Jan Wedekind <jan@wedesoft.de>
29 @author Haifeng Gong <hfgong at users.sourceforge.net>
30 @date Thu Jul 06 14:13:05 UTC 2006
31 @todo Does not compile under Solaris at the moment.
32 @todo Add documentation of all available array operators.
33 @todo Define boost-mpl templates for mimas::rgba to allow the use of
34 expression templates (how?)
35 @{ */
36 ///
37 template<
38 typename T1, typename T2, size_t NumDims, typename Allocator,
39 template< typename, size_t, typename > class MultiArray
40 >
41 boost::multi_array< T1, NumDims, Allocator > empty_clone
42 ( const MultiArray< T2, NumDims, Allocator > &x )
43 {
44 boost::array< size_t, NumDims > shape;
45 std::copy( x.shape(), x.shape() + NumDims, shape.begin() );
46 boost::multi_array< T1, NumDims > retVal( shape );
47 return retVal;
48 }
49  
50 ///
51 template<
52 typename T, class F
53 >
54 boost::detail::multi_array::sub_array< T, 1 > multi_apply
55 ( boost::detail::multi_array::sub_array< T, 1 > a, F f ) {
56 for ( typename boost::detail::multi_array::sub_array< T, 1 >::iterator i = a.begin();
57 i != a.end(); i++ )
58 f( *i );
59 return a;
60 }
61  
62 ///
63 template<
64 typename T, class F,
65 template< typename, size_t > class MultiArray
66 >
67 MultiArray< T, 1 > &multi_apply
68 ( MultiArray< T, 1 > &a, F f ) {
69 for ( typename MultiArray< T, 1 >::iterator i = a.begin();
70 i != a.end(); i++ )
71 f( *i );
72 return a;
73 }
74  
75 #if (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 2))
76 ///
77 template<
78 typename T, class F, typename Allocator,
79 template< typename, size_t, typename > class MultiArray
80 >
81 MultiArray< T, 1, Allocator > &multi_apply
82 ( MultiArray< T, 1, Allocator > &a, F f ) {
83 for ( typename MultiArray< T, 1, Allocator >::iterator i = a.begin();
84 i != a.end(); i++ )
85 f( *i );
86 return a;
87 }
88 #endif
89  
90 ///
91 template<
92 typename T, size_t NumDims, class F
93 >
94 boost::detail::multi_array::sub_array< T, NumDims > multi_apply
95 ( boost::detail::multi_array::sub_array< T, NumDims > a, F f ) {
96 for ( typename boost::detail::multi_array::sub_array< T, NumDims >::iterator i = a.begin();
97 i != a.end(); i++ )
98 multi_apply( *i, f );
99 return a;
100 }
101  
102 ///
103 template<
104 typename T, size_t NumDims, class F,
105 template< typename, size_t > class MultiArray
106 >
107 MultiArray< T, NumDims > &multi_apply
108 ( MultiArray< T, NumDims > &a, F f ) {
109 for ( typename MultiArray< T, NumDims >::iterator i = a.begin();
110 i != a.end(); i++ )
111 multi_apply( *i, f );
112 return a;
113 }
114  
115 #if (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 2))
116 ///
117 template<
118 typename T, size_t NumDims, class F, typename Allocator,
119 template< typename, size_t, typename > class MultiArray
120 >
121 MultiArray< T, NumDims, Allocator > &multi_apply
122 ( MultiArray< T, NumDims, Allocator > &a, F f ) {
123 for ( typename MultiArray< T, NumDims, Allocator >::iterator i = a.begin();
124 i != a.end(); i++ )
125 multi_apply( *i, f );
126 return a;
127 }
128 #endif
129  
130 ///
131 template<
132 typename T1, typename T2, class F, typename Allocator,
133 template< typename, size_t, typename > class MultiArray
134 >
135 boost::detail::multi_array::sub_array< T1, 1 > multi_apply
136 ( boost::detail::multi_array::sub_array< T1, 1 > a,
137 const MultiArray< T2, 1, Allocator > &b,
138 F f ) {
139 typename MultiArray< T2, 1, Allocator >::const_iterator j = b.begin();
140 for ( typename boost::detail::multi_array::sub_array< T1, 1 >::iterator i = a.begin();
141 i != a.end(); i++, j++ )
142 f( *i, *j );
143 return a;
144 }
145  
146 ///
147 template<
148 typename T1, typename T2, class F,
149 typename Allocator2,
150 template< typename, size_t > class MultiArray1,
151 template< typename, size_t, typename > class MultiArray2
152 >
153 MultiArray1< T1, 1 > &multi_apply
154 ( MultiArray1< T1, 1 > &a,
155 const MultiArray2< T2, 1, Allocator2 > &b, F f ) {
156 typename MultiArray2< T2, 1, Allocator2 >::const_iterator j = b.begin();
157 for ( typename MultiArray1< T1, 1 >::iterator i = a.begin();
158 i != a.end(); i++, j++ )
159 f( *i, *j );
160 return a;
161 }
162  
163 #if (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 2))
164 ///
165 template<
166 typename T1, typename T2, class F,
167 typename Allocator1, typename Allocator2,
168 template< typename, size_t, typename > class MultiArray1,
169 template< typename, size_t, typename > class MultiArray2
170 >
171 MultiArray1< T1, 1, Allocator1 > &multi_apply
172 ( MultiArray1< T1, 1, Allocator1 > &a,
173 const MultiArray2< T2, 1, Allocator2 > &b, F f ) {
174 typename MultiArray2< T2, 1, Allocator2 >::const_iterator j = b.begin();
175 for ( typename MultiArray1< T1, 1, Allocator1 >::iterator i = a.begin();
176 i != a.end(); i++, j++ )
177 f( *i, *j );
178 return a;
179 }
180 #endif
181  
182 ///
183 template<
184 typename T1, typename T2, size_t NumDims, class F, typename Allocator,
185 template< typename, size_t, typename > class MultiArray
186 >
187 boost::detail::multi_array::sub_array< T1, NumDims > multi_apply
188 ( boost::detail::multi_array::sub_array< T1, NumDims > a,
189 const MultiArray< T2, NumDims, Allocator > &b,
190 F f ) {
191 typename MultiArray< T2, NumDims, Allocator >::const_iterator j = b.begin();
192 for ( typename boost::detail::multi_array::sub_array< T1, NumDims >::iterator i = a.begin();
193 i != a.end(); i++, j++ )
194 multi_apply( *i, *j, f );
195 return a;
196 }
197  
198 ///
199 template<
200 typename T1, typename T2, size_t NumDims, class F,
201 typename Allocator2,
202 template< typename, size_t > class MultiArray1,
203 template< typename, size_t, typename > class MultiArray2
204 >
205 MultiArray1< T1, NumDims > &multi_apply
206 ( MultiArray1< T1, NumDims > &a,
207 const MultiArray2< T2, NumDims, Allocator2 > &b, F f ) {
208 typename MultiArray2< T2, NumDims, Allocator2 >::const_iterator j = b.begin();
209 for ( typename MultiArray1< T1, NumDims >::iterator i = a.begin();
210 i != a.end(); i++, j++ )
211 multi_apply( *i, *j, f );
212 return a;
213 }
214  
215 #if (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 2))
216 ///
217 template<
218 typename T1, typename T2, size_t NumDims, class F,
219 typename Allocator1, typename Allocator2,
220 template< typename, size_t, typename > class MultiArray1,
221 template< typename, size_t, typename > class MultiArray2
222 >
223 MultiArray1< T1, NumDims, Allocator1 > &multi_apply
224 ( MultiArray1< T1, NumDims, Allocator1 > &a,
225 const MultiArray2< T2, NumDims, Allocator2 > &b, F f ) {
226 typename MultiArray2< T2, NumDims, Allocator2 >::const_iterator j = b.begin();
227 for ( typename MultiArray1< T1, NumDims, Allocator1 >::iterator i = a.begin();
228 i != a.end(); i++, j++ )
229 multi_apply( *i, *j, f );
230 return a;
231 }
232 #endif
233  
234 ///
235 template<
236 typename T1, typename T2, typename T3, class F,
237 typename Allocator2, typename Allocator3,
238 template< typename, size_t, typename > class MultiArray2,
239 template< typename, size_t, typename > class MultiArray3
240 >
241 boost::detail::multi_array::sub_array< T1, 1 > multi_apply
242 ( boost::detail::multi_array::sub_array< T1, 1 > a,
243 const MultiArray2< T2, 1, Allocator2 > &b,
244 const MultiArray3< T3, 1, Allocator3 > &c,
245 F f ) {
246 typename MultiArray2< T2, 1, Allocator2 >::const_iterator j = b.begin();
247 typename MultiArray3< T3, 1, Allocator3 >::const_iterator k = c.begin();
248 for ( typename boost::detail::multi_array::sub_array< T1, 1 >::iterator i = a.begin();
249 i != a.end(); i++, j++, k++ )
250 f( *i, *j, *k );
251 return a;
252 }
253  
254 ///
255 template<
256 typename T1, typename T2, typename T3, class F,
257 typename Allocator2, typename Allocator3,
258 template< typename, size_t > class MultiArray1,
259 template< typename, size_t, typename > class MultiArray2,
260 template< typename, size_t, typename > class MultiArray3
261 >
262 MultiArray1< T1, 1 > &multi_apply
263 ( MultiArray1< T1, 1 > &a,
264 const MultiArray2< T2, 1, Allocator2 > &b,
265 const MultiArray3< T3, 1, Allocator3 > &c, F f ) {
266 typename MultiArray2< T2, 1, Allocator2 >::const_iterator j = b.begin();
267 typename MultiArray3< T3, 1, Allocator3 >::const_iterator k = c.begin();
268 for ( typename MultiArray1< T1, 1 >::iterator i = a.begin();
269 i != a.end(); i++, j++, k++ )
270 f( *i, *j, *k );
271 return a;
272 }
273  
274 #if (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 2))
275 ///
276 template<
277 typename T1, typename T2, typename T3, class F,
278 typename Allocator1, typename Allocator2, typename Allocator3,
279 template< typename, size_t, typename > class MultiArray1,
280 template< typename, size_t, typename > class MultiArray2,
281 template< typename, size_t, typename > class MultiArray3
282 >
283 MultiArray1< T1, 1, Allocator1 > &multi_apply
284 ( MultiArray1< T1, 1, Allocator1 > &a,
285 const MultiArray2< T2, 1, Allocator2 > &b,
286 const MultiArray3< T3, 1, Allocator3 > &c, F f ) {
287 typename MultiArray2< T2, 1, Allocator2 >::const_iterator j = b.begin();
288 typename MultiArray3< T3, 1, Allocator3 >::const_iterator k = c.begin();
289 for ( typename MultiArray1< T1, 1, Allocator1 >::iterator i = a.begin();
290 i != a.end(); i++, j++, k++ )
291 f( *i, *j, *k );
292 return a;
293 }
294 #endif
295  
296 ///
297 template<
298 typename T1, typename T2, typename T3, size_t NumDims, class F,
299 typename Allocator2, typename Allocator3,
300 template< typename, size_t, typename > class MultiArray2,
301 template< typename, size_t, typename > class MultiArray3
302 >
303 boost::detail::multi_array::sub_array< T1, NumDims > multi_apply
304 ( boost::detail::multi_array::sub_array< T1, NumDims > a,
305 const MultiArray2< T2, NumDims, Allocator2 > &b,
306 const MultiArray3< T3, NumDims, Allocator3 > &c,
307 F f ) {
308 typename MultiArray2< T2, NumDims, Allocator2 >::const_iterator j = b.begin();
309 typename MultiArray3< T3, NumDims, Allocator3 >::const_iterator k = c.begin();
310 for ( typename boost::detail::multi_array::sub_array< T1, NumDims >::iterator i = a.begin();
311 i != a.end(); i++, j++, k++ ) multi_apply( *i, *j, *k, f );
312 return a;
313 }
314  
315 ///
316 template<
317 typename T1, typename T2, typename T3, size_t NumDims, class F,
318 typename Allocator2, typename Allocator3,
319 template< typename, size_t > class MultiArray1,
320 template< typename, size_t, typename > class MultiArray2,
321 template< typename, size_t, typename > class MultiArray3
322 >
323 MultiArray1< T1, NumDims > &multi_apply
324 ( MultiArray1< T1, NumDims > &a,
325 const MultiArray2< T2, NumDims, Allocator2 > &b,
326 const MultiArray3< T3, NumDims, Allocator3 > &c, F f ) {
327 typename MultiArray2< T2, NumDims, Allocator2 >::const_iterator j = b.begin();
328 typename MultiArray3< T3, NumDims, Allocator3 >::const_iterator k = c.begin();
329 for ( typename MultiArray1< T1, NumDims >::iterator i = a.begin();
330 i != a.end(); i++, j++, k++ )
331 multi_apply( *i, *j, *k, f );
332 return a;
333 }
334  
335 #if (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 2))
336 ///
337 template<
338 typename T1, typename T2, typename T3, size_t NumDims, class F,
339 typename Allocator1, typename Allocator2, typename Allocator3,
340 template< typename, size_t, typename > class MultiArray1,
341 template< typename, size_t, typename > class MultiArray2,
342 template< typename, size_t, typename > class MultiArray3
343 >
344 MultiArray1< T1, NumDims, Allocator1 > &multi_apply
345 ( MultiArray1< T1, NumDims, Allocator1 > &a,
346 const MultiArray2< T2, NumDims, Allocator2 > &b,
347 const MultiArray3< T3, NumDims, Allocator3 > &c, F f ) {
348 typename MultiArray2< T2, NumDims, Allocator2 >::const_iterator j = b.begin();
349 typename MultiArray3< T3, NumDims, Allocator3 >::const_iterator k = c.begin();
350 for ( typename MultiArray1< T1, NumDims, Allocator1 >::iterator i = a.begin();
351 i != a.end(); i++, j++, k++ )
352 multi_apply( *i, *j, *k, f );
353 return a;
354 }
355 #endif
356  
357 ///
358 template< typename T1, typename T2, class F >
359 struct _multi_help1
360 {
361 _multi_help1( F _f ): f(_f) {}
362 T1 &operator()( T1 &x, const T2 &y ) const
363 { x = f( y ); return x; }
364 F f;
365 };
366  
367 ///
368 template< typename T1, typename T2, typename T3, class F >
369 struct _multi_help2
370 {
371 _multi_help2( F _f ): f(_f) {}
372 T1 &operator()( T1 &x, const T2 &y, const T3 &z ) const
373 { x = f( y, z ); return x; }
374 F f;
375 };
376  
377 ///
378 template<
379 typename T1, typename T2, size_t NumDims, class F, typename Allocator,
380 template< typename, size_t, typename > class MultiArray
381 >
382 boost::multi_array< T1, NumDims, Allocator > multi_func
383 ( const MultiArray< T2, NumDims, Allocator > &a, F f ) {
384 boost::multi_array< T1, NumDims > retVal( empty_clone< T1 >( a ) );
385 return multi_apply( retVal, a, _multi_help1< T1, T2, F >( f ) );
386 }
387  
388 ///
389 template<
390 typename T1, typename T2, typename T3, size_t NumDims, class F,
391 typename Allocator1, typename Allocator2,
392 template< typename, size_t, typename > class MultiArray1,
393 template< typename, size_t, typename > class MultiArray2
394 >
395 boost::multi_array< T1, NumDims > multi_func
396 ( const MultiArray1< T2, NumDims, Allocator1 > &a,
397 const MultiArray2< T3, NumDims, Allocator2 > &b, F f ) {
398 boost::multi_array< T1, NumDims > retVal( empty_clone< T1 >( a ) );
399 return multi_apply( retVal, a, b, _multi_help2< T1, T2, T3, F >( f ) );
400 }
401  
402 ///
403 template<
404 typename T1, typename T2, size_t NumDims, typename Allocator,
405 template< typename, size_t, typename > class MultiArray
406 >
407 boost::multi_array< T1, NumDims, Allocator > multi_cast
408 ( const MultiArray< T2, NumDims, Allocator > &a ) {
409 return multi_func< T1 >( a, std::_Identity< T2 >() );
410 }
411  
412 ///@}
413  
414 }
415  
416 #define __MIMASINTERNALARRAYFUNC operator*=
417 #define __MIMASEXTERNALARRAYFUNC operator*
418 #define __MIMASFUNCTIONOBJECT std::multiplies
419 #include "multi_array_op_help.h"
420  
421 #define __MIMASINTERNALARRAYFUNC operator/=
422 #define __MIMASEXTERNALARRAYFUNC operator/
423 #define __MIMASFUNCTIONOBJECT std::divides
424 #include "multi_array_op_help.h"
425  
426 #define __MIMASINTERNALARRAYFUNC operator+=
427 #define __MIMASEXTERNALARRAYFUNC operator+
428 #define __MIMASFUNCTIONOBJECT std::plus
429 #include "multi_array_op_help.h"
430  
431 #define __MIMASINTERNALARRAYFUNC operator-=
432 #define __MIMASEXTERNALARRAYFUNC operator-
433 #define __MIMASFUNCTIONOBJECT std::minus
434 #include "multi_array_op_help.h"
435  
436 #define __MIMASEXTERNALARRAYFUNC absolute
437 #define __MIMASINTERNALARRAYFUNC absoluteIt
438 #define __MIMASFUNCTIONOBJECT _abs
439 #include "multi_array_op_help2.h"
440  
441 #define __MIMASEXTERNALARRAYFUNC conj
442 #define __MIMASINTERNALARRAYFUNC conjIt
443 #define __MIMASFUNCTIONOBJECT _conj
444 #include "multi_array_op_help2.h"
445  
446 #define __MIMASEXTERNALARRAYFUNC sqr
447 #define __MIMASINTERNALARRAYFUNC sqrIt
448 #define __MIMASFUNCTIONOBJECT _sqr
449 #include "multi_array_op_help2.h"
450  
451 #define __MIMASEXTERNALARRAYFUNC logarithm
452 #define __MIMASINTERNALARRAYFUNC logarithmIt
453 #define __MIMASFUNCTIONOBJECT _log
454 #include "multi_array_op_help2.h"
455  
456 #define __MIMASEXTERNALARRAYFUNC squareRoot
457 #define __MIMASINTERNALARRAYFUNC squareRootIt
458 #define __MIMASFUNCTIONOBJECT _sqrt
459 #include "multi_array_op_help2.h"
460  
461 #define __MIMASEXTERNALARRAYFUNC sumSquares
462 #define __MIMASFUNCTIONOBJECT _sumsquares
463 #include "multi_array_op_help3.h"
464  
465 #define __MIMASEXTERNALARRAYFUNC orientation
466 #define __MIMASFUNCTIONOBJECT _orientation
467 #include "multi_array_op_help3.h"
468  
469 namespace mimas {
470 /** @addtogroup arrayOp
471 @{ */
472 ///
473 template <
474 typename T1, typename T2, size_t NumDims, typename Allocator,
475 template< typename, size_t, typename > class MultiArray
476 >
477 boost::multi_array< T1, NumDims, Allocator > norm( const MultiArray< T2, NumDims, Allocator > &a )
478 {
479 return multi_func< T1 >( a, _norm< T1, T2 >() );
480 }
481  
482 ///
483 template <
484 typename T1, typename T2, size_t NumDims, typename Allocator,
485 template< typename, size_t, typename > class MultiArray
486 >
487 boost::multi_array< T1, NumDims, Allocator > arg( const MultiArray< T2, NumDims, Allocator > &a )
488 {
489 return multi_func< T1 >( a, _arg< T1, T2 >() );
490 }
491  
492 ///
493 template <
494 typename T, size_t NumDims, typename Allocator,
495 template< typename, size_t, typename > class MultiArray
496 >
497 boost::multi_array< int, NumDims, Allocator > fastSqr( const MultiArray< T, NumDims, Allocator > &a )
498 {
499 return multi_func< int >( a, _fastsqr< T >() );
500 }
501  
502 ///@}
503 }
504  
505 #endif